Néel to dimer transition in spin-S antiferromagnets: Comparing bond operator theory with quantum Monte Carlo simulations for bilayer Heisenberg models
aa r X i v : . [ c ond - m a t . o t h e r] D ec N´eel to dimer transition in spin- S antiferromagnets: Comparing bond operator theorywith quantum Monte Carlo simulations for bilayer Heisenberg models R. Ganesh, Sergei V. Isakov, and Arun Paramekanti
1, 3 Department of Physics, University of Toronto, Toronto, Ontario, Canada M5S 1A7 Theoretische Physik, ETH Zurich, 8093 Zurich, Switzerland Canadian Institute for Advanced Research, Toronto, Ontario, Canada M5G 1Z8 (Dated: June 28, 2018)We study the N´eel to dimer transition driven by interlayer exchange coupling in spin- S Heisenbergantiferromagnets on bilayer square and honeycomb lattices for S =1 /
2, 1, 3 /
2. Using exact stochasticseries expansion quantum Monte Carlo (QMC) calculations, we find that the critical value of theinterlayer coupling, J ⊥ c [ S ], increases with increasing S , with clear evidence that the transition is inthe O (3) universality class for all S . Using bond operator mean field theory restricted to singlet andtriplet states, we find J ⊥ c [ S ] ∝ S ( S + 1), in qualitative accord with QMC, but the resulting J ⊥ c [ S ]is significantly smaller than the QMC value. For S =1 /
2, incorporating triplet-triplet interactionswithin a variational approach yields a critical interlayer coupling which agrees well with QMC. Forhigher spin, we argue that it is crucial to account for the high energy quintet modes, and show thatincluding these within a perturbative scheme leads to reasonable agreement with QMC results for S =1,3 /
2. We discuss the broad implications of our results for systems such as the triangular lattice S =1 dimer compound Ba Mn O and the S = 3 / Mn O (NO ). I. INTRODUCTION
Spin dimer compounds provide the simplest realiza-tion of a magnetically disordered ground state — onewhere strongly coupled pairs of spins entangle to formsinglets. Such systems are also of great interest since theyundergo magnetic field induced spin ordering via a quan-tum phase transition which is analogous to Bose-Einsteincondensation.
There are many well-known spin dimercompounds and well studied model Hamiltonians exhibiting such physics for S =1 / S = 1 triangular lattice dimer compound Ba Mn O ,point to a need to better understand higher spin gener-alizations and instabilities of such dimer states driven byinter-dimer interactions.Here, we explore this issue using a simple model whichexhibits such a dimerized ground state - the bilayerHeisenberg antiferromagnet with a Hamiltonian given by H = J ⊥ X i S i, · S i, + J X h ij i X ℓ =1 , S i,ℓ · S j,ℓ . (1)Here, i labels sites in one layer, ℓ = 1 , h ij i represents nearest neighbor pairs of spins withineach layer. For J ⊥ ≫ J , the first term in H dominatesand the ground state is composed of isolated interlayersinglets with S i, + S i, = 0 for every i . If J ⊥ ≪ J , thesystem will order magnetically provided the second (in-tralayer) term in the Hamiltonian is not too frustrated bythe lattice geometry. Here, we restrict our attention tocases where each layer is itself a bipartite lattice so thatthe ground state for J ⊥ ≪ J has long-range N´eel order.This model Hamiltonian has been extensively studied forthe S = 1 / . Effects of dis-order, induced by site dilution, have also been explored. However, there has been relatively little work on under-
FIG. 1: The interlayer dimer state on square and honeycombbilayers with singlet correlations between layers. The in-planeexchange J and the interplane exchange J ⊥ act as shown. standing the higher spin generalizations of the Hamil-tonian in Eq. 1. The spin-S square lattice bilayer hasbeen studied using Schwinger boson mean field theory and series expansions. Variants of this spin-S modelhave been argued to support novel spin solid phases inthree dimensions. In this paper, we study this Hamil-tonian for S = 1 / , , / and approximate analyses based on abond operator method generalized to arbitrary spin. Our main results are as follows. (i) Using exactstochastic series expansion quantum Monte Carlo (QMC)calculations, we find a N´eel to dimer transition with in-creasing J ⊥ that is in the O (3) universality class for allthe models we have studied. (ii) The critical value ofthe interlayer coupling, J ⊥ c [ S ], for the N´eel to dimertransition is found to increase for higher spin. Usinga bond operator mean field theory restricted to singletand triplet states, we find J ⊥ c [ S ] ∝ S ( S + 1), in qual-itative accord with QMC results. However, there is aquantitative discrepancy between the mean field J ⊥ c [ S ]and its QMC value, which becomes more significant forhigher spin. (iii) For S = 1 /
2, we show that takinginto account triplet-triplet interactions within a vari-ational approach brings the J ⊥ c [ S ] value close to theQMC result. For higher spin, we show that the domi-nant corrections to the critical point arise from the highenergy quintet modes and direct triplet-triplet interac-tions are less important. Incorporating the quintet ex-citations within a perturbative treatment is shown toyield a critical interlayer coupling which is in good agree-ment with QMC results for S = 1,3 /
2. We discuss thebroad implications of our results for high spin antifer-romagnets such as the triangular lattice S = 1 dimercompound Ba Mn O and the S = 3 / Bi Mn O (NO ).This paper is organized as follows. Section II containsresults from the QMC simulations on the phase diagramof the honeycomb and square lattice bilayer models for S = 1 /
2, 1, 3 /
2. In Section III, we outline the bond oper-ator formalism generalized to the case of spin-S. SectionIV gives bond operator mean field theory results for thesquare and honeycomb lattice models. Section V dis-cusses the variational approach that we use to take intoaccount corrections beyond mean field theory. Section VIanalyses the S = 1 / S > / II. QUANTUM MONTE CARLO SIMULATIONS
The bilayer honeycomb and square lattices are bipar-tite lattices which can be split into two sublattices A andB with every lattice bond being a link between sites be-longing to different sublattices. This ensures that there isno sign problem, so that the model in Eq. 1 is amenableto quantum Monte Carlo simulations. We perform quan-tum Monte Carlo simulations for S = 1 / , , / L = 12 , , , , L = 12 , , , ,
36) latticesusing the Stochastic Series Expansion algorithm. For S = 1 and S = 3 /
2, simulations are performed with mod-ified worm weights, which lead to a slightly more efficientalgorithm, as in Ref. 26. At large enough ratio J ⊥ /J ,the system undergoes a quantum phase transition from aN´eel state to a dimerized paramagnetic state. To locatequantum critical points, we perform finite size scalinganalysis of the superfluid density and the staggered mag-netization density squared. We measure the superfluiddensity ρ s by measuring winding number fluctuations ρ s = T h W + W i J , | m s | L β / ν J ⊥ / J L =12 L =16 L =24 L =32 L =40 | m s | L β / ν [ K - K c ] L ν ρ s L J ⊥ / J L =12 L =16 L =24 L =32 L =40 ρ s L [ K - K c ] L ν FIG. 2: Scaling of the superfluid density (upper panel) andof the staggered magnetization density squared (lower panel)for the S = 1 antiferromagnet on the bilayer square lattice.The curves cross at a distinct point around J ⊥ /J = 7 . z = 1, ν = 0 . β = 0 . J ⊥ c /J = 7 .
15. Lines guide theeye. The error bars are smaller than the symbol size if notvisible. where W , are the winding numbers in two spatial di-rections and T is the temperature. The staggered mag-netization density squared is given by | m s | = 3 " N N X i ( − p S zi , where ( − p = 1 for lattice sites from sublattice A,( − p = − N is thenumber of lattice sites. In the vicinity of a continuousphase transition the superfluid density scales as ρ s = L − d − z F ([ K − K c ] L /ν , T L z ) , where L is the linear system size, d = 2 is the dimension-ality of the system, T is the temperature, [ K − K c ] ≡ ρ s L J ⊥ / J L =12 L =16 L =24 L =32 L =40 ρ s L [ K - K c ] L ν FIG. 3: Scaling of the superfluid density for the S = 3 / J ⊥ /J = 13 . z = 1, ν = 0 . J ⊥ c /J = 13 . [( J ⊥ − J ⊥ c ) /J ] is the distance from the critical point J ⊥ c /J , and ν is the correlation length critical exponent.The staggered magnetization density squared scales as | m s | = L − β/ν M ([ K − K c ] L /ν , T L z ) , where β is the critical exponent. If one plots ρ s L z asa function of J ⊥ /J at large enough and fixed value of1 / ( T L z ) then the curves for different system sizes shouldcross at the critical point J ⊥ c /J . If one plots ρ s L z asa function of [ K − K c ] L /ν , with appropriately chosenvalues of the critical exponents and K c , the curves fordifferent systems sizes should collapse onto the universalcurve given by the function F . Similarly, | m s | L β/ν asa function of J ⊥ /J should have a distinct crossing pointat the critical point and | m s | L β/ν as a function of [ K − K c ] L /ν should collapse onto the universal curve given bythe function M . We perform simulations at fixed aspectratio T = J / L . A. Square lattice
The quantum critical point for the S = 1 / J ⊥ c /J = 2 . J ⊥ c /J = 7 . S = 1, and at J ⊥ c /J = 13 . S = 3 /
2. The data scale very well with the criticalexponents ν = 0 . β = 0 . O (3) univer-sality class for any value of spin. The crossing pointsand data collapse for S = 1 and S = 3 / | m s | L β / ν J ⊥ / J L =12 L =18 L =24 L =30 L =36 | m s | L β / ν [ K - K c ] L ν ρ s L J ⊥ / J L =12 L =18 L =24 L =30 L =36 ρ s L [ K - K c ] L ν FIG. 4: Scaling of the superfluid density (upper panel) and ofthe staggered magnetization density squared (lower panel) forthe S = 1 antiferromagnet on the bilayer honeycomb lattice.The curves cross at a distinct point around J ⊥ /J = 4 . z = 1, ν = 0 . β = 0 . J ⊥ c /J = 4 . the magnetization density squared for S = 3 / B. Honeycomb lattice
We find that that for the honeycomb lattice the quan-tum critical points are located at J ⊥ c /J = 1 . S = 1 / J ⊥ c /J = 4 . S = 1, and J ⊥ c /J =9 . S = 3 /
2. The data scale very well with thecritical exponents ν = 0 . β = 0 . O (3)universality class for any value of spin. The crossingpoints and data collapse for S = 1 and S = 3 / S = 3 / ρ s L J ⊥ / J L =12 L =18 L =24 L =30 L =36 ρ s L [ K - K c ] L ν FIG. 5: Scaling of the superfluid density for the S = 3 / J ⊥ /J = 9 . z = 1, ν = 0 . J ⊥ c /J = 9 . the data points are too noisy. ( J ⊥ / J ) c [ S ] S ( S +1) QMC (square)bond operator MFT (square)QMC (honeycomb)bond operator MFT (honeycomb) FIG. 6: J ⊥ c /J [ S ] as a function of S ( S + 1) for the bilayersquare and honeycomb lattices. Lines are linear fits. Notethat the curves cross approximately at S ( S + 1) = 0. C. Critical point as a function of spin
In Fig. 6, we show the quantum critical points J ⊥ c /J as functions of S ( S + 1) for the bilayer square and honey-comb lattices. We find that J ⊥ c /J [ S ] is a linear functionof S ( S + 1). In the following sections, we will attempt tomake sense of these QMC results using an extension ofthe bond operator theory of the dimerized state and itsinstability to N´eel order. III. BOND OPERATOR REPRESENTATION
An elegant approach which allows us to understand thephysics of the dimer ground state and its magnetic or-dering instabilities is the bond operator formalism whichwas first proposed for S = 1 / Inthis scheme, the spin operators are represented in a newbasis consisting of singlet and triplet states on the inter-layer bonds ( i, i, J = 0, the ground state consists of localizedsinglets on these bonds, with a gap J ⊥ to the tripletexcitations. A nonzero J allows a pair of neighbor-ing bonds ( i, i,
2) and ( j, j,
2) to exchange theirsinglet/triplet character. Such a ‘triplet hopping’ pro-cess converts the localized triplet modes into dispersing‘triplons’, with three-fold degenerate bands due to theunderlying SU (2) symmetry of the Hamiltonian. In thispicture, the dimer to N´eel transition is an O (3) transi-tion driven by the condensation of triplon modes at acertain wavevector where the dispersion minimum hitszero. Generalizations of this approach to spin-1 magnetshave been proposed earlier. Here, we adopt a recentgeneralization of the bond operator method to arbitraryspin to study bilayer Heisenberg antiferromagnets.In a spin- S bilayer system, in the limit J ⊥ ≫ J , wehave isolated interlayer bonds. The bond can be in one ofthe following states: a singlet, a 3-fold degenerate triplet,a 5-fold quintet, etc. We introduce one boson for each ofthese states: | s i i ≡ s † i | i , | t i,m ∈{− , , } i ≡ t † i,m | i , | q i,m ∈{− , ··· , } i ≡ q † i,m | i , ...The index i here runs over all interlayer bonds, and m labels the S z -component of the total spin on the inter-layer bond. These boson operators form the basis for abond operator representation. To restrict to the physicalHilbert space of spins, every interlayer bond should haveexactly one boson, s † i s i + X m = − , , t † i,m t i,m + X n = − , ··· , q † i,n q i,n + · · · = 1 . (2)In terms of bond operators, the exchange interaction onan interlayer bond is given by J ⊥ S i, · S i, = ε s s † i s i + ε t X m = − , , t † i,m t i,m + ε q X m = − , ··· , q † i,m q i,m + · · · (3)where ε s = − J ⊥ S ( S + 1), ε t = J ⊥ { − S ( S + 1) } , and ε q = J ⊥ { − S ( S + 1) } .Bond operator theory re-expresses the spin operatorsand their interactions in terms of these bond bosons. Inthe limit J ⊥ ≫ J , the singlets, triplets, quintets, etc.form a hierarchy with the energy spacing between eachtier of order J ⊥ . In this paper, we restrict our analysis tothe low energy subspace of singlets, triplets and quintetson a bond, and neglect higher spin states as they aremuch higher in energy.We first turn to the usual bond operator mean fieldtheory retaining only singlet and triplet modes, ignoringtriplet interactions and higher excited states and impos-ing the constraint in Eq. 2 on average. We then con-sider, in turn, the effect of triplet-triplet interactions for S = 1 / S > /
2. Forconvenience of notation, we henceforth set J = 1, thusmeasuring J ⊥ in units of J . IV. SINGLET-TRIPLET MEAN FIELD THEORY
At mean field level, the interlayer dimer state is de-scribed by a uniform condensate of the singlet bosons,with h s i i = h s † i i = ¯ s . Retaining only triplet excitations,the spin operators at each site are given by S + i,ℓ = ( − ℓ r S ( S + 1)3 ¯ s { t i, − − t † i, } + 1 √ { t † i, t i, + t † i, t i, − } , (4) S zi,ℓ = ( − ℓ r S ( S + 1)3 ¯ s { t i, + t † i, } + 12 { t † i, t i, − t † i, − t i, − } . (5)Using these expressions, the Hamiltonian takes theform H mf = ε s N ⊥ ¯ s + ε t X i,m t † i,m t i,m − µ X i X m t † i,m t i,m + ¯ s − ! + 2 S ( S + 1)3 ¯ s X h i,j i h { t i, + t † i, }{ t j, + t † j, } + (cid:16) { t i, − − t † i, }{ t † j, − − t j, } + h.c. (cid:17)i , (6)where µ is a Lagrange multiplier which enforces the con-straint in Eq. 2 on average. N ⊥ is the number of in-terlayer bonds. We have dropped quartic terms in thetriplet operators (which corresponds to ignoring triplet-triplet interactions).In the rest of this paper, we use the following two ba-sis sets to represent triplet states: {| t − i i , | t i i , | t i i } or {| t x i i , | t y i i , | t z i i } . The former basis labels states by thez-projection of spin. The latter labels each state by thedirection in which its spin projection is zero. We cango from one basis to another using | t i i = | t z i i and | t ± i i = ( ∓| t x i i − i | t y i i ) / √
2. Below, we will use theindex m to represent an element of the first basis and u to represent an element of the second. FIG. 7: Top view of bilayers. (Top) Square lattice with prim-itive lattice vectors ˆ x and ˆ y shown. (Bottom) Honeycomb lat-tice. The shaded region is the unit cell composed of two sites.Sites marked with a red circle belong to the A sublattice. Un-marked sites belong to the B sublattice. The primitive latticevectors ˆ a and ˆ b are shown. A. Square lattice bilayer
A top view of the square lattice bilayer with the rele-vant primitive lattice vectors is shown in Fig. 7. At meanfield level, the Hamiltonian of Eq. 1 may be written as H (0) = − J ⊥ N ⊥ S ( S + 1)¯ s − µ ¯ s N ⊥ + µN ⊥ − N ⊥ A X k ,u ∈{ x,y,z } ′ ψ † k ,u (cid:18) A + 2 ǫ k ǫ k ǫ k A + 2 ǫ k (cid:19) ψ k ,u , (7)where ψ k ,u = [ t k ,u t †− k ,u ] T . The primed summationindicates that if k is included in the sum, then − k isexcluded. The coefficients in the Hamiltonian matrix are A = J ⊥ { − S ( S + 1) } − µ, (8) ǫ k = 2 S ( S + 1)3 ¯ s (cos( k x ) + cos( k y )) . (9)Diagonalizing this Hamiltonian matrix by a bosonic Bo-goliubov transformation (see Appendix A), we obtaineigenvalues λ k = p A ( A + 4 ǫ k ) for the energies of the in-dependent ‘triplon’ modes. Each of these modes adds azero point contribution to the ground state energy, yield-ing E (0) = − J ⊥ N ⊥ S ( S + 1)¯ s − µ ¯ s N ⊥ + µN ⊥ − N ⊥ A X k ′ λ k . (10)We minimize this ground state energy with respect to µ and ¯ s , via ∂E (0) /∂µ = 0 and ∂E (0) /∂ ¯ s = 0. This yieldsthe two equations¯ s = 52 − N ⊥ X k ′ A + 2 ǫ k λ k , (11) µ = − J ⊥ S ( S + 1) + 6 N ⊥ X k ′ Aǫ k ¯ s λ k . (12)Using the values of ¯ s and µ thus obtained, we may calcu-late the gap to triplet excitations. The dimer-N´eel tran-sition occurs when the triplon gap vanishes at J ⊥ = J ⊥ c .We have explicitly checked that triplon condensationat k = ( π, π ) yields N´eel order on the bilayer. UsingEqns. 11,12 above, we arrive at the following two resultsat the critical point. (i) The value ¯ s at the dimer-N´eelcritical point is independent of spin and is given by¯ s = 52 − N ⊥ X k ′ k x + cos k y ) p k x + cos k y ) . (13)A numerical evaluation shows ¯ s c ≈ . J ⊥ c = S ( S +1) " − N ⊥ X k ′ p k x + cos k y ) . (14)A numerical evaluation yields J ⊥ c ≈ . S ( S + 1). For S = 1 /
2, this mean field result, J ⊥ c [ S = 1 / ≈ . and is slightly smaller thanthe QMC value. For higher spin, the mean field esti-mates, J ⊥ c [ S = 1] ≈ .
095 and J ⊥ c [ S = 3 / ≈ . J ⊥ c ∼ S ( S +1) has been suggested in Ref. 21 on the basisof series expansion calculations. Remarkably, as shownin Fig.6, this scaling relation derived from mean-field the-ory seems to be reasonably accurate even for exact QMCresults. B. Honeycomb lattice bilayer
The honeycomb lattice is composed of two interpene-trating triangular lattices, as shown in Fig.7. Operatorstherefore come with an additional sublattice index whichdistinguishes A and B sublattices. The mean field Hamil-tonian is given by H (0) = − N ⊥ J ⊥ S ( S + 1)¯ s − N ⊥ µ ¯ s + N ⊥ µ − N ⊥ C X k ,u ′ ψ † k ,u M k ψ k ,u , (15)where C = ( J ⊥ { − S ( S + 1) } − µ ). N ⊥ denotes thenumber of interlayer bonds in the honeycomb bilayer.The operator ψ k ,u and the Hamiltonian matrix M k are given by ψ k ,u = t k ,A,u t k ,B,u t †− k ,A,u t †− k ,B,u , M k = C β k β k β ∗ k C β ∗ k β k C β k β ∗ k β ∗ k C , (16)where β k = 2 S ( S +1)3 ¯ s γ k , with γ k = 1+ e − ik b + e − ik a − ik b ,and we have defined k a ≡ k · ˆ a and k b ≡ k · ˆ b . Diago-nalizing this Hamiltonian (see Appendix B), we obtaintwo eigenvalues for every k . The eigenvalues are givenby λ k , / = p C ∓ C | β k | . The mean field ground stateenergy is given by E (0) = − N ⊥ J ⊥ S ( S + 1)¯ s − N ⊥ µ ¯ s + N ⊥ µ − N ⊥ C X k ′ ( λ k , + λ k , ) . (17)As before, we demand ∂E (0) /∂µ = 0 and ∂E (0) /∂ ¯ s = 0.This leads to the two mean field equations¯ s = 52 − N ⊥ X k ′ (cid:20) C − | β k | λ k , + C + | β k | λ k , (cid:21) , (18) µ = − CS ( S + 1) N ⊥ X k ′ | γ k | (cid:20) λ k , − λ k , (cid:21) − J ⊥ S ( S + 1) . (19)Using the values of ¯ s and µ thus obtained, we calculatethe gap to triplet excitations. The dimer-N´eel transitionoccurs when the triplon gap vanishes at J ⊥ = J ⊥ c . Usingthe above equations, we arrive at the following two resultsat the critical point. (i) The value ¯ s at the dimer-N´eelcritical point is independent of spin and given by¯ s c = 52 + 32 N ⊥ X k ′ " | γ k | − p − | γ k | − | γ k | + 6 p | γ k | . (20)A numerical evaluation shows ¯ s c ≈ . J ⊥ c S ( S + 1) = 10 − N ⊥ X k ′ " p − | γ k | + 1 p | γ k | . (21)A numerical evaluation yields J ⊥ c ≈ . S ( S + 1). For S = 1 /
2, the mean field result, J ⊥ c [ S = 1 / ≈ . J ⊥ c [ S = 1] ≈ .
496 and J ⊥ c [ S = 3 / ≈ . J ⊥ c ∼ S ( S + 1) from mean field theory appears to bevalid even for the exact QMC results on the honeycomblattice. We have also explicitly checked that triplon con-densation of the mode with energy λ k , at momentum k = (0 ,
0) yields N´eel order on the honeycomb bilayer.
V. BEYOND MEAN FIELD THEORY:VARIATIONAL ANALYSIS
Corrections to the mean field Hamiltonian arise fromtriplet-triplet interactions, and coupling to higher spinobjects such as quintets and heptets. As a function of S,we find two regimes where two different correction termsdominate. For S = 1 /
2, the only correction stems fromtriplet-triplet interactions since higher spin states are ab-sent. For
S > /
2, the dominant correction arises fromcoupling to higher spin (quintet) states. Ordinarily, suchquintet terms can be ignored as the energy cost of ex-citing quintets is large; however, these terms scale as S as opposed to the S scaling of the triplet-triplet inter-actions and they play an increasingly important role forlarger S . These two correction terms are separately dis-cussed in the following two sections.Specifically, for the two regimes S = 1 / S > / H (0) → H (0) +∆ H and H (0) → H (0) +∆ H . We treat ∆ H as aperturbation acting upon the states of H (0) . As a varia-tional ansatz, we assume that the effect of the correctionterms is entirely accounted for by a renormalization of theparameters ¯ s and µ which enter the mean field Hamilto-nian, H (0) or H (0) . We choose µ to enforce single bosonoccupancy per site on average. The perturbations ∆ H ,for both regimes, preserve total boson number. Thus, itsuffices to evaluate total boson number using H (0) . Thisgives us the constraint¯ s + X i,m h t † i,m t i,m i = N ⊥ , (22)where the expectation value is evaluated with respect to H (0) . (For the honeycomb lattice case, there is an addi-tional sum over the sublattice degree of freedom in theabove equation). This leads precisely to the mean fieldnumber constraint in Eq. 11 or Eq. 18, which can nowbe used to determine µ . The parameter ¯ s is chosen tominimize the ground state energy, evaluated to leadingorder in perturbation theory. For S = 1 /
2, we find thatthe leading correction is obtained within first order per-turbation theory in ∆ H . For S > /
2, the dominantperturbing terms require us to go to second order in per-turbation theory. In the next two sections, we discussthese correction terms in detail.
VI. TRIPLET INTERACTION CORRECTIONSA. Triplet Interactions on square lattice
Staying within the singlet-triplet sector, the term wehave ignored in the mean field treatment is the triplet-triplet interaction term. For S = 1 /
2, there are no higherspin bosons beyond the singlet-triplet sector, so this is the only correction. For
S > /
2, this constitutes one termin a slew of correction terms. For any spin S , the tripletinteraction terms are given by∆ H (t) = − X h ij i X u, v, w, v ′ , w ′ ∈ { x, y, z } ǫ uvw t † i,v t i,w ǫ uv ′ w ′ t † j,v ′ t j,w ′ . (23)We note that there are no cubic terms in triplet opera-tors. As described in Ref. 33, this makes our bilayer prob-lem qualitatively different from other dimerized statessuch as the spin-1/2 staggered dimer on the square lat-tice. Typically, triplet-triplet interactions such as thoseof Eq. 23 are taken into account within a self-consistentHartree-Fock approximation. Here, we take the in-teractions in Eq. 23 to be a perturbation acting on H (0) and evaluate the first order correction to ground stateenergy. To this end, we decouple ∆ H (t) using bilinearsthat possess finite expectation values at the level of meanfield theory: h t † i,v t i + δ,w i ≡ δ v,w ρ, h t † i,v t † i + δ,w i ≡ δ v,w ∆ . (24)Here, i and i + δ are nearest neighbours on the squarelattice. Explicit expressions for ρ and ∆ are given inAppendix A. We note that ρ and ∆ are functions of thevariational parameters ¯ s and µ . The first order energycorrection due to triplet interactions is given by∆ E (t) = h ∆ H (t) i = 6 N ⊥ (cid:2) ρ − ∆ (cid:3) . (25)Thus, the variational energy of the ground state uponincluding the triplet interaction term is given by E (t) , var (¯ s, µ ) = E (0) + ∆ E (t) , (26)where E (0) is as defined in Eq. 10. The parameter ¯ s ischosen to minimize this energy. We find that the tripletinteractions reduce the stability of the dimer phase andshift J ⊥ c to larger values. For S = 1 /
2, this leads toa renormalized transition point J ⊥ c ≈ .
58, very closeto the QMC result. For
S > /
2, the renormalizationis too weak to account for the discrepancy between theearlier mean field result and the QMC data. These tripletcorrected results for the square lattice are summarized inTable I.
B. Honeycomb
The interaction between triplets on the honeycomb lat-tice is given by∆ H (t) = − X i,δ X u, v, w, v ′ , w ′ ∈ { x, y, z } ǫ uvw t † i,A,v t i,A,w × ǫ uv ′ w ′ t † i + δ,B,v ′ t i + δ,B,w ′ . (27)The operators δ are such that the sites ( i, A ) and ( i + δ, B )are nearest neighbours. This interaction term contributesto the ground state energy at first order in perturbationtheory. To evaluate this correction, we quadratically de-compose the interaction using the following two bilinears: h t † i,A,v t i + δ,B,w i ≡ δ v,w ρ, h t † i,A,v t † i + δ,B,w i ≡ δ v,w ∆ , (28)with the expectation values to be evaluated using the un-perturbed Hamiltonian H (0) . Explicit expressions for ρ and ∆ are given in Appendix B. The first order correctionto ground state energy is given by∆ E (t) = 92 N ⊥ [ ρ − ∆ ] . (29)The parameter ¯ s is chosen to minimize the energy E (t) , var (¯ s, µ ) = E (0) + ∆ E (t) . (30)As on the square lattice, we find that the triplet interac-tions reduce the stability of the dimer phase and shift J ⊥ c to larger values. For S = 1 /
2, this leads to a renormal-ized transition point J ⊥ c ≈ .
59, which is in reasonableagreement with the QMC result J ⊥ c = 1 . S > /
2, however, the renormalization is again too weakto account for the QMC data. These triplet corrected re-sults for the critical point on the honeycomb lattice aresummarized in Table II.
VII. QUINTET CORRECTIONS
In the previous section, we have seen that triplet cor-rection terms lead to a reasonably good agreement withQMC results for the dimer-N´eel quantum critical pointfor S = 1 /
2. However, they fail to account for the signifi-cant discrepancy between QMC and bond operator meanfield theory for
S > /
2. This leads to us to suspect thathigher order spin excitations on the dimer bonds mustbe responsible for this difference. Upon including quin-tet terms, the spin operators at a site i contain a largenumber of terms as given in Eq. (20) and Eq. (21) ofRef. 24 and reproduced in Appendix C for convenience.Using these spin expressions to rewrite the Hamilto-nian in Eq. 1, we find that correction terms beyond meanfield theory, including those involving quintet states, maybe grouped as∆ H = ˆ D tttt + ¯ s ˆ R ttq ( S )+ ˆ F ttqq ( S )+ ˆ G qqqq ( S ) . (31)The subscripts indicate the composition of the terms interms of bond operators. The scaling of each term with S is indicated in parentheses. For example, ˆ R ttq ( S ) iscomposed of terms which involve two triplet operatorsand one quintet operator, and the coefficients of theseterms scale as S . The term which we have accountedfor in the previous section is ˆ D tttt , which scales as S and contains four triplet operators. Na¨ıvely, terms involvingquintets should be less important due to the energy costof exciting quintets. However, we see from the aboveclassification of terms that the coefficients of ˆ R ttq ( S )and ˆ F ttqq ( S ) increase rapidly with increasing spin. Wefind that ˆ R ttq ( S ) is, in fact, the dominant contributionfor all S > /
2. (For the case of S = 1, we have explic-itly checked that this term dominates over triplet-tripletinteractions encoded by ˆ D tttt - see Table I). The termˆ F ttqq ( S ) is suppressed because it involves two quintetoperators which act on different sites. In our variationalscheme, this term will contribute to the ground state en-ergy at second order in perturbation theory. However,as the quintets are taken to be localized excitations, thisterm will involve intermediate states with two quintetexcitations. Therefore, it will contribute much less thanˆ R ttq ( S ).In the vicinity of the dimer-N´eel transition, we assertthat ˆ R ttq ( S ) will remain the dominant correction termfor any S > / ± not contributeat second order in perturbation theory, but will only ap-pear at higher order. As an illustration, upon includingheptets, the Hamiltonian can have a term of the form h † i,m q i,n t † j,m ′ t j,n ′ . Clearly, this term does not contributeto ground state energy at first or second order. In ad-dition, at whichever order it contributes, the energy de-nominators will involve large heptet excitation energieswhich will further suppress the energy contribution.In summary, in the vicinity of the dimer-N´eel transitionfor any value of S > /
2, the leading correction to bondoperator mean field theory comes from ¯ s ˆ R ttq ( S ). Wewrite ∆ H ( S> / ≈ ∆ H (q) ≡ ¯ s ˆ R ttq ( S ) . (32)Note that if we were to use a path integral approach to in-tegrate out the quintet excitations at this stage, we wouldbe led to an effective triplet-triplet interaction which isenhanced by a factor of S compared to the bare triplet-triplet term discussed in the previous section (althoughit would be divided by the quintet gap which scales as S near the N´eel to dimer transition). Here we follow adifferent route, similar in spirit, and treat this term per-turbatively assuming the quintet states to be local exci-tations. The energy cost of creating a quintet is given byEq. 3. We measure this energy cost from the Lagrangemultiplier µ , to get ε q − µ = J ⊥ { − S ( S + 1) } − µ (33)as the energy cost of a quintet excitation. A. Quintet corrections on the square lattice
The terms in ¯ s ˆ R ttq ( S ) may be organized as¯ s ˆ R ttq ( S ) = ¯ s X i X n = − , ··· , " q † i,n X δ ˆ T [ n ] i,i + δ + h.c. . (34)The operator ˆ T [ n ] i,i + δ is composed of triplet bilinears. Theindex δ sums over the four nearest neighbour vectors onthe square lattice. The explicit form of these operatorsis given in Appendix D. The operator ˆ R ttq ( S ) does notcontribute to ground state energy at first order, as it islinear in quintet operators. The energy correction, atsecond order, is given by∆ E (q) = ¯ s X σ =0 X i, n, δi ′ , n ′ , δ ′ h | q i ′ ,n ′ ( ˆ T [ n ] i ′ ,i ′ + δ ′ ) † | σ ih σ | q † i,n ˆ T [ n ] i,i + δ | i E − E σ . (35)The index σ sums over all excited states of H (0) , thevariational Hamiltonian. The only intermediate statesthat contribute are those with a single quintet. In ourvariational formalism, we take the quintets to be localexcitations. This constrains us to ( i = i ′ ), ( n = n ′ ). Thisleaves us with∆ E (q) =¯ s X ν =0 X i,n h | P δ ′ ( ˆ T [ n ] i,i + δ ′ ) † | ν ih ν | P δ ˆ T [ n ] i,i + δ | i E − E ν . (36)The intermediate states | ν i which contribute involve asingle quintet excitation. Within the triplet sector, atzero temperature, the intermediate states can have either(a) no triplon quasiparticles, or (b) two triplon quasi-particles. The contribution from states with no triplonquasiparticles vanishes due to global spin-rotational sym-metry of the Hamiltonian. The energy correction fromtwo triplon intermediate states is evaluated to obtain theenergy correction, ∆ E (q) . The complete expression isgiven in Appendix D.Being second order in ˆ R ttq ( S ), the energy correctionfrom quintet coupling na¨ıvely scales as S for large S .However, the energy denominator involves the energy ofquintet states which is proportional to J ⊥ . Close to the S QMC MFT MFT + MFT +triplet interactions quintet coupling0.5 2.5220(1) J ⊥ c on the square lattice from differentmethods for different values of S. MFT stands for mean fieldTheory. The QMC data for S=1/2 is from Ref. 11. Thecolumn ‘MFT+Triplet interactions’ gives variational resultsappropriate for S=1/2. The column ‘MFT+quintet coupling’gives variational results appropriate for S > / dimer-N´eel transition, at mean field level, J ⊥ approxi-mately scales as S for large S (see Eq. 14). We expectperturbative corrections to preserve this scaling of J ⊥ c with S . Thus, near the dimer-N´eel transition, ∆ E (q) scales as S /S ∼ S . The ground state energy to lead-ing order in perturbation theory is thus given by E ( S> / , var (¯ s, µ ) = E (0) + ∆ E (q) . (37)This energy is a function of ¯ s and µ . As discussed earlier, µ is tuned to enforce single boson occupancy per site,while ¯ s is chosen to minimize E ( S> / , var (¯ s, µ ).Having determined ¯ s and µ variationally, we can findthe gap to triplon excitations as a function of J ⊥ . TheDimer-N´eel transition is indicated by the vanishing ofthe triplon gap in the variationally obtained state. Assummarized in Table I, the renormalized critical pointsobtained in this manner are within 2 .
5% of the QMCresults. While the precise quantitative agreement is per-haps fortuitous, and will certainly change depending onthe nature of the approximations made, the importantproblem we have resolved is to show that the large dis-crepancy between QMC and simple bond operator meanfield theory for
S > / B. Quintet corrections on the honeycomb lattice
On the honeycomb lattice, the terms in ¯ s ˆ R ttq ( S ) maybe written as¯ s ˆ R ttq = ¯ s X i X n = − , ··· , " q † i,A,n X δ ˆ A [ n ] i,i + δ + h.c + q † i,B,n X δ ˆ B [ n ] i,i − δ + h.c. . (38)The operators ˆ A i,i + δ and ˆ B i,i − δ are triplet bilinears cen-tred on nearest neighbour bonds. We give their explicitforms in momentum space in Appendix E. The terms inˆ R ttq ( S ) contribute to ground state energy only at sec-ond order in perturbation theory. The energy correction0 S QMC MFT MFT + MFT +triplet interactions quintet coupling0.5 1.645(1) 1.312 1.588 -1 4.785(1) 3.498 3.774 4.80(9)1.5 9.194(3) 6.559 6.837 9.58(18)TABLE II: Value of J ⊥ c on the honeycomb lattice from dif-ferent methods for different values of S. MFT stands for meanfield Theory. may be written as∆ E (q) = ¯ s X σ =0 X i,n h | " q i,A,n X δ ′ ˆ A [ n ] i,i + δ ′ | σ i ×h σ | " q † i,A,n X δ ˆ A [ n ] i,i + δ | i / { E − E σ } + ( A → B ) , (39)where the index σ sums over all excited states of H (0) .As the terms in ˆ R ttq ( S ) involve one quintet operator,only intermediate states with a single occupied quintetstate will contribute.∆ E (q) = ¯ s X ν =0 X i,n h | P δ ′ ( ˆ A [ n ] i,i + δ ′ ) † | ν ih ν | P δ ˆ A [ n ] i,i + δ | i{ E − E ν } +¯ s X ν =0 X i,n h | P δ ′ ( ˆ B [ n ] i,i − δ ′ ) † | ν ih ν | P δ ˆ B [ n ] i,i − δ | i{ E − E ν } . (40)We evaluate these overlaps in momentum space, as de-scribed in Appendix E. The intermediate state | ν i couldhave either (i) no triplon quasiparticles, or (ii) two triplonquasiparticles. However, the contribution from stateswith no triplons vanishes due to global spin rotationalsymmetry. The explicit expression for ∆ E (q) is given inAppendix E.Thus, the energy of the ground state to leading orderin quintet coupling, is given by E ( S> / , var = E (0) + ∆ E (q) . (41)We choose ¯ s to minimize this energy. The vanishing ofthe triplet gap in the variationally determined state sig-nals the dimer-N´eel transition. Our results for J ⊥ c onthe honeycomb lattice are shown in Table II. The renor-malized critical points for S = 1 , / VIII. DISCUSSION
In this paper, we have studied the N´eel to dimer tran-sition in Heisenberg antiferromagnets on bilayer squareand honeycomb lattices for different spin values usingQMC and bond operator approaches. The critical bilayer exchange J ⊥ c scales as S ( S + 1) within, both, bond oper-ator mean field theory and QMC simulations. However,there is a systematic deviation between bond operatormean field theory and QMC, with the deviation itselfscaling as ∼ S . Our variational extension of bond oper-ator theory to include the dominant triplet and quintetexcitations successfully captures this systematic devia-tion and gives a more precise estimate of J ⊥ c .Bi Mn O (NO ) provides an example of a bilayerhoneycomb antiferromagnet with S = 3 /
2, wherestrong interlayer exchange couplings ∼ J have been in-ferred from electronic structure calculations . Despitethis strong bilayer coupling, our study indicates thatthis material would be deep in the N´eel ordered phaseif there are no other frustrating interactions. We arethus forced to attribute the observed lack of magneticorder in Bi Mn O (NO ) to frustration effects arisingfrom further neighbor couplings; such further neighborinteractions have been shown to drive a variety of newphases on the honeycomb lattice. One recent exam-ple of a dimer system with S = 1 is the triangular dimermaterial Ba Mn O . Our approach could be appliedto understand the triplon spectrum and the effect of quin-tet corrections in this material. In particular, our workshows that extracting exchange couplings from fitting ex-perimental data to bond operator mean field theory maynot yield precise estimates. In summary, our work pro-vides a starting point to think about the physics of highspin Heisenberg antiferromagnets in a variety of modelsystems and materials. Acknowledgments
This research was supported by the Canadian NSERC(RG, AP), an Ontario Early Researcher Award (RG,AP), and the Swiss HP C initiative (SVI). The quantumMonte Carlo simulations were performed on the Brutuscluster at ETH Zurich. RG thanks J. Rau, V. VijayShankar and C. M. Puetter for discussions.
Appendix A: Square bilayer: bosonic Bogoliubovtransformation
The MFT Hamiltonian of Eq. 7 is diagonalized by apseudounitary matrix, U k = cosh θ k sinh θ k sinh θ k cosh θ k ! . (A1)Imposing tanh 2 θ k = − ǫ k / ( A + 2 ǫ k ), we get ψ † k ,u A + 2 ǫ k ǫ k ǫ k A + 2 ǫ k ! ψ k ,u = φ † k ,u λ k λ k ! φ k ,u . (A2)1We have defined new quasiparticle operators given by ψ k ,u = U k φ k ,u so that t k ,u t †− k ,u ! = cosh θ k sinh θ k sinh θ k cosh θ k ! τ k ,u τ †− k ,u ! . (A3)The τ operators are the triplon quasiparticles. The bi-linears defined in Eq. 24, may be evaluated using theelements of U as follows: ρ = 14 N ⊥ X k ,δ h h t † k ,v t k ,v i e i k .δ i = 14 N ⊥ X k ′ (2 cos k x + 2 cos k y ) A + 2 ǫ k λ k , (A4)∆ = 14 N ⊥ X k ,δ h h t † k ,v t †− k ,v i e i k .δ i = 14 N ⊥ X k ′ (2 cos k x + 2 cos k y ) ( − ǫ k ) λ k . (A5) Appendix B: Honeycomb bilayer: bosonicBogoliubov transformation
The mean field Hamiltonian of Eq. 15 may be diago-nalized by the matrix, P k = 1 √ − b k b k − b k b k C k , S k , C k , S k , S k , C k , S k , C k , . Here, we have defined b k ≡ β ∗ k / | β k | . We take theother entries to be hyperbolic functions given by C k ,n =cosh κ k ,n and S k ,n = sinh κ k ,n , with n = 1 ,
2. With thisdefinition, this matrix P k satisfies the pseudo-unitaritycondition P k σP † k = σ , where σ = Diag { , , − , − } . Todiagonalize the Hamiltonian matrix M k , we settanh 2 κ k , = β k / ( C − β k );tanh 2 κ k , = − β k / ( C + β k ) . (B1)With this choice, the matrix P k diagonalizes the Hamil-tonian, P † k M k P k = Diag { λ k , , λ k , , λ k , , λ k , } . (B2)where λ k , / are as defined in the main body. We trans-form the triplet operators defined in Eq. 16 into newquasiparticle operators using t k ,A,u t k ,B,u t †− k ,A,u t †− k ,B,u = P k ϑ k , ,u ϑ k , ,u ϑ †− k , ,u ϑ †− k , ,u . (B3) The ϑ operators are the triplon quasiparticles. Com-pared to the square lattice case, the quasiparticle opera-tors have an additional index on account of the sublatticedegree of freedom. We can express our original tripletoperators as follows: t k ,A,u = X f =1 , (cid:16) C k ,f ϑ k ,f,u + S k ,f ϑ †− k ,f,u (cid:17) ,t − k ,B,u = X f =1 , ( − f b ∗ k (cid:16) C k ,f ϑ − k ,f,u + S k ,f ϑ † k ,f,u (cid:17) . (B4)The bilinears defined in Eq. 28 can be evaluated as ρ = 23 N ⊥ X k h t † k ,A,v t k ,B,v i γ k = 16 N ⊥ X k | γ k | (cid:20) − C − | β k | λ k , + C + | β k | λ k , (cid:21) , (B5)∆ = 23 N ⊥ X k h t † k ,A,v t †− k ,B,v i γ k = − N ⊥ X k | γ k | (cid:20) | β k | λ k , + | β k | λ k , (cid:21) . (B6) Appendix C: Spin operator expressions includingquintet terms
Including triplet and quintet operators, the completeexpression for the spin operators on the two layers of thebilayer are S + i,ℓ =1 , = ( − ℓ r S ( S + 1)3 ¯ s ( t i, − − t † i, )+ ( − ℓ r (2 S − S + 3)5 (cid:2) ( t † i, − q i, − − q † i, t i, )+ r
12 ( t † i, q i, − − q † i, t i, ) + r
16 ( t † i, q i, − q † i, t i, − ) (cid:3) + r
12 ( t † i, t i, + t † i, t i, − ) + r
32 ( q † i, q i, + q † i, q i, − )+ q † i, q i, + q † i, − q i, − , (C1)with S − i,ℓ =1 , being its Hermitian conjugate, while S zi,ℓ =1 , = ( − ℓ r S ( S + 1)3 ¯ s ( t i, + t † i, )+ ( − ℓ r (2 S − S + 3)5 (cid:2)r
13 ( t † i, q i, + q † i, t i, )+ 12 ( t † i, q i, + q † i, t i, + t † i, − q i, − + q † i, − t i, − ) (cid:3) + 12 ( t † i, t i, − t † i, − t i, − + q † i, q i, − q † i, − q i, − )+ q † i, q i, − q † i, − q i, − . (C2)2 Appendix D: Square Bilayer: inclusion of quintets
The spin operators with the inclusion of quintets aregiven in Eq. 21 of Ref. 24. Using this reference, we nowgive explicit expressions for ˆ R ttq ( S ). In the main text,we defined ˆ R ttq ( S ) in terms of triplet bilinears ˆ T [ n ] i,i + δ .Here, we give expressions for ˆ T [ n ] i,i + δ in momentum space.We use the Fourier transform convention t i,u ∈{ x,y,z } = 1 √ N ⊥ X k t k ,u e i k . r i . (D1)The operator ˆ T [ n ] i,i + δ is composed of bilinears ofthe form t i,u ( t i + δ,v ± t † i + δ,v ). Using the Fouriertransform, this generic bilinear may be written as(1 /N ⊥ ) P k , p t − k + p ,u ( t k ,u ± t †− k ,u ) e i k .δ e i p . r i .Thus, we may write X δ ˆ T [ n ] i,i + δ = MN ⊥ X k , p ˆ T [ n ] − k + p , k e i p . r i η k , (D2)where η k = P δ e i k · δ = 2(cos k x + cos k y ) and the coef-ficient M = q S ( S +1)(2 S − S +3)30 . The explicit forms ofˆ T [ n ] − k + p , k are:ˆ T [ − − k + p , k = ˜ t − k + p ,x ( t k ,x + t †− k ,x ) − ˜ t − k + p ,y ( t k ,y + t †− k ,y )+ i ˜ t − k + p ,x ( t k ,y + t †− k ,y ) + i ˜ t − k + p ,y ( t k ,x + t †− k ,x )ˆ T [ − − k + p , k = ˜ t − k + p ,z ( t k ,x + t †− k ,x ) + ˜ t − k + p ,x ( t k ,z + t †− k ,z )+ i ˜ t − k + p ,z ( t k ,y + t †− k ,y ) + i ˜ t − k + p ,y ( t k ,z + t †− k ,z )ˆ T [0] − k + p , k = r h − ˜ t − k + p ,x ( t k ,x + t †− k ,x ) − ˜ t − k + p ,y ( t k ,y + t †− k ,y )+2˜ t − k + p ,z ( t k ,z + t †− k ,z ) i ˆ T [ − − k + p , k = − ˜ t − k + p ,z ( t k ,x + t †− k ,x ) − ˜ t − k + p ,x ( t k ,z + t †− k ,z )+ i ˜ t − k + p ,z ( t k ,y + t †− k ,y ) + i ˜ t − k + p ,y ( t k ,z + t †− k ,z )ˆ T [2] − k + p , k = ˜ t − k + p ,x ( t k ,x + t †− k ,x ) − ˜ t − k + p ,y ( t k ,y + t †− k ,y ) − i ˜ t − k + p ,x ( t k ,y + t †− k ,y ) − i ˜ t − k + p ,y ( t k ,x + t †− k ,x )(D3) We have denoted some triplet operators as t and someas ˜ t . For the purposes of the square lattice, this distinc-tion can be ignored. We will use these same expressionsin the context of the honeycomb lattice also. For thehoneycomb case, t and ˜ t operators will act on differentsublattices.The energy correction due to coupling to quintets isgiven in Eq. 36. Using the Fourier transformed expressionin Eq. D2, we rewrite the energy as∆ E S> / = M ¯ s N ⊥ X m = − , ··· , X p E [ m ] p (D4)where p is the momentum of the intermediate state. Thequantity E [ m ] p is given by E [ m ] p = X ν =0 |h ν | P k ˆ T [ n ] − k + p , k η k | i| E − E ν . (D5)Here, ( − p ) is the momentum of the intermediate state | ν i . As described in the Section VII A, the intermediatestates | ν i that contribute have two triplon quasiparticleexcitations and one quintet excitation. Within the tripletsector, an intermediate state with momentum ( − p ) maybe represented as | ν − triplon i = τ † q − p ,u ′ τ †− q ,v ′ | i . (D6)With this parametrization, the sum over intermediatestates | ν i may be written as X ν =0 −→ X q X u ′ ,v ′ ∈{ x,y,z } . (D7)Evaluating the matrix elements using this parametriza-tion of the intermediate state, we find that the energycontribution E [ m ] p is the same from every m-sector, i.e., E [ m ] p = E p for all m . The quantity E p is given by E p = − X q h sinh ( θ q ) η p − q { cosh(2 θ p − q ) + sinh(2 θ p − q } + sinh ( θ p − q ) η q { cosh(2 θ q ) + sinh(2 θ q ) } i ε q − µ + λ − q + λ − p + q (D8) Appendix E: Honeycomb bilayer: inclusion ofquintets
In the main text, we defined ˆ R ttq ( S ) in terms of tripletbilinears ˆ A [ n ] i,i + δ and ˆ B [ n ] i − δ,i . Here, we give expressions for ˆ A [ n ] i,i + δ and ˆ B [ n ] i − δ,i in momentum space. We use theFourier transform convention t i,α ∈{ A,B } ,u ∈{ x,y,z } = 1 p N ⊥ / X k t α, k ,u e i k . r i . (E1)3(i) The terms in ˆ A [ n ] i,i + δ are of the form t i,A,u ( t i + δ,B,v + t † i + δ,B,v ). Using our Fourier transform convention, wemay write X δ ˆ A [ n ] i,i + δ = MN ⊥ / X k , p ˆ A [ n ] − k + p , k e i p · r i γ k , (E2)where γ k = P δ e i k · δ = 1 + e − ik b + e − ik a − ik b and thecoefficient M = q S ( S +1)(2 S − S +3)30 is the same as thatdefined for the square lattice case. The explicit forms ofˆ A [ n ] − k + p , k are the same of those of ˆ T [ n ] − k + p , k given in Eq. D3with the following redefinition:˜ t k ,u ≡ t A, k ,u t k ,u ≡ t B, k ,u (E3)(ii) The terms in ˆ B [ n ] i,i − δ are of the form t i,B,u ( t i − δ,A,v ± t † i − δ,A,v ). Using our Fourier transform convention, wewrite X δ ˆ B [ n ] i,i − δ = MN ⊥ / X k , p ˆ B [ n ] − k + p , k e i p · r i γ − k (E4)Explicit expressions for ˆ B [ n ] − k + p , k are the same as thoseof ˆ T [ n ] − k + p , k given in Eq. D3 but with the following redef-inition: ˜ t k ,u ≡ t B, k ,u t k ,u ≡ t A, k ,u (E5)The quintet energy correction on the honeycomb latticemay be rewritten as∆ E (q) = M ¯ s N ⊥ / X p X m h ( A [ m ] p ) + ( B [ m ] p ) i , (E6) where ( A [ m ] p ) = X ν =0 |h ν | P k ˆ A [ m ] − k + p , k γ k | i| E − E ν ( B [ m ] p ) = X ν =0 |h ν | P k ˆ B [ m ] − k + p , k γ − k | i| E − E ν (E7)The only intermediate states | ν i that contribute to theenergy are states with two triplon quasiparticle excita-tions. A generic intermediate state with momentum ( − p )may be characterized as | ν − triplon i = ϑ †− q ,f,u ϑ † q − p ,g.v | i . (E8)Using this parametrization of a generic state, the sumover intermediate states in Eq. E7 becomes X ν =0 −→ X q X f,g ∈{ , } X u,v ∈{ x,y,z } , (E9)There is a factor of 1 / q ′ = p − q , f ′ = g, g ′ = f ) corresponds to the same stateas ( q , f, g ). Evaluating the necessary overlaps, we findthat the contribution from each m is the same ( A [ m ] p ) =( B [ m ] p ) = E p for m = − , · · · ,
2. The quantity E p isgiven by E p = − X q ,f,g h S q ,f ( − g ( S − p + q ,g + C p − q ,g ) | γ p − q | + S p − q ,g ( − f ( S − q ,f + C q ,f ) | γ q | i ε q − µ + λ − q ,f + λ q − p ,g (E10)By plugging these expressions into Eq.E6, the correction to ground state energy may be computed. T. Nikuni, M. Oshikawa, A. Oosawa, and H. Tanaka, Phys.Rev. Lett. , 5868 (2000). T. M. Rice, Science , 760 (2002). S. Sachdev,
Quantum Phase Transitions , (Cambridge Uni-versity Press, Cambridge, 2001). T. Giamarchi, C. R¨uegg, and O. Tchernyshyov, NaturePhysics , 198 (2008) A. Oosawa, M. Ishii, and H. Tanaka, J. Phys. Condens.Matter , 265 (1999). M. Jaime, V. F. Correa, N. Harrison, C. D. Batista, N.Kawashima, Y. Kazuma, G. A. Jorge, R. Stern, I. Hein-maa, S. A. Zvyagin, Y. Sasago, and K. Uchinokura, Phys.Rev. Lett. , 087203 (2004). M. B. Stone, C. Broholm, D. H. Reich, O. Tchernyshyov, P. Vorderwisch, and N. Harrison, Phys. Rev. Lett. , 257203(2006). M. Kofu, J.-H. Kim, S. Ji, S.-H. Lee, H. Ueda, Y. Qiu,H.-J.Kang, M. A. Green, and Y. Ueda, Phys. Rev. Lett. , 037206 (2009). A. J. Millis and H. Monien, Phys. Rev. Lett. , 2810(1993). A. V. Chubukov and D. K. Morr, Phys. Rev. B , 3521(1995). A. W. Sandvik and D. J. Scalapino, Phys. Rev. Lett. ,2777 (1994). T. Dodds, B.-J. Yang, and Y.-B. Kim, Phys. Rev. B ,054412 (2010). M. B. Stone, M. D. Lumsden, S. Chang, E. C. Samulon, C.D. Batista, and I. R. Fisher, Phys. Rev. Lett. , 237201(2008). H. Liao and T. Li, arXiv:1102.1819 (unpublished). A. Abendschein, and S. Capponi, Phys. Rev. B , 064413(2007). A. Collins and C. J. Hamer, Phys. Rev. B , 054419(2008). Y. Matsushita, M. P. Gelfand, and C. Ishii, J. Phys. Soc.Jpn. , 247 (1999). D.-k. Yu, Q. Gu, H.-T. Wang, and Jue-Lian Shen, Phys.Rev. B , 111 (1999. T. Roscilde and S. Haas, Phys. Rev. Lett. , 207206(2005). K.-K. Ng, F. C. Zhang, and M. Ma, Phys. Rev. B ,12196 (1996) M. P. Gelfand, Z. Weihong, C. J. Hamer, and J. Oitmaa,Phys. Rev. B , 392 (1998). K. Gregor, and O. I. Motrunich, Phys. Rev. B , 174404(2007). A. W. Sandvik, Phys. Rev. B , R14157 (1999); O. F.Sylju˚asen and A. W. Sandvik, Phys. Rev. E , 046701(2002). B. Kumar, Phys. Rev. B , 054404 (2010). O. Smirnova, M. Azuma, N. Kumada, Y. Kusano, M. Mat-suda, Y. Shimakawa, T. Takei, Y. Yonesaki, and N. Kino- mura, J. Am. Chem. Soc., , 8313 (2009). F. Alet, S. Wessel, and M. Troyer, Phys. Rev. E , 036706(2005). E. L. Pollock and D. M. Ceperley, Phys. Rev. B , 8343(1987). L. Wang, K. S. D. Beach, and A. W. Sandvik, Phys. Rev.B , 014431 (2006). M. Campostrini, M. Hasenbusch, A. Pelissetto, P. Rossi,and E. Vicari, Phys. Rev. B , 144520 (2002). S. Sachdev and R. N. Bhatt, Phys. Rev. B , 9323 (1990). A. V. Chubukov, Phys. Rev. B , 3337 (1991). H.-T. Wang, H. Q. Lin, and J.-L. Shen, Phys. Rev. B ,4019 (2000). L. Fritz, R. L. Doretto, S. Wessel, S. Wenzel, S. Burdin,and M. Vojta, Phys. Rev. B , 174416 (2011). S. Gopalan, T. M. Rice, and M. Sigrist, Phys. Rev. B ,8901 (1994). H. C. Kandpal and J. van den Brink, Phys. Rev. B ,140412 (2011). R. Ganesh, D. N. Sheng, Young-June Kim and A.Paramekanti, Phys. Rev. B , 144414 (2011). H. Mosadeq, F. Shahbazi, and S. A. Jafari, J. Phys.: Cond.Matt. , 226006 (2011). A. F. Albuquerque, D. Schwandt, B. Het´enyi, S. Capponi,M. Mambrini, A. M. La¨uchli, Phys. Rev. B , 024406(2011). B. K. Clark, and D. A. Abanin, and S. L. Sondhi, Phys.Rev. Lett. , 087204 (2011). F. Wang, Phys. Rev. B , 024419 (2010). Y.-M. Lu, and Y. Ran, Phys. Rev. B , 024420 (2011). S. Okumura, H. Kawamura, T. Okubo, and Y. Motome, J.Phys. Soc. Japan , 114705 (2010). A. Mulder, R. Ganesh, L. Capriotti, and A. Paramekanti,Phys. Rev. B , 214419 (2010). E. C. Samulon, Y. Kohama, R. D. McDonald, M. C.Shapiro, K. A. Al-Hassanieh, C. D. Batista, M. Jaime, I.R. Fisher, Phys. Rev. Lett.103