Quantum-classical crossover in the spin-1/2 Heisenberg-Kitaev kagome magnet

The spin-1/2 Heisenberg kagome antiferromagnet is one of the paradigmatic playgrounds for frustrated quantum magnetism, with an extensive number of competing resonating valence bond (RVB) states emerging at low energies, including gapped and gapless spin liquids and valence bond crystals. Here we revisit the crossover from this quantum RVB phase to a semiclassical regime brought about by anisotropic Kitaev interactions, and focus on the precise mechanisms underpinning this crossover. To this end, we introduce a simple parametrization of the classical ground states (GSs) in terms of emergent Ising-like variables, and use this parametrizaton: i) to construct an effective low-energy description of the order-by-disorder mechanism operating in a large part of the phase diagram, and ii) to contrast, side by side, exact diagonalization data obtained from the full basis with that obtained from the restricted (orthonormalized) basis of classical GSs. The results reveal that fluctuation corrections from states outside the restricted basis are strongly quenched inside the semiclassical regime (due to the large anisotropy spin gaps), and that the RVB phase survives up to a relatively large value of Kitaev anisotropy $K$. We further find that the pure Kitaev model admits a subextensive number of one-dimensional symmetries, which explains naturally the absence of classical and quantum order by disorder reported previously.

Here we revisit the spin-1/2 Heisenberg-Kitaev (or JK-) model on the kagome lattice, which may be relevant for some rare-earth based compounds of the type A 2 RE 3 Sb 3 O 14 , where A=Mg or Zn and RE is a rare-earth ion [34][35][36][37][38].This model interpolates between the spin-1/2 kagome Heisenberg antiferromagnet (KHAF) -one of the paradigmatic playgrounds of competing resonating valence bond (RVB) states, including valence bond crystals as well as gapped and gapless spin liquids [11, -and the compass-like, Kitaev model in which the coupling in spin space is tied to the orientation of the bonds [60].
The JK model on the kagome lattice has been studied previously by Kimchi and Vishwanath [22] and by Morita et al [31,32], and a lot of results are already known, including most of the aspects of the classical ground state (GS) phase diagram [22,31], a numerical demonstration of an order by disorder mechanism operating in a large part of the parameter space, and the absence of this mechanism in the pure Kitaev model [31].However, the microscopic origin of the order by disorder mechanism and its absence in the pure Kitaev model has not been understood.More importantly, the question of whether the RVB phase of the KHAF remains robust in an extended parameter space (and, in particular, whether the answer depends on the nature of the ground state of the KHAF or the presence of a gap) has not been settled.
The main results from our study can be summarized as follows.First, we provide a simple parametrization of the classical GS manifold in terms of emergent Ising-like variables.This parametrization offers a convenient platform for the analysis of the quantum model.Second, we show that the pure Kitaev model admits a subextensive number of one-dimensional symmetries, which explain naturally the reported absence of classical and quantum order by disorder [31,32].Third, we perform a semiclassical perturbative expansion and derive an effective Hamiltonian in terms of the emergent Ising variables.This Hamiltonian provides a simple picture for the order by disorder effect observed numerically by Morita et al [31].Fourth, we demonstrate explicitly the quantum-classical crossover between the RVB physics of the KHAF and the regime stabilized by Kitaev anisotropy.This is achieved by contrasting, side by side, ED results in the full basis and in the restricted orthonormalized basis of classical GSs.The comparison shows that the regime stabilized by the Kitaev coupling has a robust semiclassical character, meaning that the fluctuation (e.g., spin-wave) corrections from states outside the restricted basis are heavily quenched by the large anisotropy spin gaps.In turn, this shows that the RVB phase remains stable in an extended range of parameters, irrespective of the actual nature of the ground state of the KHAF or the presence of a gap.In particular, our finite-size results suggest that the RVB range can extend up to relatively large values of |K| ∼ J, where J is the Heisenberg coupling.
The remaining part of the paper is organized as follows.We begin in Sec.II with a general discussion of the model, its symmetries and duality transformations.In Sec.III we revisit the classical phase diagram using the Luttinger-Tisza approach [61][62][63].This approach leads naturally to a simple parametrization of the classical GS manifold in terms of Ising-like variables.In Sec.IV we present a semiclassical perturbative expansion that reveals the order by disorder mechanism operating in a large part of the phase diagram.In Sec.V we present our extensive numerical study of the quantum spin S = 1/2 model, obtained from exact diagonalizations in two different bases, one in the full basis and the other in the restricted orthonormalized basis of classical GSs.In Sec.VI we provide a discussion and a broader perspective of our study.Technical details and auxiliary information are relegated to three Appendices (App.A-C).

A. Model
We consider interacting spins S i residing at the vertices i of the 2D kagome lattice, a portion of which is shown in Fig. 1.The kagome has a triangular Bravais lattice and a basis of three sites, A, B and C (shaded triangles).The nearest neighbour (NN) bonds of the lattice can be divided into three types, 'xx' (red), 'yy' (green) and 'zz' (blue), depending on their orientation, see Fig. 1.The Heisenberg-Kitaev or JKmodel is described by the spin Hamiltonian Here i j denotes NN lattice sites, S i and S j are the associated pseudospins-1/2 degrees of freedom residing on these sites, and J and K denote the Heisenberg and Kitaev exchange couplings, respectively.The Cartesian components γ i j appearing in the Kitaev coupling equals x, y or z, depending on whether i j belongs to the 'xx', 'yy' or 'zz' bond type.In the following, we measure energy in units of J 2 + K 2 = 1 and parametrize

B. Global symmetries
For half-integer spins, the Hamiltonian (1) is in general invariant under the global symmetry group T × C 3v × D 2 , which consists of the following operations: i) The translation group T generated by the primitive translation vectors a 1 and a 2 shown in Fig. 1.
ii) The double cover C 3v of the group C 3v ⊂ SO(3) in spinorbit space, with the three-fold axis going through one of hexagon centers (see Fig. 1).This is a three-fold rotation around [111], which maps 'xx', 'yy' and 'zz' bonds into 'yy', 'zz' and 'xx' bonds in real space, and (S x , S y , S z ) → (S y , S z , S x ) in spin space.The reflection planes (1 10), (01 1) and ( 101), are shown by dashed (brown) lines in Fig. 1.In spin space, a reflection through (1 10) maps (S x , S y , S z ) → (−S y , −S x , −S z ), and similarly for the other planes.
iii) The double cover D 2 of the point group D 2 ⊂ SO(3) in spin space alone.This consists of the three π-rotations C 2x , C 2y , and C 2z , which map (S x , S y , S z ) to (S x , −S y , −S z ), (−S x , S y , −S z ) and (−S x , −S y , S z ), respectively.FIG. 1.The kagome lattice has a basis of three sites (A, B and C) and a triangular Bravais lattice with a 1 and a 2 denoting two primitive translations.Red, green and blue bonds between NN sites (i, j) carry three distinct Kitaev interactions, S x i S x j , S y i S y j and S z i S z j , respectively.With our choice of reference frame, the lattice sits on the (111) plane and a where a is a lattice constant.
The Hamiltonian has additional, hidden symmetries (selfdualities) at special points in parameter space, see discussion at the end of Sec.II C and Sec.II D.

C. Three-sublattice dualities
Similarly to the JK-model on other lattices, [2,22,24,26,[64][65][66][67] the Hamiltonian (1) supports duality transformations -also referred to as the Klein duality [22] -which preserve the form of the Hamiltonian but alter the value of the parameters J and K.As shown in Ref. [22], the main difference with the other lattices is that here the duality transformations involve three sublattices A, B and C (see Fig. 1) instead of four.There are three such transformations, which we denote by D A , D B and D C , and consist of a combination of π-rotations around the x-, yor z-axis, depending on the sublattice index.For example, D C arises by a product of C 2y rotations for the A sites and C 2x rotations for the B sites, Under D C the various spin operators S i transform to as shown in Fig. 3 (a).Similarly, D A and D B are given by Under D A , D B or D C , the Hamiltonian maps to i.e., the general form of ( 1) is preserved but K and J map to The duality transformations allow to identify special symmetry points of the parameter space and relate different regions of the phase diagram of Fig. 2 (which will be discussed in detail in Sec.III below).Specifically, the region 'IA' of Fig. 2 maps to the region 'IB', while the region 'IIA' maps to 'IIB'.Moreover, the Heisenberg points ψ = 0 and π, where K = 0, map to the dual Heisenberg points ψ = π − arctan 2 0.6475π and − arctan 2 1.6475π, respectively, where K = 0.These dual points have therefore a hidden SU(2) rotation symmetry.The Kitaev points ψ = ±π/2 are self-dual points, since at these points J = J = 0, K = K and H = H, i.e., the transformations become symmetries.

D. Self-dualities at the Kitaev points
As it turns out, the two Kitaev points have many more selfdualities (than D A , D B , D C ), in fact, a subextensive number of them.The presence of these symmetries has not been recognized in previous studies, and explains naturally the absence of classical and quantum order by disorder at these special points [31,32].
More specifically, there are in total 2 3L self-duality transformations that map (J = 0, K) → (J = 0, K), where 3L is the total number of lines in the lattice.The general form of the transformed spin operators under these operations is shown in Fig. 3 (b).Each transformation is characterized by a set of 3L Ising-like variables ±1, one for each line of the lattice.In Fig. 3  Similar subextensive symmetries appear in the quantum compass model on the square [68,69], cubic [60] and honeycomb [70] lattices, although the analogous operations of flipping individual line variables amount to products of C 2 rotations around a fixed axis (and not around alternating axes as here).One should be cautious to differentiate here from the case of subextensive operations appearing in certain classical models but are absent from their quantum counterparts.This includes, for example, the Kitaev model on the triangular lattice [26] and the K 1 -K 2 model on the honeycomb lattice [66].In these cases, the subextensive operations involve flipping one Cartesian spin component, which is a valid operation only for classical spins [71], whereas here the transformed spin components shown in Fig. 3 (b) always correspond to flipping two components.Now, the important point about the above subextensive quantum symmetries is that they affect only a line of spin sites, and are therefore intermediate between global operations and local (single-site) operations.A generalization of Elitzur's theorem [72] by Batista and Nussinov [73] then asserts that in 2D such symmetries can be broken spontaneously only at zero temperature.This implies that the kagome Kitaev model can host a zero-temperature long-range ordered phase, and the latter can be anyone among 2 3L different degenerate GSs.Diagnosing this order in finite-size calculations, however, is far from straightforward, as we explain in detail in Sec.V C.
We note finally that the 3-sublattice dualities D A , D B and D C are special members of the subextensive family of selfdualities.Specifically, they correspond to the following choice of χ , η and ζ ( = 1 • • • L): A simple modification of the above self-dualities gives rise to a sub-extensive family of transformations that map one Kitaev point to another, namely (J = 0, K) → (J = 0, −K).The corresponding form of the transformed spin operators in this family of dualities is shown in Fig. 3 (c).Compared to the selfdualities shown in Fig. 3 (b), here the signs of the 3L numbers χ , η and ζ ( = 1 • • • L) alternate between +1 and −1 along their corresponding lines.Doing this for all lines results in effectively reversing the sign of the Kitaev coupling for all NN bonds.We learn therefore that the physics of the AF Kitaev model can be mapped to the physics of the ferromagnetic (FM) Kitaev model.

III. Classical phase diagram
Some aspects of the classical phase diagram have been discussed by Kimchi and Vishwanath [22], and a more complete characterization has been given more recently by Morita, Kishimoto and Tohyama [31] starting by the classical minimum of a triangular unit cell and tiling the solution to the lattice.We shall give here a complementary picture based on the Luttinger-Tisza (LT) approach [61][62][63]74], which leads naturally to a parametrization of the GSs that will prove convenient for the discussion of the order by disorder effect.

A. General setting of LT approach
In the LT approach one replaces the problem of minimizing the energy E under N 'hard' spin length constraints (S 2 i = S 2 , i = 1-N, where N is the total number of sites) by the much simpler problem of minimizing E under a single 'soft' constraint i S 2 i = NS 2 .The latter problem then reduces to the diagonalization of a coupling matrix Λ, whose lowest eigenvalue corresponds to the lowest energy of the soft problem.It then follows [61][62][63]74] that if an eigenvector corresponding to the minimum eigenvalue can be used to construct a configuration that satisfies the N original constraints then that configuration must be one of the GSs of the hard problem.As we will now show, this approach works successfully in the present model.
We begin by labelling the spin sites by (R, µ), where R = n 1 a 1 +n 2 a 2 (where n 1 and n 2 are integers) denotes the positions of the unit cell, and µ ∈ {A, B, C} is the sublattice index.We next go to momentum space and define where k belongs to the first Brillouin zone of the triangular Bravais lattice.We then rewrite E/N in a matrix form where S k is the 9 × 1 vector where the 9 × 9 coupling matrix Λ k has a block-diagonal form , and the 3 × 3 matrices Λ (α) k (α = x, y, z) are provide in App. A. The coupling matrix is hermitian and has therefore a complete set of orthonormal eigenvectors V k,ν (with ν = 1-9) satisfying We can then decompose where the coefficients c k,ν are defined via Equation ( 13) then tells us that we can saturate the minimum of the energy by using the modes at the special points (k * , ν * ) associated with the minimum eigenvalue λ min , i.e., to replace The energy then becomes where in the second step we used the soft constraint in momentum space, k * ,ν * |c k * ,ν * | 2 = 3S 2 .As we show next, the coefficients c k * ,ν * can be chosen in such a way as to satisfy the spin length constraints.

B. Classical ground states inside the regions IA and IB
Inside the regions IA and IB the minimum eigenvalue is achieved along the three special lines that connect opposite M points of the Brillouin zone, see representative case shown in Fig. 4.These lines are denoted by l x , l y and l z in the upper inset (yellow) hexagon of Fig. 2. The corresponding eigenvectors take the following simple form We can then build GSs by linearly combining these eigenvectors, leading to where i.e., the coefficients χ R , η R and ζ R are constrained to +1 or −1.
Note further that from the extensive set of different choices of χ R as we vary R (and similarly for η R and ζ R ), only a subextensive subset are independent.This stems from the fact that for δ 'xx' lines : since the line l x of the Brillouin zone is vertical to the direction of the 'xx' lines (see Figs.
Flipping the sign of one of these coefficients amounts to flipping the associated component for all sites on the corresponding line.Clearly, the total number of states is 2 3L .This type of sub-extensive degeneracy is accidental, except for the AF Kitaev point where the degeneracy arises from the sub-extensive self-dualities discussed in Sec.II D and shown in Fig. 3

(b).
A few comments on the general structure of the GSs are in order here.First, the total spin in each unit cell of the lattice vanishes, i.e., as can be seen from Eq. ( 18).This is precisely the condition that minimizes the classical energy of the KHAF [40,75,76].Therefore, the subextensive set of classical GSs of the regions IA and IB (including the AF Kitaev point) is a subset of the extensive classical GSs of the KHAF.The three spins in each given unit cell are coplanar (with (χ R , −η R , ζ R ) being the spin plane), and form an angle of 120 • relative to each other.However the GSs are globally non-coplanar, in general.
Second, the GSs remain non-collinear even at the AF Kitaev point.This is qualitatively different from what happens in the square compass [60,69] and triangular Kitaev model [24,26], where spins align along one of the three Cartesian axes leading to non-zero spin correlations along the corresponding type of lines and vanishing correlations between different lines.Here such a decoupling between lines does not occur because of the special corner-sharing-triangle topology of the kagome.
Third, the spin components along a given Cartesian axis are correlated antiferromagnetically along the corresponding 'xx', 'yy' or 'zz' line, see for example the alternation of zcomponents along the line associated with ζ 1 in Fig. 5 (a).
Finally, the expression for the classical GS energy inside the regions IA and IB is which can be easily veriified by looking at the general structure of the GSs in Fig. 5 (a) and by noting that each bond contributes an energy of −(K + J)S 2 /2.

C. Classical ground states inside the regions IIA and IIB
Let us now turn to the GSs inside the regions IIA and IIB of Fig. 2. We will first discuss what happens away from the special points ψ = 0.6475π, 3π/2 and 2π.Here the minima of the eigenvalue spectrum of Λ k sits at the Γ point k = 0 of the Brillouin zone, with energy per site region IIA-IIB: The minimum eigenvalue is 3-fold degenerate, with corresponding eigenvectors: where ξ = 1/ √ 2+λ 2 and λ = cos ψ+sin ψ+ 5+4 cos 2ψ+sin 2ψ / (2 cos ψ) .(26) A linear combination of these modes with coefficients χ 1 , η 1 and ζ 1 (or c Γ,1 , c Γ,2 and c Γ,3 in the above notation) delivers the following structure for the GSs: These expressions agree with the results of Morita et al [31] (where λ is denoted by f ).
Note that, unlike the regions IA and IB, here the GSs evolve with ψ, because λ changes with ψ, see Fig. 6 (blue curve).Note also that the states are uniform (they do not depend on R), they comprise three magnetic sublattices (A, B and C) and the angle between NN spins is common for all bonds, The evolution of these correlations with ψ is shown in Fig. 6 (red curve).They are large and positive (i.e., FM) in most of the region IIA, they turn AF at tan ψ 2 = 2− √ 85 9 (ψ 1.5695π), and eventually approach the value −1/2 (i.e., the angle between NN spins tends to 120 • ) as ψ → 2π.
Let us next discuss the choice of the coefficients χ 1 , η 1 and ζ 1 .Imposing the spin length constraint on A, B and C spins gives the conditions which in turn lead to span the two-dimensional S 2 surface of the unit sphere, while elsewhere they can only take the values ±1, i.e., they become Ising variables.So the LT method delivers the degeneracy expected for the FM and dual FM Heisenberg points.Away from these points the LT method gives eight GSs (except at the AF and dual AF Heisenberg points, see below).These eight states are related to each other by the symmetry operations of the dihedral group D 2 and time reversal.
D. Classical ground states at special points (ψ = 3π 2 , 0, 0.6475π) The FM Kitaev point is special in that the minimum of the eigenvalue spectrum of Λ k is achieved on the lines l x , l y and l z of the BZ, and not just on the Γ point (see hexagon in Fig. 2).This feature leads to a subextensive number of GSs, eight of which are the ones resulting from Eq. ( 27) at λ = 0.
The whole set of GSs can arise from the subextensive GSs of the AF Kitaev point by simply applying the K → −K duality transformations of Sec.II E. According to Fig. 3 (c), we must simply flip the sign of the x-, yand z-components on every second site along each 'xx', 'yy' and 'zz' line, respectively.The resulting 2 3L states are shown in Fig. 5 (b).These states are again non-coplanar but the angle between any two NN spins is now 60 • and the spin components along a given Cartesian axis are correlated ferromagnetically along the corresponding 'xx', 'yy' or 'zz' line.As in the case of the AF Kitaev point, the 2 3L degeneracy is not accidental but stems from the Kitaev self-dualities of Sec.II D, which amount to flipping the sign of any of the numbers χ , η or ζ in Fig. 5 (b).The end points of the region IIA-IIB are also special in that the minimum of the eigenvalue spectrum of Λ k is achieved on the whole BZ and not just on the Γ point.For the KHAF this feature reflects the presence of an extensive manifold of GSs [75].These states satisfy the condition that the total spin in each triangle vanishes, i.e., S R,A + S R,B + S R,C = 0 , ∀R .
The resulting manifold includes an infinite subset of coplanar states as well as an infinite subset of non-coplanar states [40,75,76].The eight uniform states resulting from Eq. ( 18) for χ R , η R and ζ R independent of R, as well as the eight states of Eq. ( 27) for λ = 2 belong to this subset.For ψ = 2π, for example, the sublattices A, B and C approach the directions , respectively, which corresponds to a uniform 120 • coplanar configuration, see last inset graphic of Fig. 6.The situation for the dual AF Heisenberg point follows by the dualities D A , D B and D C .
IV. Semiclassical analysis: Quantum order-by-disorder in the regions IA-IB We now move onto the quantum case, which we first try to approach by a semiclassical 1/S approach.In particular, let us first consider the question of the lifting of the accidental sub-extensive degeneracy inside the regions IA and IB (except the AF Kitaev point where the degeneracy is symmetry related as mentioned above).To address this question we follow the so-called real space perturbation theory (RSPT) approach [77][78][79][80][81][82].In this approach one introduces a local axes frame along the classical spin directions of each spin site i, and then splits the Hamiltonian H into a diagonal part H 0 that includes fluctuations in the local field, and a perturbation V = H − H 0 , which couples off-diagonal fluctuations on different sites (see App. B).
In the present case, these fluctuations give rise to effective couplings between the Ising-like variables χ , η and ζ that parametrize the classical GSs of the regions IA and IB.As it turns out, the leading effective couplings appear in fourthorder perturbation theory, with lower order terms giving rise to global energy shifts.The fourth-order virtual processes involve five-site clusters, like the one shown in Fig. 7.In this cluster, a coupling between the Ising-like variables ζ 1 and ζ 2 arises from virtual processes generated by the parts of V living on the bonds a, b, c and d of Fig. 7.With the initial state being the same as the final state, the virtual process must involve either the bonds {a, c}, or {a, b}, or {d, b}, or {d, c}.The individual contributions from these four types of processes (evaluated for S = 1/2) are : {a, b} or {d, c} : δE (4)  1 = − 7 288 where 'cst' denote global constants.Disregarding these constants (and the ones arising from second-and third-order processes) and adding all contributions gives the total effective coupling involving ζ 1 and ζ 2 : where Similar effective couplings arise for all NN lines of any type, leading to the order-by-disorder effective Hamiltonian In addition, the sign of J eff is negative everywhere inside the regions IA and IB of the classical phase diagram.Therefore, fourth-order virtual processes tend to align the Ising-like variables ferromagnetically.So the leading quantum fluctuations lift the 2 3L subextensive accidental degeneracy of the regions IA and IB (except for the AF Kitaev point) and select the eight states with χ = χ 1 , η = η 1 and ζ = ζ 1 for all = 1 • • • L. These states are uniform (i.e., S R,A = S A , S R,B = S B and S R,C = S C for all R) and globally coplanar.The full quantum mechanical calculations presented below reveal that this leading orderby-disorder effect remains robust to all orders of perturbation theory, except near the AF Heisenberg point and its dual, where the physics is different.
The tendency of quantum fluctuations to select the above eight uniform states has been noticed previously in the numerical results of Morita et al [31,32].The RSPT analysis presented here unveils the nature of the leading virtual processes responsible for this order by disorder effect, and the effective Hamiltonian of Eq. ( 34) encapsulates this effect in a particularly simple form in terms of effective exchange couplings between emergent Ising-like variables.
An important comment for the AF Kitaev point is in order here.According to Eq. ( 33), the effective coupling J eff vanishes at ψ = π/2.This feature will in fact survive to all orders of perturbation theory, and is a direct consequence of the subextensive self-dualities of Sec.II D, which prohibit any coupling between the Ising variables.These symmetries therefore explain naturally the reported [31,32] absence of classical and quantum order by disorder at the Kitaev points.Still, this statement does not imply absence of long-range order in the thermodynamic limit, as a spontaneous breaking of the selfdualities is still possible at zero temperature [73], see also discussion in Sec.V C.

V. Quantum spins S = 1/2: Exact Diagonalization study
We now turn to the study of the quantum S = 1/2 version of the model (1).To this end, we have performed exact diagonalizations (ED) on finite-size clusters with periodic boundary conditions, and have examined the symmetry structure of the low-energy spectrum and the GS spin-spin correlation patterns as we vary ψ.The interpretation of the resulting data is unveiled by contrasting the results from two independent types of ED calculations: i) ED in the full basis, and ii) ED in the restricted variational basis of the orthonormalized classical GSs of the regions IA-IB and IIA-IIB.Taken together, the results show that the qualitative semiclassical 1/S picture presented above remains robust down to S = 1/2, except near the AF and dual AF Heisenberg points.

A. Preliminaries
The results presented here are obtained for a 24-site and a 27site cluster, with spanning vectors (T 1 , T 2 ) = (a 1 + 2a 2 , −3a 1 + 2a 2 ) and (T 1 , T 2 ) = (3a 1 , 3a 2 ), respectively.The 27-site cluster has the full point group C 3v symmetry of the infinite system, whereas the 24-site cluster has a lower (C 2v ) point group.The nine allowed momenta of the 27-site cluster are shown in the bottom right corner of Fig. 8, and include the Γ point, the two corners ±K of the Brillouin zone (labeled by K * ), and six momenta inside the Brillouin zone (labeled by q * ), which are related to each other by C 3v .Similarly, the eight allowed momenta of the 24-site cluster are shown in the side inset of Fig. 13 and include the Γ point, a single M point, and three pairs ±q 1 , ±q 2 and ±q 3 (not related to each other by any symmetry) inside the Brillouin zone.
In our exact diagonalizations we have implemented the symmetries under translations, real space inversion as well as spin inversion (i.e., global π-rotation around the x-axis in spin space).The presented spectra therefore carry quantum numbers associated with the momentum k, the parity (even or odd, denoted by 'e' and 'o', respectively) under real space inversion and the parity under spin inversion (even or odd, denoted by 'Sze' and 'Szo', respectively).For the 27-site cluster, the spin inversion leads to an extra twofold degeneracy due to Kramers theorem (i.e., the sectors 'Sze' and 'Szo' are degenerate).
The ED calculations in the restricted basis of classical GSs are performed as follows.We first generate the set of relevant classical states (depending on the region of ψ, see below Here α runs from one up to the number of classical GSs, N is the number of spins in the cluster, and (θ i,α , φ i,α ) are the spherical angles parametrizing the direction Ω i,α of the i-th spin in the α-th state.The restricted basis {|α } can then be orthonormalized using the overlap matrix O, whose matrix elements are given by O αβ = α|β .We have checked numerically that the rank of this matrix equals its dimensionality (i.e., the states {|α } are linearly independent).The variational problem then reduces to diagonalizing the effective Hamiltonian where the matrix elements of H inside the basis {|α } can be found by splitting H into individual bond terms H i j and For the 27-site cluster (L = 3) and for the region inside IA-IB, the restricted classical basis includes all 2 3L = 512 states of Fig. 5 (a).For the region inside IIA-IIB, the basis includes the eight uniform states of Eq. ( 27), as well as the 512 states of Fig. 5 (b) which become relevant close to the FM Kitaev point.

B. Low-energy spectra
We are now ready to examine the ED spectra.The data from the two clusters are qualitatively consistent with each other, so here we focus on the spectrum of the larger, 27-site cluster, shown in Fig. 8 (the spectrum of the 24-site cluster is shown in Fig. 13

of App. C).
First of all, due to the dualities D A , D B or D C , the spectrum in the regions IB and IIB can be mapped to the spectrum in the regions IA and IIA, respectively, by a rescaling [see Eq. ( 7)] Furthermore, since the duality transformations are uniform and do not break real space inversion, the momentum and parity quantum numbers of the spectra are retained in the mapping.Likewise, the spectra at the isolated points ψ = ±π/2 are identical to each other, as can be clearly seen in Fig. 8.This is due to the duality that maps K → −K, mentioned in Sec.II E.
Second, apart from the vicinities of the AF and dual AF Heisenberg points (see below), the spectra in the full and the restricted basis are qualitatively the same, both in the multiplicities of the levels and the symmetry quantum numbers, up to a relatively high excitation energy.This is the first strong spectral evidence that the semiclassical picture survives down to S = 1/2, except near the AF and dual AF Heisenberg points.
Let us examine some individual subregions, starting with the vicinity of the AF Kitaev point.Here we find a characteristic sequence of highly degenerate levels consisting of both zero-and finite-momenta states.Exactly at ψ = π/2, we find, above the fourfold degenerate GS (decomposing into 2Γ.e), a 36-fold first-excited level (decomposing into 6Γ.e ⊕ 2q * ), followed by a 108-fold level (decomposing into 6Γ.e ⊕ 6q * ⊕ 6K * ), followed by another 108-fold level (decomposing into 4Γ.e⊕2Γ.o⊕6q* ⊕6K * ), and so on.The qualitative agreement [84] with the variational ED data shown in the left panel of Fig. 8 (b) demonstrates that this characteristic pattern is a manifestation of expected tunnelling among the 512 classical GSs due to the finite value of S = 1/2 (the tunnelling goes to zero and the levels collapse to each other in the limit S → ∞).Moreover, the observed high multiplicities of the levels arise from the fact that the Kitaev point features a nontrivial symmetry group, consisting of a subextensive number of operations (512 for the 27-site cluster), which in turn leads to irreducible representations of very high dimensionality.
For completeness, we have carried out an independent symmetry analysis of the 512 classical GSs of the 27-site cluster Real-space pattern of GS spin-spin correlations S z i S z j in the 27-site cluster (enclosed by dashed lines) for four representative points of the region IA (the corresponding data for S x i S x j and S y i S y j can be obtained from S z i S z j by three-fold rotations).The two rows of panels show data obtained using ED in the full basis (a-d) and ED in the restricted classical basis (e-h) of the 512 states of Fig. 5 (a).The black open circle at each panel denotes the reference site i.Positive (negative) correlations are shown by filled blue (filled red) circles, whose radius scales with the magnitude of the correlation.The strongest correlation, corresponding to NN sites in (d), is about -0.137; the classical value is −S 2 /2 = −0.125,see Eq. ( 41).Incidentally, this number is very close to the exact result -0.13 for the NN correlations in the Kitaev Honeycomb model in the thermodynamic limit [83].From these 512 states, only 256 are visible in the energy range shown in the full basis ED results of Fig. 8 (a), exactly half from each sector.The remaining half reside at higher energies, and some of them mix with states outside the classical manifold.
The FM Kitaev point shows an identical pattern of highly degenerate levels, but these levels go up in energy very fast as we depart from ψ = 3π/2, in contrast to what happens around ψ = π/2.This qualitative difference originates from the fact that the respective 512 states are classical GSs only at ψ = 3π/2, and away from this point these get replaced by the 8 states of Eq. ( 27).Moreover, this difference is captured by the variational energy spectra shown in right panel of Fig. 8 (b).
Next, in the regions extending roughly between ψ ∼ 0.7π and 1.45π and between ψ ∼ 1.55π and 1.85π, the low-energy spectrum consists entirely of zero-momentum states which are well separated in energy from finite-momentum states.Furthermore, all low-lying states have even parity with respect to real space inversion.These features suggest that the system does not break the translation and inversion symmetry in these regions, in agreement with the predictions of the semiclassical analysis (given that the states of Eq. ( 27) as well as the 28 states of the S = 27/2 GS manifold of the FM Heisenberg point are all uniform and even).
Finally, let us discuss what happens in the vicinities of the AF and dual AF Heisenberg points.Figure 8 (a) shows a rapid rearrangement of the low-energy spectrum, with a very large number of states coming down in energy near these points, leading to a very dense spectrum in all possible momentum sectors.These states are clearly not related to the 512 states around the AF Kitaev point, as can be seen by a direct comparison to the left panel of Fig. 8 (b).The dense excitations are in fact not unexpected, as it is well known that the lowlying spectrum of the KHAF features an extensive number of low-lying excitations, in both the singlet and the higher spin sectors [42,43,45,59,85,86].It is further known that the origin of these low-lying states is actually not classical [49,86] but strongly quantum [41,43,46,58,59,87,88].
The competition between the quantum low-lying states of the KHAF and the semiclassical states favoured by the Kitaev Hamiltonian (in particular, the eight states stabilized via the order-by-disorder mechanism in IA and IB, see Sec.IV) appears to give rise to a phase transition between the two, at some critical point K/J (and a similar transition around the dual point ψ = 0.6475π).A rough estimate for this point can be deduced from the observation that the rapid rearrangement of the spectrum happens around ψ ∼ 0.25π in both 27-and 24-site clusters.The spin-spin correlation data from these two clusters alone give a similar rough estimate, as we show in the next section.This issue requires further discussion, however, and we will return to it in Sec.VI.

C. Spin-spin correlations
Further insights arise by examining the GS spin-spin correlations.We will focus on results obtained for the 27-site cluster, and we will present correlations of the type S z i S z j , since the corresponding patterns for S x i S x j and S y i S y j arise by threefold rotations.Real-space pattern of spin-spin correlations S z i S z j for the 27-site cluster (enclosed by dashed lines) for four representative points of the region IIB.The two rows of panels show data obtained from ED in the full basis (a-d) and ED in the restricted classical basis of the 512 states of the FM Kitaev point plus the 8 states of Eq. ( 27) (e-h).The black open circle at each panel denotes the reference site i.Positive (negative) correlations are shown by filled blue (filled red) circles, whose radius scales with the magnitude of the correlation.The strongest correlation, corresponding to NN sites in (a), is about 0.137; the classical value is S 2 /2 = 0.125, see Fig. 5 (b).

Regions IA-IB
We first examine the correlation patterns at four representative points inside the region IA, ψ = 0.1π, 0.3π, 0.4π and 0.5π, see Fig. 9 (a-d).The corresponding correlations inside IB can be deduced via the dualities D A , D B , and D C .
We begin by analyzing the results for the Kitaev point [Fig.9 (d)].Here we observe a strongly anisotropic profile with correlations being nonzero only along the single zz-line containing the reference site i (black open circle) and alternating in sign from one site to the next.This characteristic pattern is similar to what happens in the compass model in the square lattice [60,69], and its origin is the following.According to Fig. 5 (a), for any classical GS, the correlations S z i S z j are nonzero along all zz lines and not just along a single zz line.Take, for example, two spin sites in Fig. 5 (a) one (S i ) sitting on the horizontal line labeled by ζ 1 and the other (S j ) on the line labeled by ζ 2 .From the general form of the GSs, we find that the classical S z i S z j correlation is given by: where d i j is the distance between i and j.So, for fixed i and j, the overall sign of the correlation depends on the relative signs of ζ 1 and ζ 2 .As such, it is positive for half the states in the classical GS manifold, and negative for the other half.Now, the quantum GS of the 27-site cluster at ψ = π/2 is equal, to a very good approximation, to a symmetric superposition of all 512 classical states of Fig. 5 (a).This is corroborated by the agreement between Fig. 8 (a) and (b, left panel), as well as the agreement between Fig. 9 (d) and (h) for the correlations (Note, in particular, the striking agreement in the NN correlations S z i S z j = −0.136999from ED in the full basis and −0.133779 from ED in the variational basis).As a result, the above type of correlations (i.e., between spins belonging to different zz lines) average out in this superposition.This cancellation does not occur when the two sites sit on the same zz line, since in that case the same for all members of the classical GS manifold.
It is important to emphasize that the absence of correlations between different zz lines is only a finite-size (or finite-S ) effect, as the system can still develop long-range order in the limit N → ∞ (or S → ∞, respectively) at zero temperature [73].The actual correlation pattern will then depend on the particular member of the manifold that is selected spontaneously.This is similar to what happens in the square lattice compass model at zero temperature [69].
Moving away from the AF Kitaev point, the correlations on zz lines not including the reference site become nonzero even for finite-sizes, as can be seen in Figs. 9 (c, g), (b, f) and (e).This happens because, away from the Kitaev point, the 512 states are classically degenerate by accident and not by symmetry.As such, the degeneracy is lifted by quantum fluctuations, leading to the eight uniform states discussed in Sec.IV.The fact that the pattern seen in Figs. 9 (c, g), (b, f) and (e) is qualitatively consistent with that of the eight uniform states is therefore a numerical confirmation of the order by disorder effect of Sec.IV.Similar numerical confirmation has been reported by Morita et al [31,32].
We should further note that the spatially anisotropic profiles shown in Figs. 9 (c,g) and (b,f), i.e., the fact that the correlations along the zz line of the reference site are much stronger than those on the remaining zz lines, reflects the presence of the  42), evaluated in the GS of the 27-site cluster and for all allowed momenta of the cluster.There are two sets of data, one from ED in the full basis (symbols) and the other from ED in the orthonormalized classical basis of each region.The maximum eigenvalue in (a) at ψ = π is 0.75 which is consistent with the value 3S 2 expected for the FM GS [see definition of M k in Eq. ( 42)].
remaining members of the classical manifold at low energies, at an energy scale ∝ J eff given by Eq. (33).Given that the latter grows with L, one expects that the influence of these remaining members will diminish with L, and the correlations to become eventually uniform in strength throughout the bulk of the system for L → ∞.In other words, the spatial anisotropy in the profiles of Figs. 9 (c,g) and (b,f) is a finite-size effect.The pattern seen in Fig. 9 (a) for ψ = 0.1π is different from this picture [as also seen by the contrast with Fig. 9 (e)], and is more consistent with the AF short-range correlations expected for the KHAF.The crossover to the long-range pattern begins around ψ ∼ 0.25π as seen from the spectra reorganization discussed above, as well as in the evolution of the eigenvalues of the correlation matrix discussed in Sec.V C 3.

Regions IIA-IIB
Next, we examine the correlation patterns of S z i S z j at four representative points inside the region IIB, ψ = 1.5π, 1.53π, 1.85π and 1.9π, see Fig. 10 (a-d).The corresponding correlations inside IIA can be deduced via the dualities D A , D B , and D C .The situation at ψ = 1.5π [Fig.10 (a)] is analogous to that at ψ = π/2, the only difference being that the amplitudes are FM along the zz line of the reference site.The correlation profiles away from the FM Kitaev point [Figs.10 (b,f), (c,g) and (h)] are all consistent with the eight uniform states of Eq. ( 27).In particular, the sign of the NN correlations switches from positive (FM) to negative from 1.53π to 1.85π, in full agreement with the classical picture, which predicts a sign change at ψ = 1.5695π.Finally, the contrast between the patterns of Figs. 10 (d) and (h) highlights once again the RVB physics in the vicinity of the AF Heisenberg point (note, in particular, the similar profiles in Figs. 9 (a) and 10 (d)).

Spin-spin correlation matrix
To shed further light into the crossover between the shortrange physics in the vicinity of the AF Heisenberg point and the long-range ordering favoured by the Kitaev interactions we examine the evolution of the eigenvalues of the GS spin-spin correlation matrix M k , defined as where µ, ν run over the three sublattices A, B and C of the kagome lattice, and k belongs to the set of the allowed momenta of the given cluster.
Figure 11 shows the evolution of the eigenvalues of M k for the 27-site cluster, for all allowed k points and for two sets of GSs, one obtained in the full basis (symbols) and one in the restricted classical basis (lines).Quite generally, the maximum eigenvalue of M Γ is much larger than that of M q * and M K * .But most importantly, the comparison between the ED data in the full basis and the ones in the variational classical basis show almost perfect agreement everywhere, except in the vicinities of the AF and dual AF Heisenberg points.This is consistent with the general picture obtained from the lowenergy spectra that the semiclassical physics remains robust down to S = 1/2, as long as we are sufficiently away from the AF and dual AF Heisenberg points.Furthermore, the crossover between the RVB physics of the KHAF and the long-range ordering favoured by the Kitaev couplings seems to occur again around ψ ∼ 0.25π, mirroring the onset of the low-lying spectral rearrangement discussed above.

VI. Discussion
In conjunction with previously known results [22,31,32], the present study concludes a rather consistent picture for the zero-temperature physics of the spin-1/2 Heisenberg-Kitaev model on the kagome lattice.Our ED results show a clear quantum-classical crossover from the RVB physics of the Heisenberg antiferromagnet to the regime stabilized by Kitaev anisotropy.The striking agreement between the ED results in the full basis and in the restricted classical GS basis in the Kitaev regime shows that this regime has a strong semiclassical character.By this we do not only mean that the GSs in this regime are qualitatively captured by the semiclassical limit, but also that the quantum (e.g., spin-wave) corrections are largely quenched.This is based on the fact that the ED in the restricted basis captures only the quantum tunneling between the different members of the basis, and not the fluctuation corrections from states outside the basis.The strong quenching of these corrections is essentially a manifestation of the large spin gaps generated by the Kitaev anisotropy.This aspect is also demonstrated in Fig. 12, which shows that, inside the Kitaev regime, the classical ground state energies are very close to the corresponding quantum energies found from ED on finite sizes.
On a broader perspective, this crossover originates in the qualitatively different degree of frustration between the AF Heisenberg model on one hand and the Kitaev model on the other.Indeed, unlike what happens in the honeycomb lattice [21], the Kitaev anisotropy in the kagome lattice does not come with any local conservation laws (explicit or emergent), but only with one-dimensional symmetries, much like what happens in the compass model on the square lattice [68,69].
The presence of the one-dimensional symmetries has not been recognized before and explains naturally a number of features, including the reported absence of classical and quantum order by disorder at the AF Kitaev point in finite-size calculations [31,32].We re-emphasize here that this absence of order by disorder does not imply absence of long-range order.Such an order is allowed by the generalized Elitzur's theorem by Batista and Nussinov [73] at zero temperature, but its diagnosis in finite-size calculations is not straightforward as discussed in Sec.V C.
Importantly, the above one-dimensional operations cease to be symmetries as soon as we depart from the Kitaev points.As a result, the subextensive classical GS degeneracy inside the region IA-IB of the parameter space (away from ψ = 0 and π/2) is accidental, and is therefore lifted by fluctuations, a result that has been highlighted by Morita et al [31].The perturbative analysis of Sec.IV provides a rather intuitive picture for this order by disorder effect, in terms of emergent interactions between collective, Ising-like variables describing whole lines of spins.This effective description leads naturally to the uniform coplanar states found in [31] and also revealed in our spin-spin correlation results.
Finally, we return to the important question of whether the above quantum-classical crossover will show up in the thermodynamic limit with a phase transition at a nonzero value of K. Namely, whether the RVB phase of the KHAF will survive in an extended region around K = 0. Our ED results for the low-energy spectra and spin-spin correlations show a qualitative change around |K| ∼ J for both 24-and 27-size clusters.By itself, this observation is not conclusive and a more systematic system-size dependence is needed to address this question.We believe however that energetic considerations alone provide a clear insight that there will be a transition at a nonzero value of K, irrespective of the nature of the GS at K = 0. On one hand, the strong quenching of fluctuation corrections mentioned above implies a weak renormalization of the GS energy in the semiclassical regime.Moreover, some of the classical GSs of the region IA-IB remain members of the classical GS manifold of the KHAF.On the other hand, we know that the low-energy states in the vicinity of K = 0 are not related to the semiclassical limit [49,86] and have a strong quantum character [41,43,46,58,59,87,88].This is explicitly demonstrated in Fig. 12 (b) which shows that the classical states become energetically very costly around the KHAF and its dual.The classical orders of the regime IA or IIB must therefore overcome this large quantum energy cost before they become ground states deep inside the Kitaev regime.Importantly, this large energy cost is not related to the energy gap (if any) above the quantum ground state of the KHAF.Hence, the RVB phase of the KHAF will survive in a finite parameter range irrespective of whether this phase has a gap or not, and irrespective of the actual nature of the ground state.
Acknowledgments: This work was supported by the U.S. Department of Energy, Office of Science, Basic Energy Sciences under Award No. DE-SC0018056.We also acknowledge the support of the Minnesota Supercomputing Institute (MSI) at the University of Minnesota.
Appendix A: Coupling matrices Λ (α)  k (α = x, y, z) As mentioned in the main text, the 9 × 9 coupling matrix Λ k of the Luttinger-Tisza method is block diagonal and has the following general form  42), evaluated in the GS of the 24-site cluster and for all allowed momenta of the cluster.There are two sets of data, one from ED in the full basis (symbols) and the other from ED in the orthonormalized classical basis of each region.with and f ν,k = 1 + e ik•a ν , ν = 1-3, a 3 = a 2 − a 1 .By symmetry, the three matrices Λ (α) k (α = x, y, z) are related to each other by threefold rotations, namely As such, these matrices have the same overall spectrum of eigenvalues.The eigenvalues of Λ (z) k can be written in the following analytical form We remind here the general setting of the real space perturbation expansion [77][78][79][80][81][82].In this expansion one calculates the energy renormalization of a given classical GS from quantum fluctuations, in the following way.First we fix the classical GS, and parametrize the direction of each spin in that state by a local axis e z i .We then choose two perpendicular axes e x i and e y i and define the combinations e ± i = 1 2 (e x i ± ie y i ).Rewriting spin operators in the local frame, S i = S z i e z i + S + i e − i + S − i e + i . (B1) gives the following general form of the Hamiltonian A zz i j S z i S z j + A z+ i j S z i S − j + A +z i j S − i S z j + A ++ i j S − i S − j + A +− i j S − i S + j + h.c., (B2) where A is a second-rank tensor containing the Heisenberg and Kitaev interactions with A zz i j = e z i • A i j • e z j , A z+ i j = e z i • A i j • e + j , etc. (B3) Next, we define the deviation operator n i = S − S z i and rewrite where E cl = S 2 /2 i j A zz i j is the classical energy, B j = −S i e z i • A i j = B j e z j is the local exchange field on site j (with magnitude B j = −S i A zz i j ), and we have also used the relation A z+ i j S − j = −B j e z j • e + j = 0 .(B5) We then set where stand for the the double spin-flip processes (V 1 ), the single spin-flip hopping (V 2 ), and the correlated, single spin-flip processes (V 3 ).The latter are analogous to the cubic magnon terms in the standard Holstein-Primakof spin-wave expansion.Equation (B6) form the basis for the RSPT, which proceeds via a perturbation theory in powers of V.In the present problem, the leading contributions to the order by disorder effect arise in fourth order, as we discussed in Sec.IV.
Appendix C: ED results for the 24-site cluster Figure 13 shows the low-energy spectrum of the 24-site cluster as obtained from ED in the full basis.Apart from certain unimportant details (related to the different symmetry sectors of this cluster), the results show the same qualitative picture with that from the 27-site cluster discussed in the main text.Similar qualitative agreement arises for the eigenvalues of the spin-spin correlation matrixof the 24-site cluster, which are shown in Fig. 14.

FIG. 2 .
FIG. 2. Classical phase diagram of the JK-model (1) parametrized by the angle ψ of Eq. (2).The inset hexagons show the k-points of the first Brillouin zone associated with the minimum eigenvalue of the coupling matrix Λ k .
(b), these variables are denoted by {χ 1 , χ 2 , • • • , χ L } for 'xx' lines (red), {η 1 , η 2 , • • • , η L } for 'yy' lines (green), and {ζ 1 , ζ 2 , • • • , ζ L } for 'zz' lines (blue).Note that changing the sign of any of these variables corresponds to flipping two spin components for every site residing on the corresponding line.For the horizontal line of Fig. 3 (b) associated with the number ζ 1 , for example, changing ζ 1 → −ζ 1 amounts to a product of alternating C 2y and C 2x rotations along this line.One can check that the Hamiltonian remains invariant under this operation.
1 and 2).Similarly, η R+δ = η R for δ 'yy' lines, and ζ R+δ = ζ R for δ 'zz' lines.Hence, for any given GS configuration, the coefficients χ R , η R and ζ R are fixed along individual 'xx', 'yy' or 'zz' lines, respectively.The general structure of the resulting states is shown in Fig.5 (a).Each state is characterized by a set of 3L Isinglike variables χ , η and ζ

FIG. 4 .
FIG. 4. Evolution of the three eigenvalues of Λ (z)k at ψ = 0.3π (region IA) along the first Brillouin zone symmetry paths shown in the inset.

FIG. 5 .
FIG. 5. General form of classical GSs inside the regions IA and IB of Fig. 2 (a) and at the FM Kitaev point (b).The numbers χ , η and ζ (where = 1 • • • L) are constrained to +1 or −1.The overall normalization prefactors of S / √ 2 have been omitted.

FIG. 8 .
FIG. 8.The low-energy spectrum of the 27-site cluster measured from the GS energy E 0 (ψ).The symmetry sectors associated with the different symbols are shown in the right along with their dimensionalities (which here include an extra factor of 2 coming from time reversal).Note that the horizontal axis is nonlinear in order to highlight the details of the spectrum near specific regions.(a) Upper panel: ED spectra in the full basis using the Lanczos algorithm.From the spectrum shown only the lowest 5 levels in each symmetry sector have converged to the requested accuracy of 10 −12 in absolute energy.(b) Lower panels: ED spectra in the variational classical basis consisting of the 512 states of the region IA-IB (left panel), or the 512 states of the FM Kitaev point plus the 8 states of Eq. (27) (right panel).
FIG.9.Real-space pattern of GS spin-spin correlations S z i S z j in the 27-site cluster (enclosed by dashed lines) for four representative points of the region IA (the corresponding data for S x i S x j and S y i S y j can be obtained from S z i S z j by three-fold rotations).The two rows of panels show data obtained using ED in the full basis (a-d) and ED in the restricted classical basis (e-h) of the 512 states of Fig.5 (a).The black open circle at each panel denotes the reference site i.Positive (negative) correlations are shown by filled blue (filled red) circles, whose radius scales with the magnitude of the correlation.The strongest correlation, corresponding to NN sites in (d), is about -0.137; the classical value is −S 2 /2 = −0.125,see Eq.(41).Incidentally, this number is very close to the exact result -0.13 for the NN correlations in the Kitaev Honeycomb model in the thermodynamic limit[83].
and have found the following symmetry decomposition 512 states of IA-IB → 36Γ.e ⊕ 4Γ.o ⊕ 28q * ⊕ 24K * .(39) FIG.10.Real-space pattern of spin-spin correlations S z i S z j for the 27-site cluster (enclosed by dashed lines) for four representative points of the region IIB.The two rows of panels show data obtained from ED in the full basis (a-d) and ED in the restricted classical basis of the 512 states of the FM Kitaev point plus the 8 states of Eq. (27) (e-h).The black open circle at each panel denotes the reference site i.Positive (negative) correlations are shown by filled blue (filled red) circles, whose radius scales with the magnitude of the correlation.The strongest correlation, corresponding to NN sites in (a), is about 0.137; the classical value is S 2 /2 = 0.125, see Fig.5 (b).

FIG. 11 .
FIG. 11.Evolution of the eigenvalues of the correlation matrix M k of Eq. (42), evaluated in the GS of the 27-site cluster and for all allowed momenta of the cluster.There are two sets of data, one from ED in the full basis (symbols) and the other from ED in the orthonormalized classical basis of each region.The maximum eigenvalue in (a) at ψ = π is 0.75 which is consistent with the value 3S 2 expected for the FM GS [see definition of M k in Eq. (42)].

FIG. 12 .
FIG. 12. (a) Evolution of ground state energy per site, E 0 /N, with ψ.The different curves correspond to: the classical result (black line), results from ED in the full basis of the 24-and the 27-site clusters (blue and dark yellow), and results from ED in the variational basis of the 27-site cluster (dashed red).(b) The difference between the classical energy per site and the energy obtained from ED in the full basis of the 27-site cluster.

FIG. 13 .FIG. 14 .
FIG.13.The low-energy spectrum of the 24-site cluster measured from the GS energy E 0 (ψ).The symmetry sectors associated with the different symbols are shown on the right along with their dimensionalities.Note that the horizontal axis is nonlinear in order to highlight the details of the spectrum near specific regions.The data shown are obtained from ED in the full basis using the Lanczos algorithm.From the spectrum shown only the lowest 5 levels in each symmetry sector have converged to the requested accuracy of 10 −12 in absolute energy.