A novel nonlinear spin wave theory for the spin 1/2 antiferromagnetic Heisenberg model on a triangular lattice
aa r X i v : . [ c ond - m a t . s t r- e l ] O c t APS/123-QED
A novel nonlinear spin wave theory for the spin antiferromagnetic Heisenberg modelon a triangular lattice Lihua Wang and Sung Gong Chung , Computational Condensed-Matter Physics Laboratory, RIKEN, Wako, Saitama 351-0198, Japan Department of Physics and Nanotechnology Research and Computation Center,Western Michigan University, Kalamazoo, Michigan 49008, USA and Asia Pacific Center for Theoretical Physics, Pohang, Gyeonbuk 790-784, South Korea (Dated: July 3, 2018)We extend the nonlinear spin wave theory (NLSWT) for the spin 1/2 antiferromagnetic Heisenbergmodel on a triangular lattice (TAFHM). This novel NLSWT considers the corrections one orderhigher in 1 /S than the linear spin wave theory (LSWT). It also distinguishes in which circumstancethe negative energy excitation, the sign of the breakdown of LSWT, shall be renormalized to bepositive both by a boson normal ordering and a self-consistent iteration. We draw a phase diagramby testing the stability of various magnetic orders for different parameters. In particular, theincommensurate configuration is found unstable by our study. The new phase transition point(PTP) of the collinear configuration agrees well with various previous studies. PACS numbers: 75.30.Ds , 75.10.Jm , 75.40.Mg, 75.30.Kz
I. INTRODUCTION S = 1 / CuCl , κ -(BEDT-TTF) Cu (CN) , and EtMe Sb[Pd(dmit) ] for whichexperiments have indicated a magnetically disorderedMott insulator, i.e., a possible spin liquid state, at lowtemperatures. Various theories are applied to understandhow the geometric frustration and low dimensionality af-fect the physics. The diagonalization method is exactbut limited in the lattice size, while the finite size ef-fect is so strong to even change the phase diagram .There are also nearly-exact numerical simulations suchas DMRG . They have encountered problems with over-whelmingly large entanglement when the dimensionalityincreases to 2D, although our own numerical work usingthe Entanglement Perturbation Theory(EPT), moves asmall step further in this direction. The numerical sim-ulation is powerful as we already see very rich detailsof the physics which is believed to be closely related tothe geometric frustration, low dimensionality and strongquantum fluctuation. An example is the phase changedue to the finite-size effect mentioned above.However, currently, the most efficient theoreticalframework to understand the model in the thermody-namic limit is the spin wave theory (SWT), which so farunfortunately only gives approximate description. Thereare two well known questions on the validity of SWT, i.e.,if the assumed magnetic order is physical and if its ap-proximation is sufficient when the first answer is positive,because the conversion to the bosonic Hamiltonian in anexpansion of Holstein-Primakoff transformation restoresthe spin statistics only either in large S limit or whenhigher order corrections are properly taken into account.To get an insight into the answers, one usually examinesthe excitation spectrum. First, there must be Goldstonemodes in the K-space to reflect the symmetry of both the lattice and the assumed order. It is automatically satis-fied in most SWT including LSWT because the unitarytransformation, that is the starting point of SWT, natu-rally inherits those symmetries. Special cares, however,should be taken if it is involved with a self-consistentiteration or when there is a singularity cancelling issuein the perturbation. Second, the excitation spectrummust be positive. The emergence of the negative energyis a sign of the overwhelming quantum fluctuation overthe assumed classical configuration, implying the break-down of SWT. Furthermore the interplay of the abovetwo factors can sometimes lead to subtle situations. Par-ticularly for a spin 1 / that yield a le-gal spectrum, but are unphysical at all. In contrast, thespectrum from LSWT has areas of negative energy forthe collinear configuration when 1 . < µ ≡ J ′ /J < . µ being the anisotropy defined below. But a variety ofprevious studies suggest that the collinear order isrealized in this regime. It is believed that the instabilityof LSWT manifests the insufficiency of approximation.Therefore one of the objectives of this paper is to seeif a NLSWT can fill this gap, with the negative energyexcitation being considered in a novel way.Now let us clarify the model. J ′ mentioned above is theinter-chain spin exchange strength and J intra-chain. SeeFig.1 for the illustration. The nearest neighbors alongthe three orientations, 0 o , 60 o and 120 o , are defined as( i, j ) δ . δ and δ respectively correspond to the primitivetranslation vectors ~τ and ~τ shown in Fig.1. ~x directionis along ~τ . The Hamiltonian is H = X δ J δ X ( i,j ) δ (cid:0) S xi S xj + S yi S yj + S zi S zj (cid:1) (1)Here we consider the anisotropy J = J , J = J = J ′ .The expected magnetic orders are signaled by the con-tour peaks at wave number ~Q ’s in the spin structure fac- τ → τ → JJ' J'
FIG. 1: Spin 1/2 TAFHM with the anisotropic nearest neigh-bor couplings J and J ′ ( J ′ /J ≥
0) indicated in the figure.The primitive translation vectors are denoted by τ and τ .The periodic boundary conditions are assumed. (a) µ = 1 . µ = 1 . -1 0 1-101 / π y q q x /π -1 0 1-101 / π y q q x /π FIG. 2: Contour plots of the pin structure factor by EPT for a6 × µ = 1 . µ = 1 .
3. (a) refers to the 120 o spiral spin configuration, (b)collinear. tor, Fig.2, from our numerical simulation . Typically (a) ~Q = (cid:0) π , (cid:1) , etc, for the 120 o spiral configuration and(b) ~Q = (cid:16) , π √ (cid:17) , etc, collinear. They are used to rotatethe laboratory coordinate frame to unify the direction ofthe local spin in SWT.Since many works have been done to formulateSWT for the spin 1/2 TAFHM, we simply start from theunitarily transformed Hamiltonian H ′ = X δ J δ X ( i,j ) δ (cid:2) S yi S yj + cos (cid:16) ~Q · ~δ (cid:17) (cid:0) S zi S zj + S xi S xj (cid:1) + sin (cid:16) ~Q · ~δ (cid:17) (cid:0) S zi S xj − S xi S zj (cid:1)(cid:3) (2)It is the Holstein-Primakoff transformation that trans-forms the above spin Hamiltonian to the bosonic operator(magnon) version. In practice it is approximated by ex-pansion in 1 /S . For example, LSWT is in the order of O ((1 /S ) − / ). The expansion restores the spin statisticsin large S limit even in LSWT. However, for finite S , it needs infinite order to do so. In particular for S = 1 / /S . It takes the following expansion S x = √ S (cid:18) a † + a − a † a † a + a † aa S (cid:19) S y = i √ S (cid:18) a † − a − a † a † a − a † aa S (cid:19) (3) S z = S − a † a The Hamiltonian is rewritten in the K-space after theFourier transformation as H ′ = H + H + H + H (4)where H is a constant contribution to the ground stateenergy (GSE) E GS ≡ H = LJ cosQ x + LJ ′ cos Q x cos √ Q y L being the lattice size. On the other hand, H isthe quantum fluctuation bilinear with respect to magnonoperators. It is H = X k (cid:20) A k a † k a k − B k (cid:16) a † k a †− k + a − k a k (cid:17)(cid:21) (6)with A k = − JcosQ x − J ′ cos Q x cos √ Q y J cosQ x + 12 cosk x + J ′ cos Q x cos √ Q y ! cos k x cos √ k y B k = − J cosQ x − cosk x − J ′ cos Q x cos √ Q y − ! cos k x cos √ k y H and H , which arethe only terms considered by LSWT. H and H are leftfor the later discussion. There is a handy transformationnamed after Bobolyubov to diagonalize H . Its usualdefinition is as follows a k = u k b k + v k b †− k a † k = u k b † k + v k b − k (9)The coefficients, u k and v k , are the roots of the followingequations u k + v k = S ( A k ) A k p | A K − B K | u k v k = S ( A k ) B k p | A K − B K | (10)Note that S ( A K ) ≡ Sign ( A k ) is used in (10). Since (10)is in second order, one usually complements it with thefollowing convention u k − v k = 1 (11)to preserve the same commutation relationship ofmagnon. Then H , after the Bogolyubov transformation,reads as H = − X k A k X k ǫ k (cid:18) b † k b k + 12 (cid:19) (12)where ǫ k = S ( A k ) S (cid:0) A k − B k (cid:1) q | A k − B k | (13)defines the excitation spectrum. Since the ground stateis the vacuum of Bogolyubov bosons, H contributes toGSE as follows E ′ GS = 12 X k ( ǫ k − A k ) (14)Therefore GSE calculated within LSWT is E GS = E GS + E ′ GS (15)Now look at the property of ǫ k . It can be easily shownthat A k = B k at ~k = (0 , A k = − B k at ~k = ± ~Q .Therefore they are always the zero-energy points, the so-called Goldstone modes. It is the characteristics of asystem which spontaneously breaks the symmetry in thethermodynamic limit. In addition, the fisrt derivative ofa linear dispersion at the Goldstone modes defines themagnon velocity, indicating that the excitation is indeedmassless, another characteristics of the spontaneous sym-metry breaking. It can also be shown that there mightbe areas in K-space that have the negative energy exci-tation as µ varies. For example, a collinear configurationwhen ~Q = (cid:0) , π/ √ (cid:1) has A k = J ( − µ + cosk x ) B k = J µcos k x cos √ k y ! (16) A k − B k is then expanded around ~k = ~Q as follows A k − B k = J (h ( µ − − i k x + 3 µ (cid:18) k y − π √ (cid:19) ) (17)from which it is seen that the stable region of the collinearstates is 2 . ≤ µ ≤ ∞ in LSWT. (Both A k and B k arepositive all over the K -space). LSWT starts to yield thenegative energy excitation in the area surrounding thoseGoldstone modes when µ decreases below 2 .
0. It is obvi-ously unphysical, meaning that LSWT fails to account forthe collinear spin configuration in this parameter regime.On the other hand, a variety of studies, either theoret-ical or numerical including our own work, suggest thatthe collinear configuration is realized when 1 . < µ ≤ ∞ .Therefore it is desirable to investigate it in a way dif-ferent from LSWT. To this end, we here develop a NL-SWT which extends the Bogolyubov transformation tohandle negative energy modes. We also introduce a self-consistent iteration scheme. II. A NOVEL NONLINEAR SPIN WAVETHEORY
First of all, it is possible that the negative energy exci-tation exists even for SWT in infinite order if the startingspin configuration is unphysical. The emergence of neg-ative energy excitation could also be merely due to aninsufficient approximation even if that spin configurationphysically exists. The example here is the negative en-ergy excitation in LSWT for the collinear configurationwhen µ < .
0, which can be further divided into µ ≤ µ and µ < µ < . µ is the physical PTP of thecollinear configuration. The negative energy excitationmight exist in the former regime for infinite order SWT,showing a sign of instability of such collinear spin order.While in the latter regime, one should be able to get anall positive spectrum in SWT with a sufficient approxi-mation, see below.Our idea of handling the negative energy modes is sim-ply to proceed like in the Dirac theory of positron .We apply the particle-hole (anti-particle) transformation, b → b † , with a boson normal-ordering, thereby formallyfilling up those negative energy modes and regard thefilled states as a new physical vacuum. Since this for-mally involves an infinite negative energy, it is a sort ofsingular perturbation theory. When combined with meanfield treatment of non-linear terms, we arrive at an iter-ative scheme of normal ordering and mean field approxi-mations to eliminate negative energy modes and improveon the physical vacuum which consists of zero-point mo-tion of all positive energy spectra. Or equivalently, wecan make a singular choice, u k − v k = − k , which leadsto the particle-hole transformation, h b k , b † k i = (cid:0) u k − v k (cid:1) h a k , a † k i = − H = 12 X k ( ǫ k − A k ) + X k C ∗ k ǫ k b † k b k (20)where C ∗ = 1 for positive energy modes and − H due to its relative simplicity and then H . After theformulation, we analyze the collinear spin configurationfirst. It is a self-consistent iterative treatment of H alone, since H vanishes for the collinear configuration.This is followed by the overall correction combining H and H together for a general non-collinear spin config-uration and anisotropy. A. Quartic term
The quartic term written in the magnon operator readsas H = H + H + H (21)where H = X δ J δ X ( i,j ) δ cos (cid:16) ~Q · ~δ (cid:17) a † i a i a † j a j (22) H = X δ J δ X ( i,j ) δ − cos (cid:16) ~Q · ~δ (cid:17) (cid:16) a † i a † j a † j a j + a i a † j a j a j + a † i a † i a i a † j + a † i a i a i a j (cid:17) (23)etc. After the Bogolyubov transformation with eitherregular boson or new boson with a normal ordering, itcan be written as H = δE + δH + ¯ H (24)The first term is the Hartree-Fock correction to theground state energy and the second the magnon self-energy. Since the last term describes two-particle scatter-ing processe, it only yields higher order 1 /S correctionscompared to δE and δH . So we neglect it here. Beforederiving the explicit form of δE and δH , we apply theHartree-Fock decoupling to (22). Due to Wick’s theorem, one needs to define the following mean-field parameters n ≡ D a † i a i E = X k C ∗ k v k (25) m ≡ D a † i a j E = X k C ∗ k γ k v k (26)∆ ≡ h a i a j i = X k C ∗ k γ k u k v k (27)Ω ≡ (cid:10) a i (cid:11) = X k C ∗ k u k v k (28)To include the effect of anisotropy, one also needs to de-fine m δ and ∆ δ respectively as follows m ≡ D a † i a j E δ = X k C ∗ k cosk x v k ∆ ≡ h a i a j i δ = X k C ∗ k cos k x √ k y ! u k v k (29)etc, where u k and v k are defined either to give the reg-ular bosonic commutation or the new commutation forthe Bogolyubov transformation depending on the sign ofthe excitation energy. One should notice that in the ex-tended stable parameter region all the excitation wouldbe positive after the self-consistent iteration converges.Therefore the mixture of bosons and new bosons (theydon’t interact because they are defined for different wavenumbers) will evolve during the iteration. in the con-verged state, only regular bosons survive. The mean-fieldparameters are functional of both ~Q and ~k . On the otherhand, γ k only reflects the average of the nearest neigh-bors in three orientations of a triangular lattice. It isdefined as γ k ≡ X δ e i~k · ~δ = 13 cosk x + 2 cos k x cos √ k y ! (30)Now let us look at the following term which is a constantcorrection to the ground state energy, δE = X δ J δ " cos (cid:16) ~Q · ~δ (cid:17) (cid:0) n + m δ + ∆ δ (cid:1) + 1 − cos (cid:16) ~Q · ~δ (cid:17) n ∆ δ + m δ Ω) − cos (cid:16) ~Q · ~δ (cid:17) nm δ + ∆ δ Ω) (31)It emerges after the normal-ordering of the quartic term.The left-over two-operator terms are mean-field decou-pled, multiplied by the mean-field parameters. For ex-ample, a † i a i a † j a j → D a † i a i E a † j a j + D a † j a j E a † i a i + D a † i a † j E a i a j + h a i a j i a † i a † j + D a † i a j E a † j a i + D a † j a i E a † i a j (32)The correction to the quantum fluctuation of LSW fromthe quartic term reads in the K -space as follows δH = X k δA k a † k a k − δB k (cid:16) a † k a †− k + a k a − k (cid:17) (33)where δA k (cid:18) cosQ x m + 1 − cosQ x − cosQ x n (cid:19) cosk x + µ (cid:18) f Q m ′ + 1 − f Q − (1 + f Q ) n (cid:19) f k + 1 − cosQ x + µ (1 − f Q ) ∆ ′ − cosQ x m − µ (1 + f Q ) m ′ + ( cosQ x + 2 µf Q ) n (34) δB k (cid:18) cosQ x − n + cosQ x + 14 Ω − cosQ x ∆ (cid:19) cosk x + µ (cid:20) ( f Q − n + 1 + f Q − f Q ∆ ′ (cid:21) f k + cosQ x − m + cosQ x + 14 ∆ + µ (cid:18) f Q − m ′ + f Q + 12 ∆ ′ (cid:19) (35)with the auxiliary mean-filed parameters being definedas follows m ′ = X k C ∗ k f k v k ∆ ′ = X k C ∗ k f k u k v k (36) f k = cos k x cos √ k y B. Cubic term
The cubic term H in (4) reads after the Primakofftransformation as H = X ( i,j ) δ √ SJ δ sin (cid:16) ~Q · ~δ (cid:17) h a † i a i (cid:16) a † j + a j (cid:17) − a † j a j (cid:16) a † i + a i (cid:17)i (37) Under the Bogolyubov transformation, it becomes H = X kq h Γ ( q, k − q, k ) b † q b † q − k b k + Γ ( q, − k − q, k ) b † q b †− q − k b † k + h.c. i (38)whereΓ (1 , ,
3) = iJ √ S
2! [ Z ( u + v ) ( u u + v v )+ Z ( u + v ) ( u u + v v )+ Z ( u + v ) ( u v + v u )] (39)Γ (1 , ,
3) = iJ √ S
3! [ Z ( u + v ) ( u v + v u )+ Z ( u + v ) ( u v + v u )+ Z ( u + v ) ( u v + v u )] (40)Above we wrote the expressions symmetric to the wavenumbers. 2! in Γ is from the the permutation of 1 ↔ b † b † b ( b b b † ). And 3! in Γ is from the permuta-tion among all three subscripts of b † b † b † ( b b b ). Thissymmetric form is somehow more suitable for the numer-ical integration over the first BZ, especially when there iscancellation of singular terms. Similar to the one from but (38) now includes the anisotropy as follows Z ( Q, k, µ ) = sinQ x sink x + 2 µ ( h Q h k + g Q g k ) (41) h k = sin k x cos √ k y g k = cos k x sin √ k y H eff = H + H (42)where H either refers to the LSW term H or any unper-turbed Hamiltonian. If the anisotropy is set to vanish,(39) coincides with the ones in . Following the stan-dard diagram perturbation method (See Appendix fordetails), the lowest order contributions are X a ( k, ω ) = 2 X q | Γ ( q, k ) | ω − ǫ q − ǫ k − q + i X b ( k, ω ) = − X q | Γ ( q, k ) | ω + ǫ q + ǫ k + q − i X c ( k, ω ) = 6 X q Γ ( q, k ) Γ ∗ ( q, − k ) ω + ǫ q + ǫ k + q − i X d ( k, ω ) = − X q Γ ( q, − k ) Γ ∗ ( q, k ) ω − ǫ q − ǫ k − q + i O (cid:0) S (cid:1) as the normal one’s. The new poles of themagnon Green’s function will then only reflect the quar-tic term and the normal self-energies from the cubic term.The real part of the poles determines the renormalizedspectrum while the imaginary part determines the damp-ing rate. Explicitly, the new pole is as follows ǫ k = ¯ ǫ k + iχ k = ǫ k + X q | Γ ( q, k ) | ǫ q + ǫ k − q − ǫ k − i | Γ ( q, k ) | ǫ q + ǫ k + q + ǫ k ! (47)Since Γ , are expressed in terms of u k and v k , it is de-sirable to write out Γ , explicitly in order to avoid theuncertainty of the sign of u k and v k . For example,2Γ (1 , , SJ = − X i =1 T i (48)with T Z = ( u + v ) ( u u + v v ) = S ( A ) ( A + B )2 p | A − B | S ( A A ) ( A A + B B ) p | ( A − B ) ( A − B ) | + C ∗ C ∗ ! (49)and2 T Z Z =4 ( u + v ) ( u u + v v ) ( u + v ) ( u u + v v )= C ∗ C ∗ S ( A ) ( A + B ) p | A − B | + C ∗ C ∗ S ( A ) ( A + B ) p | A − B | + C ∗ C ∗ S ( A ) ( A − B ) p | A − B | + S ( A A A ) ( A + B ) ( A + B ) ( A + B ) p | ( A − B ) ( A − B ) ( A − B ) | (50)etc. Γ is formulated likewise. Note that the sign of theexcitation energy of the unperturbed Hamiltonian is in-cluded. In the case where the unperturbed spectrum ispositive, it has the same expression as after the usual Bo-golyubov transformation. However, (50) should be usedwhenever one considers NLSWT treatment involving aproblematic spectrum from LSWT with a small area ofnegative energy. More discussion is given in the next twosections. III. SELF-CONSISTENT NSWT ANALYSIS OFCOLLINEAR SPIN CONFIGURATION
The purpose of this paper is to investigate the phasediagram of TAFHM by NLSWT, taking into account thenegative energy excitation issue with higher order cor-rections. For the collinear spin configuration, we need tostudy the spectrum renormalization due to H alone. Inthe last section, we extend the Bogolyubove transforma-tion to reformulate δH in (33), the contribution from H to the quantum fluctuation H in LSWT. For thecollinear spin configuration, the coefficients in (33) are δA k = (2 − µ ) n + 4 µ ∆ ′ − m + (2 m − n ) cosk x + µ (2Ω − m ′ ) f k (51) δB k = (Ω − ) cosk x + µ (4∆ ′ − n ) f k + ∆ − µm ′ (52)The correction from H to the ground state energy, (31),is written explicitly for the collinear configuration as well δE = (1 − µ ) (cid:0) n + m + ∆ (cid:1) + 2 µ (2 n ∆ + m Ω) − (2 nm + ∆Ω) (53)Now we look at the spectrum renormalization due to δH . One can obtain it by¯ ǫ k = ǫ k + S ( A k ) A k δA k − B k δB k | ǫ k | (54)or by ¯ ǫ k = C ∗ ′ k r(cid:12)(cid:12)(cid:12) ( A K + δA k ) − ( B K + δB k ) (cid:12)(cid:12)(cid:12) (55)where C ∗ ′ k is determined in the same way as C ∗ k . One im-mediately sees that when the correction δA k and δB k aresmall compared with A k and B k (54) is almost equiva-lent to (55) except that (54) has singularity issue at someGoldstone modes where the sign of δB k /δA k is oppositeto that of B k /A k . Physically this singularity will be can-celled either in the summation over the wave numbersin the first BZ or by the other source such as the cubicterm in the same order as we will discuss later. However,it is challenging to numerically carry it out. Anotherdrawback of (54) is that the mean-field quantities do notcontain the information of the higher order correction.Therefore we prefer to use (55). In fact, the first way canbe formulated from the Green’s function using H as theunperturbed Hamiltonian. The second way can be ob-tained by absorbing δH into the unperturbed Hamilto-nian. We believe that when the correction is small, thesetwo approximations shall both yield satisfactory results.However so far they are only a half-cooking recipe, be-cause even if they can obtain positive spectrum, neitherguarantees the non-recurrance of negative energy prob-lem. The situation would be more sound if a convergencecan be provided by a self-consistent iteration. Below we -0.70-0.65-0.60-0.55-0.50-0.45-0.40 =2.0 E G S / LSWT NLSWT WITHOUT S-I NLSWT WITH S-I VMC,EPT =1.3
FIG. 3: LSWT and NLSWT results for the collinear spinconfiguration on a square lattice by approaching µ to ∞ andsetting ~Q = (cid:16) , π √ (cid:17) . illustrate the idea.First we start the self-consistent iteration from theground state of LSWT to get the mean-field parame-ters. If it has the negative energy excitation for somewave numbers, it is not the lowest energy state but ameta-stable state. When we calculate the mean-field pa-rameters, the new commutation (18) should be used tocollect the contribution of those wave numbers associatedwith the negative energy excitation. At the end of eachiteration one uses the obtained mean-field parameter toupdate the excitation spectrum. In other words, one re-assesses which excitation is regular boson or the new bo-son in the first BZ. In the next iteration, the groundstate is chosen to be the vacuum of those updated bo-son and new boson. Still it is a meta-stable state if thenegative excitation does not vanish in the whole spec-trum. This iteration is carried on until convergence. Inthe end, there are two possibilities. First, the excitationspectrum becomes all positive, and the ground state hasthe lowest energy, a physical ground state. Second, iter-ation is converged but there still remain some negativeenergy excitations, signaling the instability of assumedmagnetic order.We perform our calculation on a 1800 × ~Q = (cid:16) , π √ (cid:17) and µ → ∞ to check the well known re-sult of a square lattice with a collinear spin configuration.In this case, there is no negative energy excitation issueeven in LSWT. Indeed, in Fig.3 one of our NLSWT calcu-lations, which uses (54) and doesn’t take a self-consistentiteration process, agrees well with . We also see inFig3 that the energy after the self-consitent iteration con-verged is obviously further from the believed ground stateenergy shown by star due to variational Monte Carlo (a) µ = 1 . µ = 1 . Kx K y -4 -2 0 2 4 -202 2.52.01.51.00.50.0 Kx K y -4 -2 0 2 4 -202 FIG. 4: Comparison between the contour plots of the exci-tation spectrum in the first BZ for (a) µ = 1 . . that is tested by our own nemerical simulation . Sinceself-consistency means to take higher order correlationsinto account, NLSWT with self-consistent iteration is ex-pected to give a better ground state energy than withoutself-consistent iteration. That it is not the case reflects afact that SWT is not variational and hence an approachto the correct value need not be monotonic.On the other hand, the self-consistent iteration worksbetter in the intermediate parameter regime in the sensethat it improves on the excitation spectrum in 1 . ≤ µ ≤ .
00 and that it avoids an artificial oscillation ofthe NLSWT results without a self-consistent iteration.Furthermore, the discontinuity at µ ≈ . µ ≥ . µ decreases below 1 . µ = 2 .
0. In contrast, it is µ ≈ . µ de-creases below their PTPs, the softening of the spin waveat Γ-point deserves an emphasis. The strongest softeninghappens when the sound velocity is calculated along thedirection from Γ-point to X-point. For example, it is zeroat µ = 2 . µ in the range of 1 . ≤ µ ≤ . -2024 X k / J k =1.2 NLSWT =1.3 NLSWT =1.3 LSWTM 024 k / J k =1.9 LSWT =2.0 LSWT =2.0 NLSWTM X FIG. 5: Comparison between the spectra of LSWT and NL-SWT around the phase transition point (PTP) at (a) µ = 1 . . PTP at Each Step PTP of LSWT
Iteration Step
FIG. 6: The star (blue online) show PTP of LSWT. It movesto µ = 1 .
52 immediately after one iteration. Eventually itconverges to µ ≈ . is an important feature, not only because it reveals thatthe new commutation relation renormalizes the spectrumsystematically (In fact, without a correct treatment of thenegative energy excitation we found no way to renormal-ize the spectrum to the physically meaningful one whilethe correct symmetry, i.e., the Goldstone mode, is in-tact), but also because this nice tendency yields a usefulinformation about the phase transition even when it doesnot converge. IV. NLSWT ANALYSIS OF NON-COLLINEARSPIN CONFIGURATION
We also want to test our method for non-collinear spinconfigurations and anisotropy by applying it to the limit-ing case which can be readily compared with the existingcalculations. We discuss in Sec.IV A the commensurate (a) (b) =2 /L = /L [11] I (150/L) I FIG. 7: The finite-size scaling of the contribution to theground state energy from (a) CT and (b) QT. This scalingis for the 120 o spiral configuration on an isotropic lattice spiral order, ~Q = (4 π/ ,
0) and µ = 1. It is a 120 o threesub-lattice spin configuration on the isotropic lattice. InSec.IV B we extend the discussion to the incommensuratespiral order. A. o Commensurate Spiral Order
Note that the contributions from the cubic term (CT)and from the quartic term (QT) both have the singularityat ~k = ± ~Q . From Fig.7 one sees that their singularitiesexactly cancel out only at the infinite limit. Since thenumerical calculation can be done only on a finite lat-tice, a finite-size scaling is necessary to extrapolate theprecise corrections. For example we show two quantitiesthat are related to the contribution to GSE, I I I ≡ − δE / J with δE being defined in (31) and I ≡ − δE / J with δE being the contribution to GSE from CT. I I | Γ ( q, k ) | ǫ q + ǫ k + q + ǫ k (56)One sees from Fig.7 that I I
4. Note that our direct numerial simulation is muchheavier compared to the finite-size extrapolation usingclusters with different aspect ratios , and gives a clearconvergence as seen in Fig.7. The extrapolated values of I I denoted by blue stars. In addi-tion, it is well known that in a direct numerical simulationa small number δ , comparable with the increment of nu-merical integration in the K space, should be inserted inthe denominator of any fractional number. It is to pre-vent the overflow in the numerical simulation. Thereforefinite-size scaling using different δ is presented in the fig-ure. Smaller δ gives a faster convergence. The standardchoice δ = π/L is also good as seen in the figure. k k
300 600 1200 2400 4800 9600 19200 38400 76800 K FIG. 8: The finite-size scaling of the spectrum for the 120 o spiral configuration on an isotropic lattice. The slice from Γpoint, ~k = (0 ,
0) to K point, ~k = ~Q = (4 π/ , We then use δ = π/L to perform the numerical simu-lation for the spectrum. A finite-size scaling is shown inFig.8. It is clear that the singularities from CT and QTat K -point indeed exactly cancel out at the thermody-namic limit but it scales very slowly. The spectra fromtop to bottom in the figure account for different latticesizes, from L = 300 to 76800. The extrapolated spectrumagrees with .It is also notable that there exists discontinuity in thespectrum, shown as X -point in Fig.7. It could be curedby some techniques. For example the damping can beabsorbed into the Dyson equation to further renormalizethe spectrum like in the off-shell approximation . How-ever we don’t go in this direction for two reasons. First,the 120 o spiral order when µ = 1 . µ varies away from 1 . µ . B. Incommensurate Sprial Order
There exist incommensurate spiral states in LSWTwhen µ varies away from 1 .
0. The relationship be-tween µ and the ordering wave number is: ~Q = (cid:0) cos − ( − µ/ , (cid:1) when 0 . < µ < . ~Q = (cid:0) π − cos − ( − µ/ , π/ √ (cid:1) when 1 . < µ < .
0. Notethat the splitting into branches is just due to the foldinginto the first BZ of an originally continuous line. A sim-ple calculation shows that this relationship is identicalto the one in the classical picture except that it extends (a) (b) k k K 01 K k k FIG. 9: The finite-size scaling of the spectrum for incommen-surate spiral configurations. A slice is shown, from Γ point to K -point of (a) (7 π/ ,
0) and (b) (cid:0) π/ , π/ √ (cid:1) into 0 < µ < .
27. The difference can be explained if thereduction, the difference between the magnetization andthe classical value (1 / < µ < .
27 in LSWT, indicating its breakdownin this regime and that the state is disordered.We study these incommensurate spiral configurationsby NLSWT in this session. A few wave numbers arechosen to carry out the numerical simulation. The cor-responding value of µ resides on both sides of 1 .
0. Andonly two typical spectra are plotted. They are (a) ofFig.9 when ~Q = (7 π/ ,
0) and (b) (cid:0) π/ , π/ √ (cid:1) . Thespectrum is plotted along the same slice as in Fig.8, withonly the vicinity of K -point, one of the Goldstone modes,being shown.The solid lines show the extrapolated spectrum by thefinite-size scaling. They exhibit both the negative energyexcitation and a discontinuity at K -point. The relation-ship between the distance from the commensurate wavenumber, ∆ Q ≡ Q x − π/ Q x being the unfoldedincommensurate wave number, and ∆, the depth of thenegative energy, i.e., half of the kink since the disconti-nuity is symmetric about 0, is plotted in Fig.10. It islinear near the commensurate spiral wave number. Al-though only a few wave numbers are checked, it seemsin Fig.10 that all the incommensurate wave numbers arecorresponding to spectra with negative energy excitationand discontinuity at K -point. Namely, the incommensu-rate spiral spin configuration is unstable in the NLSWT. C. Summary and Proposal of a Self-ConsistentIteration
From the last two sections we conclude that in thisNLSWT only the commensurate 120 o spiral order is sta-ble among all the spiral configurations. But it would beinteresting to see if special treatments like the off-shellapproximation can cure the discontinuity at the K -point,0 -0.5 0.0 0.50.00.10.20.3 Q FIG. 10: The finite-size scaling of the spectrum for a 120 o spiral configuration on an isotropic lattice. The slice from Γpoint, ~k = (0 ,
0) to K point, ~k = ~Q = (4 π/ , say, at least for some incommensurate spiral wave num-ber. If so, the depth of the negative energy will thendisappear because it is tightly associated with the discon-tinuity. However, this perspective is essentially differentfrom what off-shell approximation does for the commen-surate spiral state where the discontinuity at X -point issmoothed by it while the spectrum at K -point is not al-tered. V. CONCLUSION AND DISCUSSION
We have extended the Bogolyubov transformation byconsidering a new commutation relation for the bosonicquasi-particle which has negative energy. This new trans-formation, equivalent to a boson normal ordering, bridgesa self-consistent iteration from a spectrum with areas ofnegative energy excitation in LSWT to an expected posi-tive energy spectrum in NLSWT. The new stable regime1 . < µ < ∞ for the collinear spin configuration on aspin 1 / . Considering the clear physical picture ofboth the ground state and the excitation spectrum pro-vided by SWT, such improvement of the approximationprecision over the previous SWT studies stimulatesour interest to apply the same scheme to non-collinearspin configurations. Indeed our study shows that amongspiral configurations, only the 120 o commensurate sprialis stable at µ = 1 . ~Q = (cid:0) π/ , π/ √ (cid:1) , namely a four sub-lattice spi-ral state has been seen in a very narrow window around µ = 1 .
18 by our recent Eangtanglement PerturbationTheory study . It implies that due to the approximatenature of SWT, especially when the spin is as small as1 /
2, NLSWT still misses certain details in the phase dia-gram. Moreover, our more elaborate self-consistent itera-tion study for a non-collinear spin configuration does notgive a satisfactory result of an expected broadening of thestable region for 120 o spiral state. A further developmentof NLSWT along this line shall be highly desirable. VI. ACKNOWNLEDGEMENT
Lihua Wang would like to acknowledge the fruitful dis-cussion with Dr. Kazuo Ueda. The numerical calculationof this work utilized the RIKEN Integrated Cluster ofClusters at Advanced Center for Computing and Com-munication, RIKEN.
Appendix A: Self-Energy
In this appendix, we show the formulation of the self-energy in the lowest order for CT, for instance, P a ( k, ω ).The green function is defined as follows G ( k, τ ) = − D T τ ˜ b k ( τ ) b † k (0) E = − D T τ U ( β ) b k ( τ ) b † k (0) E (A1)with U ( β ) = ∞ X n =0 (cid:18) − h (cid:19) n n ! Z β dτ · · · dτ n h T τ H ′ ( τ ) · · · H ′ ( τ n ) i (A2)The perturbative Hamiltonian is H ′ = H + H . Oneimmediately sees that H , if one choose not to absorb itto H , contributes in (A2) as the first order, it will give˜ ǫ k , just a HF approximated term as we discussed before.Now we focus on the self-energy due to H , which arisesfrom the expansion of (A2) in the second order: U (2) ( β ) = 1 h Z β dτ dτ H ( τ ) H ( τ ) (A3)Therefore G (2) ( k, τ ) = − h Z β dτ dτ D H ( τ ) H ( τ ) b k ( τ ) b † k (0) E (A4)Using (38) we have1 H ( τ ) H ( τ ) = X q k Γ ( q , k − q , k ) X q k Γ ∗ ( q , k − q , k ) b † q ( τ ) b † q − k ( τ ) b k ( τ ) b q ( τ ) b q − k ( τ ) b † k ( τ ) + · · · (A5)The first term above will give rise to P a ( k, ω ), Now weshow the formulation step by step. Due to Wick’s the- orem, one of terms in the form of operator-pair to besummed in (A4) is (cid:10) T τ b † q ( τ ) b q ( τ ) (cid:11) D T τ b † k − q ( τ ) b k − q ( τ ) E D T τ b k ( τ ) b † k (0) E D T τ b † k ( τ ) b k ( τ ) E = (cid:10) T τ b q ( τ ) b † q ( τ ) (cid:11) δ q ,q D T τ b k − q ( τ ) b † k − q ( τ ) E δ k − q ,k − q D T τ b k ( τ ) b † k (0) E δ k ,k D T τ b k ( τ ) b † k ( τ ) E δ k,k = G ( q , τ − τ ) δ q ,q G ( k − q , τ − τ ) δ k − q ,k − q G ( k, τ ) δ k ,k G ( k, τ − τ ) δ k,k (A6)Where G is the unperturbed Green function, G ( k, τ ) = − D T τ b k ( τ ) b † k (0) E (A7)Substituting (A7) into (A6) and using the kronig-deltafunction one has G (2)1 ( k, τ ) = − h Z β dτ dτ X q | Γ ( q, k ) | G ( q, τ − τ ) G ( k − q, τ − τ ) G ( k, τ − τ ) G ( k, τ )(A8)The first subscript stands for the lowest order while thesecond for only part of terms. After the Frourier trans-formation, (A8) reads as G (2)1 ( k, iω )= 12 β X q X ω ′ | Γ ( q, k ) | G ( q, iω ′ ) G ( k, iω ) G ( k, iω ) G ( k − q, iω − iω ′ )= G ( k, iω ) X a ( k, ω ) G ( k, iω ) (A9)Thus one gets the implicit expression of part of the firstnormal self-energy X a ( k, ω ) = 12 β X ω ′ X q | Γ ( q, k ) | G ( q, iω ′ ) G ( k − q, iω − iω ′ ) (A10)If one needs its explicit expression, the unperturbedGreen function above must be integrated out with re- spect to q and ω ′ . Take the definition G ( k, iω ) = 1 iω − ω k (A11)And note that the boson distribution function is n B ( z ) = 1 e βz − β X iω → πi I c dzn B ( z ) (A13)Explicitly, one defines a contour integral as follows I = 12 X q | Γ ( q, k ) | (cid:20) πi I c dzn B ( z ) (cid:18) z − iω + ω k − q z − ω q (cid:19)(cid:21) =0 (A14)(A14) can be rewritten as follows I = R + R = 0 (A15)where R = − P a ( k, w ), the residue on poles of n B ,and R is the residue on poles of G ’s. After analyticcontinuation, the expression of P a ( k, w ) is2 X a ( k, ω ) = X q | Γ ( q, k ) | (cid:26) n B ( ω k ) 1 ω q + ω k − q − ω − iδ + [1 − n B ( ω k − q )] 1 ω + iδ − ω k − q − ω q (cid:27) (A16)At the zero temperature, n B = 0, its expression is as fol-lows X a ( k, ω ) = 12 X q | Γ ( q, k ) | ω − ǫ q − ǫ k − q + iδ (A17)The other terms due to Wick’s theorem can be obtainedsimilarly. The final expression of P a ( k, ω ) is X a ( k, ω ) = 2 X q | Γ ( q, k ) | ω − ǫ q − ǫ k − q + iδ (A18) R. Coldea, D. A. Tennant, A. M. Tsvelik, and Z. Tylczyn-ski, Phys. Rev. Lett. , 1335 (2001). Y. Shimizu, K. Miyagawa, K. Kanoda, M. Maesato, andG. Saito, Phys. Rev. Lett. , 107001 (2003). M. Yamashita, N. Nakata, Y. Senshu, M. Nagata, H. M.Yamamoto, R. Kato, T. Shibauchi, and Y. Matsuda, Sci-ence , 1246 (2010). B. Bernu, P. Lecheminant, C. Lhuillier, and L. Pierre,Phys. Rev. B , 10048 (1994). L.Wang, T.Shirakawa, H.Watanabe, and S.Yunoki, J.Phys.: Conf. Ser. , 012131 (2011), l. Wang et al. tobe published (2011). M. Q. Weng, D. N. Sheng, Z. Y. Weng, and R. J. Bursill,Phys. Rev. B , 012407 (2006). T. Holstein and H. Primakoff, Phys. Rev. , 1908 (1940). A. E. Trumper, Phys. Rev. B. , 2987 (1999). S. Yunoki and S. Sorella, Phys Rev. B , 014408 (2006). L. Capriotti, Int. J. Mod. Phys B , 2003 (2003). A. L. Chernyshev and M. E. Zhitomirsky, Phys. Rev. B. , 144416 (2009). J. Merino, R. H. McKenzie, J. B. Marston, and C. H.Chung, J. Phys. Condens. Matter , 2965 (1999). R. P. Feynman, Phys. Rev. , 749 (1949). S. J. Miyake, Prog. Theor. Phys. , 468 (1985). E. Manousakis, Rev. Mod. Phys. (1991). D. A. Huse and V. Elser, Phys. Rev. Lett. , 2531 (1988). P. Horsch and W. V. D. Linden, Z. Phys. B , 181(1988). E. Manousakis, Phys. Rev. B , 4904 (1989). S. R. White and A. L. Chernyshev, Phys. Rev. Lett.99