Magnetic order and spin excitations in layered Heisenberg antiferromagnets with compass-model anisotropies
aa r X i v : . [ c ond - m a t . s t r- e l ] M a r Magnetic order and spin excitations in layered Heisenberg antiferromagnetswith compass-model anisotropies
A.A. Vladimirov a , D. Ihle b and N. M. Plakida a a Joint Institute for Nuclear Research, 141980 Dubna, Russia and b Institut f¨ur Theoretische Physik, Universit¨at Leipzig, D-04109, Leipzig, Germany (Dated: March 9, 2015)The spin-wave excitation spectrum, the magnetization, and the N´eel temperature for the quasi-two-dimensional spin-1/2 antiferromagnetic Heisenberg model with compass-model interaction inthe plane proposed for iridates are calculated in the random phase approximation. The spin-wavespectrum agrees well with data of Lanczos diagonalization. We find that the N´eel temperature isenhanced by the compass-model interaction and is close to the experimental value for Ba IrO . PACS numbers: 74.72.-h, 75.10.-b, 75.40.Gb
Spin-orbital physics in transition-metal oxides hasbeen extensively studied in recent years. A number oftheoretical models was proposed to describe a compli-cated nature of phase transitions induced by compet-ing spin and orbital interactions as originally was con-sidered in Ref. [1]. Whereas the isotropic spin interac-tion can be treated within the conventional Heisenbergmodel, to study the orientation-dependent orbital inter-action the compass model is commonly used. The latterreveals a large degeneracy of ground states resulting in acomplicated phase diagram. In particular, quantum andthermodynamic phase transitions in the two-dimensional(2D) compass model were studied in Refs. [2–4], where afirst-order transition was found for the symmetric com-pass model. A generalized 2D Compass-Heisenberg (CH)model was introduced in Ref. [5], where an important roleof the spin Heisenberg interaction in lifting the high de-generacy of the ground state of the compass model wasstressed. In Ref. [6] a phase diagram of the CH modeland excitations within Lanczos exact diagonalization forfinite clusters on a square lattice were considered in de-tail. In particular, spin-wave excitations and column-flipexcitations in nanoclusters characteristic to the compassmodel were analyzed.A strong relativistic spin-orbital coupling reveals acompass-model type interaction in 5 d transition met-als. In particular, it was shown in Ref. [7], that astrong spin-orbit coupling in such compounds as Sr IrO and Ba IrO results in an effective antiferromagnetic(AF) Heisenberg model for the pseudospins 1 / T N = 230 K in Sr IrO and T N = 240 Kin Ba IrO (see, e.g., [8]). The spin-wave spectrum mea-sured by magnetic resonance inelastic x-ray scattering(RIXS) in Sr IrO shows a dispersion similar to that onein the undoped cuprate La CuO [9].In the present paper we calculate the spin-wave ex-citation spectrum and magnetization for a layered AFHeisenberg model with anisotropic compass-model inter-action in the plane. To take into account the finite-temperature renormalization of the spectrum and to cal-culate the N´eel temperature T N , we employ the equa- tion of motion method for the Green functions (GFs)for spin S = 1 / H = 12 X i,j n J ij S i S j + Γ xij S xi S xj + Γ yij S yi S yj o . (1)Here J ij = J ( δ r j , r i ± a x + δ r j , r i ± a y ) + J z δ r j , r i ± c , where J is the exchange interaction between the nearest neighbors(NN) in the plane with the lattice constants a x = a y = a ,and J z is the coupling between the planes with the dis-tance c . The compass model interaction is given byΓ xij = Γ x δ r j , r i ± a x , Γ yij = Γ y δ r j , r i ± a y . The ab initio many-body quantum chemistry calculations give the fol-lowing parameters for Ba IrO : J = 65 meV, Γ x = Γ y =Γ = 3 . J z & − µ eV [13]. To compare ourresults with the theoretical studies of the 2D CH modelin Ref. [5], we consider also large anisotropic compass-model interactions, Γ x > Γ y > J . In Refs. [11, 12] thespin-wave spectrum was calculated for a similar model(1) in the linear spin-wave theory (LSWT), where in ad-dition to the isotropic exchange interaction J ij betweenthe NN in (1) the next-nearest neighbor (NNN) interac-tion was also taken into account. We consider this moregeneral model in the Appendix.We adopt a two-sublattice ( A, B ) representation forthe AF LRO below the N´eel temperature. Then theHamiltonian (1) with Γ x = Γ y > x, y ) plane – is de-generate. To lift the degeneracy, we assume anisotropiccompass-model interactions Γ x > Γ y > h S xi ⊂ A i = −h S xi ⊂ B i fixed along the x axis. We can con-sider also the limiting case, Γ x = Γ y . The AF LRO canbe described by the AF wave vector Q = ( π/a, π/a, π/c ).It is convenient to write the Hamiltonian (1) in termsof the circular components S ± i = S yi ± iS zi in the form H = 12 X h i,j i n J xij S xi S xj + J yij (cid:2) S + i S − j + S − i S + j (cid:3) + 14 Γ yij (cid:2) S + i S + j + S − i S − j (cid:3) o , (2)where J xij = J ij + Γ xij , J yij = J ij + (1 / yij .To calculate the spin-wave spectrum of transverse spinexcitations, we introduce the retarded two-time commu-tator GFs [14]: G α,βnm ( t − t ′ ) = − iθ ( t − t ′ ) h [ S αn ( t ) , S βm ( t ′ )] i = Z + ∞−∞ dω π e − iω ( t − t ′ ) hh S αn | S βm ii ω , (3)where α, β = ( ± ), and h . . . i is the statistical average.The indexes n, m run over N/ i ( j ) in thesublattice A ( B ).There are four types of the GFs due to the two-sublattice representation for normal and anomalous GFswhich can be written as 4 × G ( ω ) = hh S + i S − i S − j S + j | (cid:16) S − i ′ S + i ′ S + j ′ S − j ′ (cid:17) ii ω . (4)Here the lattice sites i, i ′ refer to the sublattice A whilethe lattice sites j, j ′ refer to the sublattice B .Using equations of motion for spin operators, i ( d/dt ) S ± i ( t ) = [ S ± i , H ] = ∓ P n J xin S ± i S xn ± P n [ J yin S xi S ± n + (1 / yin S xi S ∓ n ] , we obtain a sys-tem of equations for the matrix components of theGF (4). In particular, ω hh S + i | S − i ′ ii ω = 2 h S xi i δ i,i ′ − X n J xin hh S + i S xn | S − i ′ ii ω + X n [ J yin hh S xi S + n | S − i ′ ii ω + (1 / yin hh S xi S − n | S − i ′ ii ω ] ,ω hh S − j | S + j ′ ii ω = − h S xj i δ j,j ′ + X m J xjm hh S − j S xm | S + j ′ ii ω − X m [ J yjm hh S xj S − m | S + j ′ ii ω + (1 / yjm hh S xj S + m | S + j ′ ii ω ] . In the RPA [10] for all GFs the following approximationis used for the lattice sites n = i, m = j , as e.g., hh S xi S αn | S βi ′ ii ω = h S xi i hh S αn | S βi ′ ii ω = σ hh S αn | S βi ′ ii ω , hh S xn S αi | S βi ′ ii ω = h S xn i hh S αi | S βi ′ ii ω = − σ hh S αi | S βi ′ ii ω , (5)where h S xi i = σ for i ∈ A while h S xn i = − σ for n ∈ B . A similar approximation is used for the B sublattice,where h S xj i = − σ for j ∈ B while h S xm i = σ for m ∈ A .The RPA results in a closed system of equations for thecomponents of the matrix GF (4). To solve the obtained system of equations weintroduce the Fourier representation of spin op-erators for N/ S ± i = p /N P q S ± q exp( ± i qr i ) and S ± j = p /N P q ′ S ± q ′ exp( ± i q ′ r j ) , where q and q ′ run over( N/
2) wave vectors in the reduced BZ of each sublattice.Using this transformation the equation for the Fourierrepresentation of the matrix GF (4) can be written inthe from ˆ G ( q , ω ) = { ω ˆ I − ˆ V ( q ) } − × σ ˆ I , (6)where ˆ I is the unity matrix, ˆ I is a diagonal matrix withthe elements d = d = 1 and d = d = −
1, andthe interaction matrix is given byˆ V ( q ) = A B ( q ) C ( q )0 − A − C ( q ) − B ( q ) B ( q ) C ( q ) A − C ( q ) − B ( q ) 0 − A . (7)Here the interaction parameters are: A = σ J x (0) = σ [ J (0) + 2 Γ x ] ,J ( q ) = 2 J (cos q x + cos q y ) + 2 J z cos q z ,B ( q ) = σ Γ y cos q y , C ( q ) = σ [ J ( q ) + Γ y cos q y ] . (8)The spectrum of spin waves is determined from the equa-tion Det | ω ˆ I − ˆ V ( q ) | = 0 . (9)After some algebra we obtain the biquadratic equationfor the frequency ω of spin-wave excitations: ω − ω [ A + B ( q ) − C ( q )] + [ B ( q ) − C ( q )] − A [ C ( q ) + B ( q )] + A = 0 . The solution of this equation reads ω ν ( q ) = ±{ A + B ( q ) − C ( q ) + 2 νA B ( q ) } / ≡ ± σ ε ν ( q ) , (10)where ν = ±
1. The energy of excitations for “acoustic” ε − ( q ) and “optic” ε + ( q ) modes are ε − ( q ) = n J (0) − J ( q ) + 4 Γ x [ J (0) + Γ x ] −− y [ J (0) + J ( q ) + 2 Γ x ] cos q y o / , (11) ε + ( q ) = n J (0) − J ( q ) + 4 Γ x [ J (0) + Γ x ] ++ 2 Γ y [ J (0) − J ( q ) + 2 Γ x ] cos q y o / . (12)These two branches are coupled by the relation ε − ( q + Q ) = ε + ( q ) for the AF wave vector Q .For the symmetric compass-model interaction, Γ x =Γ y = Γ, for q = 0 we have the gapless acoustic modewhile the optic mode has a gap: ε − (0) = 0 , ε + (0) = 2 p Γ J (0) + 2Γ > . (13)For the wave vector q = Q we have the opposite re-sults: ε − ( Q ) = ε + (0) > , ε + ( Q ) = ε − (0) = 0 . In theanisotropic case Γ x > Γ y the spectrum of excitations hasgaps both at q = 0 and q = Q : ε − (0) = ε + ( Q ) = 2 q (Γ x − Γ y ) ( J (0) + Γ x ) . (14)For a conventional AF Heisenberg model with Γ x = Γ y =0 we have only one branch with the dispersion ε − ( q ) = ε + ( q ) = p J (0) − J ( q ) which is gapless both at q = 0and q = Q .A similar equation of motion method for the matrixGF (4) can be employed in the LSWT using the trans-formation S + i = √ S a i , S − i = √ S a † i , S xi = S − a † i a i for the sublattice A and the similar transformation forthe sublattice B ( a i → b † i ). Then keeping only linearterms in the bose-like operators ( a i , a † i ) and ( b i , b † i ) weobtain Eqs. (10), (11), (12) for the spin-wave spectrum inLSWT with the sublattice magnetization σ substitutedby spin S . The same spectrum in LSWT was obtained inRefs. [5, 6]. Note that in the RPA the energy of spin ex-citations ω ± ( q ), Eq. (10), is reduced in comparison withthe LSWT since σ < S even at zero temperature dueto zero-point fluctuations in the AF state. The spectrum(10) for the symmetric compass model, Γ x = Γ y , is sim-ilar to the spectrum of the anisotropic AF Heisenbergmodel considered in Ref. [15] and Refs. [11, 12].In Figure 1 the spectrum of spin waves ω ± ( q ) inthe plane in RPA for the parameters J = 65 meV,Γ = 3 . IrO [13] is shown at T = 0.The spectrum ω − ( q ) shows a gap at the wave vector Q given by ω − ( Q ) = 2 σ p Γ J (0) + 2Γ ≈ . J p Γ /J ≈
22 meV for σ = 0 .
37. This value is comparablewith the maximum energy of excitations ω max − ( Q /
2) =4 σ J p /J ≈ . J that gives ω − ( Q ) /ω max − ( Q / ≈ .
22. We can suggest that the spin-wave spectrum inBa IrO should be similar to that one measured by RIXSin Sr IrO [9]. The latter was fitted by a one-branchphenomenological J − J ′ − J ′′ model with J = 60 meV, J ′ = −
20 meV, and J ′′ = 15 meV. The spectrum doesnot reveal a gap in the acoustic branch ω − ( q ) at Q asfor Ba IrO . However, since the intensity of scatter-ing on magnons is proportional to 1 /ω ( q ) , strong scat-tering on the gapless branch ω + ( q ) → q → Q completely suppresses scattering on the gapped ω − ( q )branch. To distinguish scattering on the two branches,high-resolution studies are necessary. A possibility ofobservation of a two-peak structure in the RIXS experi-ments caused by the two-branch spectrum of spin waves isdiscussed in Refs. [12, 16, 17]. For the model (1) with theexchange interaction J ij only between nearest neighborsthe energy of excitations at q = ( π/ , π/
2) , ω − ( q ) = ω + ( q ), is nearly equal to ω ± ( q = π,
0) (up to ± Γ /J ),while in the RIXS experiment ω ( q ) ≈ (1 / ω ( q = π, q (0, 0)( , )( , 0)(0, 0) ( q ) / J FIG. 1: Spectrum of spin-wave excitations ω − ( q ) (bold line)and ω + ( q ) (dashed line) along the symmetry directions in theBZ for the symmetric compass model with Γ x = Γ y = Γ =0 . J and J z = 0. q (0, 0)( , 0)( , )(0, )(0, 0) ( q ) / J FIG. 2: Spectrum of spin-wave excitations ω − ( q ) (boldline) and ω + ( q ) (dashed line) along the symmetry direc-tions in the BZ for the anisotropic compass model withΓ x = 8 . J, Γ y = 4 . J , J z = 0. Circles are numericalresults from Ref. [5]. Figure 2 shows the spin-wave dispersion for largeanisotropic interaction, Γ x = 8 . J, Γ y = 4 . J usedin Ref. [5] in numerical calculations with Lanczos exactdiagonalization. Our RPA calculations give a similar for-mula for the spectrum as in LSWT except for the prefac-tor σ = 0 .
44 instead of S = 1 / / I correspondsto our (1 / J in Eq. (1), and in Fig. (4) of Ref. [5] the en-ergy unit is J c = 10 I . The spectrum reveals a large gapat all wave vectors caused by the large value of Γ x and anoticeable dispersion only along the Γ(0 , → Y (0 , π ) di-rection due to a large, in comparison with J , interactionΓ y = 4 . J .To calculate the sublattice magnetization σ = h S xi i inRPA, we use the kinematic relation S xi = (1 / − S − i S + i for spin S = 1 / σ = 12 − N/ X q h S − q S + q i . (15)The spin correlation function h S − q S + q i can be calculatedfrom the GF hh S + q | S − q ii ω which follows from the GF (6): hh S + q | S − q ii ω = 2 σ a q ( ω )[ ω − ω − ( q )][ ω − ω ( q )] , (16) a q ( ω ) = ω + A ω − [ A + B ( q ) − C ( q )] ω −− A + A [ B ( q ) + C ( q )] . Using the spectral representation for GFs, for the corre-lation function we obtain h S − q S + q i = 2 σ X µ,ν = ± I µν ( q ) N ( µω ν ( q )) , (17)where N ( ω ) = [exp( ω/T ) − − , and the contributionfrom the four poles of the GF (16) is given by I µν ( q ) = a q ( µω ν ( q ))8 µνω ν ( q ) AB ( q ) . (18)Note that I µν ( q ) does not depend on σ .Using relation (17) we perform the self-consistent solu-tion of Eq. (15) for the magnetization σ . Figure 3 shows T/J
FIG. 3: Sublattice magnetization σ = h S xi i for theparameters J z = 5 · − J, Γ x = 0 . J for Γ y / Γ x = 1 (solidline), 0.95 (dashed line), 0.5 (dotted), and for Γ y / Γ x . the sublattice magnetization for J z = 5 · − J, Γ x =0 . J for various Γ y / Γ x . For the symmetric com-pass model, Γ x = Γ y = 0 . J , the N´eel temperature T N = 0 . J = 275 K is close to T N = 240 K observedin experiment for Ba IrO . We stress that the anisotropyof the compass-model interaction, Γ y / Γ x < T N .To study the T N dependence on the parameters of themodel, we consider Eq. (15) in the limit σ →
0. In thislimit N ( ω ν ) ≈ ( T /σε ν ), and for the N´eel temperature wehave the equation:12 = 1 N/ X q X µ,ν = ± I µν ( q ) 2 T N µε ν ( q ) . (19) -5 -4 -3 -2 -1 00.00.20.40.60.81.0 T N / J Log (J z /J) FIG. 4: N´eel temperature T N as a function of J z with Γ x =Γ y = 0 . J (solid line) and Γ x = Γ y = 0 (dashed line). Therefore, T N = 14 C , C = 1 N/ X q X µ,ν I µν ( q ) µε ν ( q ) . (20)Let us study in which cases the integral over q in (20)has a finite value that results in a finite T N .At first we consider the symmetric compass model,Γ x = Γ y = Γ. In this case ε − ( q ) = 0 at q = 0 and ε + ( q ) = 0 at q = Q . Since these two branches aresymmetric, we can consider only the divergency of theintegral in (20) at q = 0 for ε − ( q ) given around q = 0by ε − ( q ) = 2( J (0) + Γ) { J q x + [ J + Γ / ( J (0) + Γ) ] q y + J z q z } . (21)The integral in (20) diverges as R d q /ε − ( q ) if any coef-ficient before q x , q y or q z in (21) is zero. In particular,for nonzero J (0) there is no LRO at finite T for J z = 0.In the limiting case Γ → I µν ( q ) =( A + µω q ) / (4 µω q ) with ω q = p A − C ( q ). FromEq. (20) we get the conventional formula for T N of theAF Heisenberg model (c.f. Ref. [18]) : T N (Γ = 0) = ( J (0) N X q J (0) − J ( q ) ) − . (22)Thus, for a symmetric 2D compass model we have noLRO at finite T . To obtain LRO, we must have finitevalues of both J and J z . The N´eel temperature T N as a function of the interplane coupling J z is shown inFig. 4 for the interaction Γ x = Γ y = 0 . J and forΓ x = Γ y = 0. We can conclude that the compass-modelinteraction enhances the N´eel temperature and, in par-ticular, the anisotropy of the compass-model interactionresults in a further increase of T N as shown in Fig. 3. Inthe anisotropic case Γ x > Γ y the spectrum of excitationshas a gap at q = 0, Eq. (14), and therefore neitherbranch of this spectrum ever reaches zero, so that wehave a finite T N even for J z = 0. Figure 5 demonstrates T N / J x /J J / T N -log ( x /J) FIG. 5: N´eel temperature T N as a function of Γ x for J z = 0,Γ y = 0 . x (solid line) and Γ y = 0 . x (dashed line). In theinset the 1 /T N dependence is shown in the logarithmic scalefor small Γ x . the dependence of T N on Γ x for J z = 0, Γ y = 0 . x and Γ y = 0 . x . For Γ x → q = 0 or at the AF wave vector Q fornonzero compass-model interactions. The calculation ofthe N´eel temperature T N shows that for the symmetriccompass-model interaction, Γ x = Γ y , and a nonzeroexchange interaction J , the AF LRO at finite T canexist only for a finite coupling J z between the planes.For the anisotropic compass-model interaction, Γ x > Γ y ,and a finite exchange interaction J in the plane, the AFLRO with finite N´eel temperature emerges even in the2D case as observed in finite cluster calculations [5, 6].In any case, T N is enhanced by the compass-modelinteraction.The authors would like to thank G. Jackeli, A. M. Ole´s,J. Richter, Yu. G. Rudoy, and V. Yushankhai for valuablediscussions. We thank T. Nagao who drews our attentionto Refs. [11, 12, 16, 17] for comments on the first versionof our paper [19]. A financial support by the Heisenberg–Landau Program of JINR is acknowledged. Appendix A: Further distant neighbors
To fit the spin-excitation spectrum observed by RIXSin Sr IrO in Ref. [9] we consider a more general exchangeinteraction which includes the NNN interaction J nnij = J ′ ij + J ′′ ij , where J ′ ij = J ′ [ δ r j , r i ± ( a x + a y ) + δ r j , r i ± ( a x − a y ) ]and J ′′ ij = J ′′ ( δ r j , r i ± a x + δ r j , r i ± a y ). Note that in thetwo-sublattice model for the Hamiltonian (1) the in-planeexchange interaction J ij acts between the NN on thetwo sublattices, while for the NNN exchange interaction J nnij the lattice sites i and j refer to the same sublattice. Therefore, in the equation for the GF (6) the exchangeinteraction J nnij gives a contribution only for the diagonalterms in the interaction matrix (7). This results in thetransformation of the parameter A = σ J x (0) to the func-tion A ( q ) = σ [ J x (0) − J nn (0)+ J nn ( q )], where J nn ( q ) =4 J ′ cos q x cos q y + 2 J ′′ (cos 2 q x + cos 2 q y ). The solution ofEq. (9) yelds the same spectrum of spin excitations (10)with the parameter A substituted by A ( q ). The energyof excitations for “acoustic” ε − ( q ) and “optic” ε + ( q )modes is given by Eqs. (11) and (12), where instead of J (0) we have to use the function [ J (0) − J nn (0)+ J nn ( q )].If we take Γ x = Γ y = 0, the spectrum of spin waves trans-forms to ε ( q ) = { [ J (0) − J nn (0) + J nn ( q )] − J ( q ) } / which is gapless both at q = 0 and q = ( π, π ). q (0, 0)( , )(0, )(0, 0) ( q ) / J FIG. 6: Spectrum of spin-wave excitations ω − ( q ) (bold line)and ω + ( q ) (dashed line) along the symmetry directions in theBZ for the model with Γ x = Γ y = Γ = 0 . J , J z = 0, andNNN interactions J ′ = − (1 / J , J ′′ = (1 / J . Taking into account the J nn ( q ) term with J ′ = − (1 / J and J ′′ = (1 / J as suggested in experiment[9] and in Refs. [11, 12] we obtain the spin-wave spec-trum in RPA shown in Fig. 6. In comparison with Fig. 1now the excitation energy ω ( π/ , π/ ≈ (1 / ω ( π, ω ( π,
0) = 2 . J is larger than in Fig. 1, where ω ( π, ≈ . J , but it is still smaller than the experimental valueof 200 meV [9] due to the renormalization of the spin-excitation energy in RPA given by the reduced magneti-zation σ = 0 .
36 in comparison with σ = S = 1 / J ′ = − (1 / J and J ′′ = (1 / J in comparison with J ′ ≈ − . J found in La CuO [20]may be explained by a mixing of the j eff = 1 / j eff = 3 / T N = 0 . J = 220 K (for J =65 meV), in comparison with T N = 0 . J = 275 Kfound for J nnij = 0 and is close to T N = 240 K observedin experiment for Ba IrO . [1] K. Kugel and D. Khomskii, Uspekhi Fiz. Nauk , 621(1982) [translated in Sov. Phys. Usp. , 231 (1982)].[2] R. Or´us, A. C. Doherty, and G. Vidal, Phys. Rev. Lett. , 077203 (2009).[3] S. Wenzel and W. Janke, Phys. Rev. B 78 , 064402(2008).[4] S. Wenzel, W. Janke, and A. M. L¨auchli, Phys. Rev.
E78 , 066702(2010).[5] F. Trousselet, A. M. Ole´s, and P. Horsch, Europhys. Lett. , 40005 (2010),[6] F. Trousselet, A. M. Ole´s, and P. Horsch, Phys. Rev. B86 , 134412 (2012).[7] G. Jackeli and G. Khaliullin, Phys. Rev. Lett. ,017205 (2009).[8] S. Boseggia, R. Springell, H. C. Walker, H. M. Rønnow,Ch. R¨uegg, H. Okabe, M. Isobe, R. S. Perry, S. P. Collins,and D. F. McMorrow, Phys. Rev. Lett. , 117207(2013).[9] J. Kim, D. Casa, M. H. Upton, T. Gog, Y.-J. Kim, J. F.Mitchell, M. van Veenendaal, M. Daghofer, J. van denBrink, G. Khaliullin, and B. J. Kim, Phys. Rev. Lett. , 177003 (2012).[10] S. V. Tyablikov,
Methods in the Quantum Theory ofMagnetism , Plenum, New York, 1967) (2-nd Edition: “Nauka”, Moscow, 1975).[11] J. I. Igarashi and T.Nagao, Phys. Rev.
B 88 , 104406(2013).[12] J. I. Igarashi and T.Nagao, Phys. Rev.
B 89 , 064410(2014).[13] V. M. Katukuri, V. Yushankhai, L. Siurakshina, J. vanden Brink, L. Hozoi, and I. Rousochatzakis, Phys. Rev.
X 4 , 021051 (2014).[14] D. N. Zubarev, Usp. Phys. Nauk, , 71 (1960) [trans-lated in Sov. Phys. Usp. , 320 (1960)].[15] V. I. Lymar’ and Yu. G. Rudoy, Theor. Math. Phys. ,86 (1974).[16] J. I. Igarashi and T.Nagao, J. Phys. Soc. Jpn. B 83 ,053709 (2014).[17] J. I. Igarashi and T.Nagao, Phys. Rev.
B 90 , 064402(2014).[18] A. A. Vladimirov, D. Ihle, and N. M. Plakida, Theor.Math. Phys. , 1540 (2013).[19] A. A. Vladimirov, D. Ihle, and N. M. Plakida, JETPLetters , 780 (2014).[20] R. Coldea, S.M. Hayden, G.Aeppli, T.G. Perring, C.DFrost, T.E. Mason, S.-W. Cheong, Z. Fisk, Phys. Rev.Lett.86