Dynamics of Thermal Effects in the Spin-Wave Theory of Quantum Antiferromagnets
DDynamics of Thermal Effects in the Spin-Wave Theory of Quantum Antiferromagnets ´Angel Rivas ∗ and Miguel A. Martin-Delgado Departamento de F´ısica Te´orica I, Facultad de Ciencias F´ısicas, Universidad Complutense, 28040 Madrid, Spain
We derive a master equation that allows us to study non-equilibrium dynamics of a quantumantiferromagnet. By resorting to spin-wave theory, we obtain a closed analytic form for the magnondecay rates. These turn out to be closely related to form factors, which are experimentally accessibleby means of neutron and Raman scattering. Furthermore, we compute the time evolution of thestaggered magnetization showing that, for moderate temperatures, the magnetic order is not spoiledeven if the coupling is fully isotropic.
PACS numbers: 42.50.Lc, 03.65.Yz, 75.30.Ds
I. INTRODUCTION
The properties of the quantum Heisenberg model playa fundamental role in the physics of many-body effectsfor models defined by quantum Hamiltonians on a lat-tice, in several spatial dimensions [1, 2]. One of the firstnon-perturbative methods devised to study the quantumHeisenberg model is known as spin-wave theory (SWT).This is a type of mean-field theory method that is es-pecially suited to study the quantum fluctuations of in-teracting spins. The basic assumption is the existenceof a ground state that spontaneously breaks the globalsymmetry of the Heisenberg Hamiltonian. In this case, itcorresponds to rotational symmetry SO(3) about an ar-bitrary axis. In SWT this symmetry is broken by fixinga preferred axis called magnetization axis of the groundstate, and excitations appear in the form of fluctuationsfrom the fixed direction. These are the Goldstone bosonsof this spontaneously breaking mechanism and repre-sent the magnon modes propagating as spin waves inthe quantum system. However, spatial dimensionalityis crucial in order to have a well-defined semiclassical ex-pansion in the parameter 1 /S , where S is the total spinat each site of the system lattice. Namely, in a quantumantiferromagnet the spatial dimension of the lattice hasto be large enough in order to sustain the assumptionof a given order in the ground state. Otherwise, strongquantum fluctuations in one-dimensional lattices breakthe long-range order and makes the SWT invalid. How-ever, many interesting systems are materials in 3D, andSWT provides very good approximations to their observ-able quantities.SWT has been extensively developed in many aspects.It has become by now a standard and reference tool inorder to have a good approximate description of quantumantiferromagnetic systems, whenever the validity of itsapplication is justified.To the best of our knowledge, there is an important as-pect of SWT that remains vaguely explored, namely, themodification of SWT in order to adapt it to describe the ∗ Electronic address: [email protected] ; Fax: +34 913945197 natural interaction of a quantum antiferromagnet with anexternal or surrounding thermal bath that is interactingwith it. A typical example is provided by the phonons ofthe lattice, where the quantum spins are located. This isa basic and fundamental problem since it entails the de-scription of both dynamical effects, i.e. time-dependent,as well as finite-temperature effects outside the state ofthermal equilibrium.Embedding thermal fluctuations in the dynamics of asystem may be approached from several points of view.For instance, in the classical domain, it is common to con-sider the effect of a noisy magnetic thermal field actingon the Heisenberg Hamiltonian [3]. However, that situa-tion is different from what we focus in this work, wherethe noise is described from a microscopic model based onthermal excitation of the surrounding environment. Thebranch of the quantum theory that deals with this kindof problems is the theory of open quantum systems [4–7]that plays a fundamental role in quantum informationtheory [8, 9]. From this point of view, the quantum mag-net is considered as an open system, which exchangesenergy with its environment.The best method to describe an open system stronglydepends on the explicit nature of each situation. Forexample, recently an approach based on the non-equilibrium functional renormalization group has beenproposed for the study of the thermalization of a magnongas in contact with a thermal phonon bath [10]. In thiswork, we have applied the Davies formalism, which is asuitable description of an open system weakly interact-ing with a large environment. One of its main features isthat it allows us to derive an evolution equation for anyspin observable of the quantum antiferromagnet coupledto a generic thermal bath at a certain temperature T .Namely, it provides us with an equation for the evolu-tion of the density matrix ρ ( t ). Furthermore, as a con-sequence of how this fundamental equation is obtained,a series of interesting results for the enlarged SWT havebeen obtained: i/ the quantum antiferromagnet thermal-izes towards the Gibbs state for long enough times; ii/ thedecay rate of this thermalization process can be obtainedin a closed analytical form as a function of the latticemomentum; iii/ the thermal bath cannot be arbitrary inorder to ensure the convergence of any observable to its a r X i v : . [ c ond - m a t . s t r- e l ] J a n thermal value, but it has to belong to the class of super-ohmic baths with specific parameters, depending on thequantum antiferromagnet; iv/ the staggered magnetiza-tion can be computed analytically and we can obtain itsbehaviour with time and temperature, thereby unveilingthe fate of the antiferromagnetic order parameter; and v/the thermal evolution of the magnon form factor can alsobe computed explicitly. These quantities are of physicalimportance and observable in inelastic neutron scattering[11–13] and Raman experiments [14] for instance.Let us emphasize that the framework of our investiga-tions is the out-of-equilibrium dynamics in a spin-wavesystem coupled to a bosonic thermal bath. The method-ology employed is the master equation formalism for openquantum systems. With this combination of dynamicsand methodology, we have found new behavior for thespin-wave decay rates at finite temperature, that havenot been treated previously. Earlier studies of dampingeffects in spin waves at finite temperature, such as [15],rely on the use of the Gibbs state at different tempera-tures. Nevertheless, let us note that, in our study, themagnons are damped while the system is approaching theGibbs state, not once the system is in equilibrium withthe environment at some temperature. It is this type ofnew physics that we can address in a different way thanthe previous investigations.As for the physical nature of the coupling between thesystem spin waves and the bosonic external bath, we mayconsider at least two possible practical realizations.a/ Quantum simulations with optical lattices: the ex-perimental realization of a controlled Mott insulator tosuperfluid transition with cold atoms in an optical trap[16] has opened the field to quantum simulations of newphysics in a range of parameters and types of couplingsthat are not easy to find in nature, but they are feasibleto engineer.b/ Interaction with phonons in a crystal lattice: al-though it is natural to think of lattice phonons in a con-densed matter system as a candidate for the bosonic cou-pling to the spin waves, this possibility comes with sev-eral caveats. First, we have employed spin-wave theoryin the first-order approximation (linear spin-wave the-ory (LSWT)). While there are theoretical studies thatsupport the use of these approximations in 3D quantumspin systems at finite temperature [17], not all materialsof this class exhibit a behavior according to LSWT. Nev-ertheless, it is possible to find several types of compoundswhose magnons behave very well as predicted by LSWT[18–21], and these are candidates for the application ofour results. However, we also point out that in addi-tion to the magnon-boson channel studied in this paper,real materials may also have other decay channels dueto magnon-electron interactions or coupled orbital-latticefluctuations [22] that are outside our current framework.This paper is organized as follows: in Sect.II, we re-view the linear spin-wave theory and establish our nota-tion. In Sect.III, we describe the microscopic couplingof the SWT Hamiltonian with a Hamiltonian bath of FIG. 1: Arrangement of a quantum antiferromagnet in theN´eel state on a square lattice in 2D. The color of the spinsdenotes the two different sublattices, A (blue arrows) and B (orange arrows). The true ground state is close to this stag-gered configuration; however, there is a slight disarrangementin the orientation of the spins due to quantum fluctuations. bosonic operators. In Sect.IV, we derive the completemaster equation for describing the evolution and ther-mal effects of the system due to its interaction with thebath. In Sect.V, we compute relevant observables underthe above conditions, such as the staggered magnetiza-tion, two-spin correlators, and form factors. Sect.VI isdevoted to conclusions. We refer to appendix A for ex-pressions of the time evolution of the first and secondmoments and to appendix B for the detailed calculationinvolving the two-spin spatial correlation functions. II. SPIN WAVE THEORY FOR QUANTUMANTIFERROMAGNETS
First, let us briefly recall the spin-wave theory forquantum antiferromagnets and for establishing our no-tation. The system consists of a lattice with a spin S on every vertex. The Hamiltonian contains only two-body terms between the first neighbors according to theHeisenberg interaction H S = J (cid:88) (cid:104) r , r (cid:48) (cid:105) S r · S r (cid:48) , (1)with J > A and B in such a way that the first neighbors of a A belong to B and vice versa, see figure 1 for an illustration of thetwo-dimensional case), the ground state is close to a stag-gered spin configuration known as N´eel state. However,if the lattice is not bipartite (e.g. triangular lattice), thesystem becomes frustrated, and no simple configurationis found to be a ground state for the diagonal part of theHamiltonian (1). For our purposes, we shall consider a3D square lattice.The diagonalization of the Hamiltonian (1) is not aneasy task, and no exact solutions are known for spatialdimensions d ≥ S ≥ d = 1. Thus,approximation methods become very useful. Probablythe most fundamental of them is based on the Holstein–Primakoff approximation [23, 24] and leads to the so-called spin-wave theory [25], which is also applicable toferromagnets [26]. This method rewrites the spin oper-ators in terms of bosonic annihilation and creation op-erators, a and a † , [ a, a † ] = 1. Concretely, for r in thesublattice A S + r = √ Sf S ( a † r a r ) a r ,S − r = √ Sa † r f S ( a † r a r ) , (2) S z r = S − a † r a r , and for r in the sublattice BS + r = √ Sb † r f S ( b † r b r ) ,S − r = √ Sf S ( b † r b r ) b r , (3) S z r = b † r b r − S, with f S ( x ) = (cid:16) − x S (cid:17) / . (4)By writing the Hamiltonian (1) at the first order in f S ( x ),we obtain the so-called linear spin-wave theory: H LSW = J (cid:2) − N dS + 2 dS (cid:80) r ( a † r a r + b † r b r )+ S (cid:80) (cid:104) r , r (cid:48) (cid:105) ( a r b r (cid:48) + a † r b † r (cid:48) ) (cid:105) . (5)This approximation is valid to describe states, where (cid:104) f S ( a † r a r ) (cid:105) = (cid:104) f S ( b † r b r ) (cid:105) (cid:39)
1, and thus, they also ver-ify (cid:104) a † r a r (cid:105) , (cid:104) b † r b r (cid:105) (cid:28) S. (6)This is the self-consistent condition characteristic of thismean-field theory method.The Hamiltonian H LSW is quadratic in boson opera-tors, so in order to diagonalize it, we take Fourier trans-form: a r = (cid:114) N A (cid:88) k e − i k · r a k , (7) b r = (cid:114) N B (cid:88) k e i k · r b k , (8)with N A = N B = N/ k = 2 π m N A,B = 4 π m N , m ∈ A, B. (9)Then, the first term of H LSW is easy to compute, giventhe orthonormalization rule N (cid:80) r ∈ A,B e i k · r = δ k , . Forthe second one, we parameterize r (cid:48) neighbor to r as r (cid:48) = r + ˆ r µ , where ˆ r µ is the unit vector in the µ direction, which in d = 3 and starting from the first site, can be(1 , , , ,
0) or (0 , , H LSW = J (cid:104) − N dS + 2 S (cid:80) k d ( a † k a k + b † k b k )+ ξ k ( a k b k + a † k b † k ) (cid:105) , (10)where ξ k = (cid:88) µ cos( k · ˆ r µ ) . Next step is to perform a Bogoliubov transformationto new boson operators α k and β k : a k = cosh( θ k ) α k − sinh( θ k ) β † k , (11) b k = − sinh( θ k ) α † k + cosh( θ k ) β k . (12)The function θ k is chosen so that the coefficient of α k β k and α † k β † k is zero: tanh(2 θ k ) = ξ k d . (13)With this choice, the Hamiltonian of the system is diag-onalized: H LSW = E + (cid:88) k ω ( k )( α † k α k + β † k β k ) . (14)Here, the energy dispersion relation is ω ( k ) = 2 JS (cid:113) d − ξ k , (15)and E is a constant E = − JN S (cid:34) dS + 2 N (cid:88) k (cid:18) d − (cid:113) d − ξ k (cid:19)(cid:35) . In summary, we have transformed the intricate Hamil-tonian (1) with interaction terms into another approxi-mate Hamiltonian, which is just a collection of uncou-pled harmonic oscillators, and hence, it is easy to writethe whole spectrum analytically. The excitations of theseharmonic oscillators are called “magnons”, because theyrepresent the minimal collective magnetic excitation ofthe spin lattice.
III. INTERACTION WITH A BOSONICENVIRONMENT
The antiferromagnetic system may be affected by adissipative dynamics due to the interaction with its en-vironment. In principle, the most common source of dis-sipation will be bosonic excitations in the lattice (e.g.phonons). Thus, the interaction Hamiltonian will begiven typically by the so-called spin-boson model [5, 27], V ∝ S · R , where S = ( S x , S y , S z ) is the spin vector and R = ( X, Y, Z ) is the position operators of the bosonicenvironment. Other types of coupling (e.g. [28]) couldeventually be taken into account. In addition, and as afirst proposal, we assume a local environmental model: V = (cid:88) j (cid:88) r g ( ω j ) (cid:104) S x r ( A x r ,j + A x † r ,j )+ S y r ( A y r ,j + A y † r ,j ) + S z r ( A z r ,j + A z † r ,j ) (cid:105) . (16)Here, A and A † stand for annihilation and creation op-erators of the environmental boson modes, and we haveassumed that the coupling function g ( ω j ) is isotropic andthe same for every member of the lattice. On the otherhand, the Hamiltonian of the environment is H E = (cid:88) j (cid:88) r ω j ( A x † r ,j A x r ,j + A y † r ,j A y r ,j + A z † r ,j A z r ,j ) , (17)which is written as H E = (cid:88) j (cid:88) k ω j ( A x † k ,j A x k ,j + A y † k ,j A y k ,j + A z † k ,j A z k ,j ) , (18)after taking Fourier transform.In linear spin-wave theory approximation, the interac-tion term reads V LSW = (cid:80) j g ( ω j ) (cid:110)(cid:80) r ∈ A (cid:104)(cid:113) S ( a r + a † r )( A x r ,j + A x † r ,j ) − i (cid:113) S ( a r − a † r )( A y r ,j + A y † r ,j )+ ( S − a † r a r )( A z r ,j + A z † r ,j ) (cid:105) + (cid:80) r ∈ b (cid:104)(cid:113) S ( b r + b † r )( A x r ,j + A x † r ,j )+i (cid:113) S ( b r − b † r )( A y r ,j + A y † r ,j )+ ( b † r b r − S )( A z r ,j + A z † r ,j ) (cid:105)(cid:111) . Now, the whole Hamiltonian has become much more in-volved than the original spin-wave theory Hamiltonian.However, we can consider a simplified version of the in-teraction Hamiltonian V LSW based on the following twofacts: • The terms a † r a r and b † r b r are negligible in compari-son to the others in the regime where the spin-wavetheory is valid (6). • We ignore the term S ( A r ,j + A † r ,j ) because it is afast oscillator, which we may neglect in the weakcoupling limit, see below.Therefore, after taking Fourier transform, we arrive at V LSW = (cid:114) S (cid:88) j (cid:88) k g ( ω j ) (cid:104) ( a k + b † k ) (19) × ( A x − k ,j + A x † k ,j − i A y − k ,j − i A y † k ,j ) + h . c . (cid:105) . Finally, the Bogoliubov transformation of Eqs. (11) and(12) leads to V LSW = (cid:114) S (cid:88) j (cid:88) k g ( ω j ) (cid:18) d − ξ k d + ξ k (cid:19) / × (cid:104) ( α k + β † k )( A x − k ,j + A x † k ,j − i A y − k ,j − i A y † k ,j )+ h . c . ] . (20) IV. MASTER EQUATION FOR A THERMALENVIRONMENT
The dynamics of the system and the environment isgiven by the von Neumann equation dρdt = − i (cid:126) [ H, ρ ] , (21)where H = H LSW + H E + V LSW . (22)We aim at writing a dynamical equation for the state ofthe system ρ S = Tr E ( ρ ), where the trace is taken over theenvironment degrees of freedom. This task is generallyquite complicated. However, we are particularly inter-ested in describing how the system evolves to the Gibbsstate because of the lack of insulation, and such a case isexpected to happen for a large environment in thermalequilibrium (a “bath”) with a small coupling constant.Under these conditions, an equation, called the masterequation, can be found by resorting to perturbation the-ory [29].The initial state of the environment is then written as ρ E = Z − e − βH E (23)= Z − e − β (cid:80) j (cid:80) k ω j ( A x † k ,j A x k ,j + A y † k ,j A y k ,j + A z † k ,j A z k ,j ) , where Z = Tr (cid:0) e − βH E (cid:1) is the partition function with β =1 /k B T . From now on, we shall use natural units (cid:126) = k B = 1.Due to the Riemann-Lebesgue lemma [4], for smallenough [30] coupling g ( ω j ), we can safely neglect thecounter-rotating terms in (20): V LSW = √ S (cid:80) j (cid:80) k g ( ω j ) (cid:16) d − ξ k d + ξ k (cid:17) / (cid:104) α k ( A x † k ,j − i A y † k ,j )+ β † k ( A x − k ,j − i A y − k ,j ) + h . c . (cid:105) . (24)Now, the problem becomes equivalent to two collectionsof uncoupled harmonic oscillators given by their opera-tors α k and β k , which are coupled to a set of independentenvironments characterized by k . The standard tools toobtain a master equation for a weak interaction with theenvironment can be found in references [4–7]. If we applythose techniques to this system, we arrive at dρdt = L ( ρ ) = − i[ H LSW , ρ ] + (cid:88) k γ k (¯ n k + 1) (cid:18) α k ρα † k − { α † k α k , ρ } + β k ρβ † k − { β † k β k , ρ } (cid:19) + γ k ¯ n k (cid:18) α † k ρα k − { α k α † k , ρ } + β † k ρβ k − { β k β † k , ρ } (cid:19) (25)Here γ k := 2 πS (cid:115) d − ξ k d + ξ k J ( ω ( k )) , (26)where J ( ω ) = (cid:80) j g ( ω ) δ ( ω − ω j ) is the so-called spec-tral density of the bath. This one, for solid-state environ-ments, is usually parameterized in the continuous limit[5, 27] as J ( ω ) = αω s ω s − c e − ω/ω c , (27)where α accounts for the strength of the coupling and ω c is the cut-off frequency of the bath. Typically three casesare distinguished: s > s = 1 (ohmic),and s < n k is themean number of phonons in the bath with frequency ω ( k ): ¯ n k := [exp( ω ( k ) /T ) − − . (28) A. Approach to the Equilibrium
By construction [29], the Gibbs state ρ th = Z − e − H LSW /T , at the same temperature T as the bath,is the steady state of equation (25), i.e. L ( ρ th ) = 0. Thisis straightforwardly verified by taking into account thate − H LSW /T α k = e ω ( k ) /T α k e − H LSW /T , (29)e − H LSW /T β k = e ω ( k ) /T β k e − H LSW /T . (30)Moreover, any initial state of the system becomes closerand closer to this Gibbs state during time evolution.We have thus constructed a dynamical equation to de-scribe the thermal relaxation process of a quantum anti-ferromagnet. Remember that for the spin-wave theory tomake sense, the number of magnons has to be small (6),so for large bath temperatures this treatment is not validin the long-time limit where the system approaches theGibbs state (which contains a large number of magnonsfor large T ). However the predictions of equation (25)should also agree reasonably well with the exact ones atshort times. B. Magnon Decay Rates
A remarkable property of this system is that every ex-ponent is not allowed in the spectral density (27) in order −2 0 2−202010203040 k x k z = 0 k y −2 0 2−202010203040 k x k z = ± πk y FIG. 2: Magnon decay rate (26) in the first Brillouin zone.The surface for k z = 0 and k z = ± π is depicted on the leftand right, respectively. to obtain finite results for many-body observables. Thisis because quantities, such as magnon decay rates (26)and the thermal number of phonons, become infinite forcertain values of k , so for those values, the spectral den-sity has to approach zero fast enough. Particularly, itrequires a super-ohmic spectral density. It is worth re-calling here that this kind of problems may also arisein simpler systems, for instance, in a single spin whensubject to a pure dephasing environment (see [4]). Theconcrete values of the rest of parameters of J ( ω ) are notvery relevant for our purposes as we always assume tobe in a sufficiently weak interaction regime [31]; we shalltake s = 3 , α = J/ , ω c = max k ω ( k ) = 2 JSd. (31)In figure 2, we have represented two sheets of themagnon decay rate (26) in the first Brillouin zone. On theone hand, we note that the magnon decay rate vanisheson the origin and on the eight corners of the Brillouinzone k = ( ± π, ± π, ± π ). Therefore, magnons with thesemomenta are not affected by the presence of the thermalbath. From a quantum information point of view, thesubspace S = span {| n k (cid:105)| k = ( ± π, ± π, ± π ) } (32)is a decoherence free subspace, where we can store infor-mation robustly. Note that this is true independently ofthe temperature T and the number of spins N .On the other hand, the magnon decay rate reaches themaximum value on the points of a sphere of radius r (cid:39) .
947 centered just at these eight minimum points k =( ± π, ± π, ± π ). Between both cases, there is a transitionthat we have tried to illustrate by taking the values k z =0 , ± π in the figure (given the symmetry of the decay rate,we can use k x,y instead of k z , leading to the same figures).From (25), it is possible to compute the evolution ofany combination of α k and β k in the Heisenberg picture.We give the result for the evolution of the first and sec-ond moments in Appendix A. Those expressions allowus to compute the evolution of any spin operator in thequantum antiferromagnet and relevant observables con-structed out of them. V. DYNAMICS OF RELEVANT OBSERVABLES
In this section, we study the time evolution of someproperties that have special interest in the description ofa quantum antiferromagnet. For concreteness, we haveselected two of them: the staggered magnetization andthe spin correlation functions.
A. Staggered magnetization
Due to the isotropy of Hamiltonian (1), one may expectthe ground state also to be symmetric under rotations.However, as we have already mentioned, the ground stateturns out to be close to the N´eel state, which has clearlya privileged orientation. This is an example of sponta-neous symmetry breaking [23] . The figure of merit tocompute this order in a quantum antiferromagnet is theexpectation value of the staggered magnetization opera-tor: ˆ m st z = 1 N (cid:88) r ( − (cid:107) r (cid:107) S z r , (33)which in the thermodynamic limit reads m st = lim N →∞ (cid:104) ˆ m † z (cid:105) = lim N →∞ N (cid:88) r ( − (cid:107) r (cid:107) (cid:104) S z r (cid:105) . (34)By using the equations (2) and (3), we may write thisoperator asˆ m st z = 1 N (cid:88) r ( S − n r )= S − N (cid:88) r ˆ n r = S − N (cid:88) k ˆ n k , (35)with ˆ n k = ˆ n ( a ) k + ˆ n ( b ) k . Note that k = 2 π m / ( N/ m varies two by two instead of one by one. So in thethermodynamic limitlim N →∞ N (cid:88) k = 12(2 π ) (cid:90) B . Z . d k . (36) Here, B.Z. stands for the first Brillouin zone, and theextra factor 1 / k on the left-hand side. Thus, thestaggered magnetization becomes m st = S − π (cid:90) B . Z . d k (cid:104) ˆ n k (cid:105) . (37)When the system is interacting with a thermal bath,the staggered magnetization approaches in time to itsthermal value. This is exactly zero for any T (cid:54) = 0 due tothe Mermin-Wagner theorem [32] in 1D and 2D; however,that is not the case in 3D. Additionally, note that theinteraction Hamiltonian (16) is also isotropic, so it is nottrivial to find also magnetic order when the quantumantiferromagnet is not isolated.From the master equation (25), we are able to visualizehow staggered magnetization varies as a function of time.For this aim, we just need to find the evolution of theobservables ˆ n k . In terms of the operators α k and β k , wehave ˆ n ( a ) k = cosh ( θ k ) α † k α k − sinh( θ k ) cosh( θ k )( α † k β † k + α k β k ) + sinh ( θ k )( β † k β k + 1) , (38)ˆ n ( b ) k = sinh ( θ k )( α † k α k + 1) − sinh( θ k ) cosh( θ k )( α † k β † k + α k β k )+ cosh ( θ k ) β † k β k . (39)Particularly, if we start from the ground state, (cid:104) α † k α k (0) (cid:105) = (cid:104) β † k β k (0) (cid:105) = (cid:104) α k β k (0) (cid:105) = 0, we find (cid:104) ˆ n ( a ) k ( t ) (cid:105) = (cid:104) ˆ n ( b ) k ( t ) (cid:105) = cosh(2 θ k )¯ n k (cid:0) − e − γ k t (cid:1) + sinh ( θ k ) . Introducing these values in (37) and using equation (13), m st = m st0 − d π (cid:90) B . Z . d k (cid:32) ¯ n k (cid:112) d − ξ k (cid:33) (cid:0) − e − γ k t (cid:1) , (40)where m st0 = S − π (cid:90) B . Z . d k (cid:32) d (cid:112) d − ξ k − (cid:33) (41)is the expectation value of the staggered magnetization inthe ground state. In 3D, for a square lattice and S = 1 / m st0 (cid:39) . m st . This is due to its dependence on t through the in-tegral of (40), which renders combinations of differentexponentials. Remarkably, there is a short period, wherethe order is lost very fast (between t = 0 and t ∼ . /J ,see the inset figure). After that, the system continues Time (1 /J ) m s t T = 1 . T = 1 K T = 0 . FIG. 3: Decay of the staggered magnetization from theground state showing the approach to the Gibbs state values m st β . The red line corresponds to T = 0 . m st β (cid:39) . T = 1 K and m st β (cid:39) . T = 1 . m st β (cid:39) . evolving slower to the Gibbs state. This suggests that ifwe want to visualize variations of m st due to the envi-ronment, the best chance is to look for them in systemswith not very small J . Other thermal initial conditions lead to similar evolution in the magnetization. B. Two-point correlation functions
It is also worthwhile to study the second moments ofangular momentum operators. For the sake of illustra-tion, we focus in this section on the transversal two-pointspatial correlation function, which is S ⊥ ( r , r , t ) = Tr[ S x r S x r ρ ( t )] . (42)Without loss of generality, we take r ∈ A . Then, for r ∈ A , S ⊥ ( r , r , t ) = S a r + a † r )( a r + a † r ) ρ ( t )] , (43)and S ⊥ ( r , r , t ) = S a r + a † r )( b r + b † r ) ρ ( t )] , (44)for r ∈ B .Details of the computation are found in Appendix B;finally in the thermodynamic limit, we obtain S ⊥ ( r , r , t ) = 2 S (2 π ) (cid:90) B . Z . d k cos[ k · ( r − r )]Θ k ( r ) (cid:34) n k (1 − e − γ k t ) + 1 (cid:112) d − ξ k (cid:35) , (45)where Θ k ( r ) = (cid:40) d, if r ∈ A, − ξ k , if r ∈ B. (46)We have plotted this correlation for some time instantsin figure 4. In addition, figure 5 shows different caseswhen the Gibbs state has been reached. C. Response function
Other interesting quantities in this system are the re-sponse functions. They are the Fourier transform of two-time correlation functions of spin operators, and for in-stance, they directly appear in cross-sections of inelasticneutron scattering, which are experimentally accessible.For an antiferromagnet with staggered magnetization inthe z direction the inelastic scattering is related to thecorrelation (cid:104) S x − k ( t + τ ) S x k ( t ) (cid:105) , where S x k = 1 √ N (cid:88) r e i k · r S x r . (47) One has to be especially careful when computingmultitime-correlation functions for non-unitary evolu-tions. This is because the evolution of the product oftwo operators, say a and b , is not equal to the product ofthe individual evolutions of a and b when the dynamicsis not unitary, i.e. ( ab )( t ) (cid:54) = a ( t ) b ( t ). However, we cancircumvent this problem by writing the correlation func-tion on the extended space where the evolution is indeedunitary: (cid:104) | S x − k ( t + τ ) S x k ( t ) | (cid:105) = Tr (cid:104) [ S x − k ( t + τ ) S x k ( t ) | (cid:105)(cid:104) | ⊗ ρ E ]= Tr (cid:104) e i H ( t + τ ) S x − k e − i Hτ S x k e − i Ht | (cid:105)(cid:104) | ⊗ ρ E (cid:105) . Here, the trace operation is taken over both the sys-tem and the environment degrees of freedom, and H = H LSW + H E + V LSW is the whole Hamiltonian of the sys-tem and the environment. Then, it is possible to obtainthat (see the detailed discussion in [7]) (cid:104) | S x − k ( t + τ ) S x k ( t ) | (cid:105) = (cid:104) | (cid:2) S x − k ( τ ) S x k (cid:3) ( t ) | (cid:105) . (48)That is, it is needed to obtain first the Heisenberg evolu-tion with respect to the parameter τ of the operator S x − k and after that the Heisenberg evolution with respect tothe parameter t of the product S x − k ( τ ) S x k . −10 −8 −6 −4 −2 0 2 4 6 8 10−2−1012345 r − r S ⊥ ( r , r ) −10 0 1000.20.40.6 t = 0 t = 10 − t = 10 − t = 10 − FIG. 4: Evolution of S ⊥ ( r , r , t ) from the ground state for abath with T = 5 K in the thermodynamic limit for differenttime instants (in units of J − ). Note the oscillating behaviortypical of antiferromagnetic systems. The inset figure illus-trates the similarity between the cases for | S ⊥ ( r , r , t ) | afternormalization. For linear spin-wave theory, we have S x r = S + r + S − r (cid:114) S (cid:40) a r + a † r , if r ∈ A,b r + b † r , if r ∈ B, (49)thus, according to (7) and (47), S x k = √ S a k + a †− k + b − k + b † k ) . (50)If we perform the Bogoliubov transformation (11) and(12), the Eqs. (A1) and (A2) lead to S x − k ( τ ) = √ S − γ k τ/ (cid:110) e − i ω ( k ) τ [cosh( θ k ) α − k − sinh( θ k ) β k + cosh( θ k ) β k − sinh( θ k ) α − k ]+ e i ω ( k ) τ (cid:104) cosh( θ k ) α † k + cosh( θ k ) β †− k − sinh( θ k ) α † k − sinh( θ k ) β †− k (cid:105)(cid:111) , (51)where we have used the fact that θ − k = θ k and ω ( − k ) = ω ( k ). Since by assumption γ k is small, for small τ , we canneglect it in comparison to the complex exponential e ± i ω ( k ) τ : S x − k ( τ ) (cid:39) √ S θ k ) − sinh( θ k )] (cid:104) e − i ω ( k ) τ ( α − k + β k ) + e i ω ( k ) τ ( α † k + β †− k ) (cid:105) . (52)Finally, by using (A3) and (A4), we compute the evolution of S x − k ( τ ) S x k with respect to t , and after simplifyingvanishing terms, the correlation function reads (cid:104) S x − k ( t + τ ) S x k ( t ) (cid:105) = S θ k ) − sinh(2 θ k )] (cid:110) e − i ω ( k ) τ [ (cid:104) α − k α †− k ( t ) (cid:105) + (cid:104) β k β † k ( t ) (cid:105) ]+ e i ω ( k ) τ [ (cid:104) α † k α k ( t ) (cid:105) + (cid:104) β †− k β − k ( t ) (cid:105) ] (cid:111) = S ( d − ξ k )2 (cid:112) d − ξ k (cid:110) e − i ω ( k ) τ (cid:2) ¯ n k (cid:0) − e − γ k t (cid:1) + 1 (cid:3) + e i ω ( k ) τ ¯ n k (cid:0) − e − γ k t (cid:1)(cid:111) . (53)The Fourier transform with respect to τ leads to the re-sponse function S ⊥ ( k , t, ω ) = S − LSW ( k , t ) δ [ ω − ω ( k )]+ S + LSW ( k , t ) δ [ ω + ω ( k )] , (54) with S − LSW ( k , t ) = S ( d − ξ k )2 (cid:112) d − ξ k (cid:2) ¯ n k (cid:0) − e − γ k t (cid:1) + 1 (cid:3) , (55) S + LSW ( k , t ) = S ( d − ξ k )2 (cid:112) d − ξ k (cid:2) ¯ n k (cid:0) − e − γ k t (cid:1)(cid:3) . (56)Therefore, at t = 0 (or T = 0) ,only the form factor S − LSQ ( k , t ) remains. On the other hand, we conclude −10 −8 −6 −4 −2 0 2 4 6 8 10−2−10123456 r − r S ⊥ ( r , r ) T = 0 K T = 1 K T = 5 K −10 0 1000.20.40.6 FIG. 5: Thermal values of S ⊥ ( r , r ). The inset representsagain its normalized absolute value. that as temperature increases, S ± LSQ ( k , t ) also increases,and they have the same geometry in momentum space asthe magnon decay rates γ k . VI. CONCLUSIONS
In this work, we have analyzed the behavior of a quan-tum antiferromagnet in contact with a boson thermalbath. Based on the spin-wave theory, we have appliedthe weak coupling procedure (Davies theory) to obtaina master equation for the dynamics. We believe thatthis is a basic and fundamental problem, which has re-mained quite unexplored so far. It is at the crossroadsof strongly correlated systems and the physics of open quantum systems that is so much rooted in quantum in-formation theory.From the open systems point of view, spin-wave the-ory provides us with a nice framework to apply the well-known techniques developed for quantum optics or quan-tum chemistry settings to quantum many-body prob-lems. Interestingly, some features, which are typicallyencountered in small systems under weak coupling limit,e.g. the exponential decay of observables, may be lostwhen computing the observables, which are relevant forthe many-body systems. We have exemplified this pointby studying the staggered magnetization, which for mod-erate temperatures and despite of the isotropic couplingto the bath does not vanish. In fact, it does not show anexponential decay either.Furthermore, we have illustrated the versatility of ourmaster equation approach to the dynamics of thermaleffects in quantum antiferromagnets by computing two-point correlation and response functions, also known asform factors. The geometry in momentum space of theseresponse functions S LSQ ( k , t ) is closely related to that ofthe decay rate function in the first Brillouin zone. Theseform factors, in turn, are directly related to differentialcross-sections in experiments of inelastic neutron scat-tering, which, we believe, that may shed light to thecurrent knowledge of a quantum antiferromagnet undernon-isolated situations. Acknowledgments
We thank the Spanish MICINN grant FIS2009-10061,CAM research consortium QUITEMAD S2009-ESP-1594, European Commission PICC: FP7 2007-2013,Grant No. 249958, UCM-BS grant GICC-910758.
Appendix A: Time-evolution of the first and second moments
The generator L (cid:93) in the Heisenberg picture is obtained by the equality Tr[ X L ( ρ )] = Tr[ ρ L (cid:93) ( X )], for any operator X . By solving the dynamical equations, we obtain α k ( t ) = e [ − i ω ( k ) − γ k / t α k (0) , (A1) β k ( t ) = e [ − i ω ( k ) − γ k / t β k (0) , (A2) α † k α k ( t ) = e − γ k t α † k α k (0) + ¯ n k [1 − e − γ k t ] , (A3) β † k β k ( t ) = e − γ k t β † k β k (0) + ¯ n k [1 − e − γ k t ] , , (A4)The terms α k α † k and β k β † k are obtained by using the commutation relations [ α k , α † k ] = [ β k , β † k ] = 1, and the remainingones are just the composition of the dynamics given in (A1) and (A2) and their Hermitian conjugate.0 Appendix B: Computation of the 2-spin correlator S ⊥ ( r , r , t ) We can compute the evolution of S ⊥ ( r , r ) in the Heisenberg picture. The terms a r a r , a r b † r and their Hermitianconjugate do not contribute to the evolution. For a r a † r , we have( a r a † r )( t ) = 1 N A (cid:88) k , k (cid:48) e − i k · r e i k (cid:48) · r ( a k a † k (cid:48) )( t ) (B1)= 1 N A (cid:88) k , k (cid:48) e − i k · r e i k (cid:48) · r { [cosh( θ k ) α k − sinh( θ k ) β † k ][cosh( θ k (cid:48) ) α † k (cid:48) − sinh( θ k (cid:48) ) β k (cid:48) ] } ( t ) . By using Eqs. (A1)–(A4) and simplifying the mean values which vanish on the ground state, the sum is left only withthe terms, where k = k (cid:48) : (cid:104) a r a † r ( t ) (cid:105) = 1 N A (cid:88) k e − i k · ( r − r ) [cosh ( θ k ) (cid:104) α † k α k ( t ) (cid:105) + sinh ( θ k ) (cid:104) β † k β k ( t ) (cid:105) + cosh ( θ k )]= 1 N A (cid:88) k e − i k · ( r − r ) [cosh(2 θ k )¯ n k (cid:0) − e − γ k t (cid:1) + cosh ( θ k )] . (B2)Similarly, for the remaining non-vanishing terms, we obtain (cid:104) a † r a r ( t ) (cid:105) = 1 N A (cid:88) k e − i k · ( r − r ) [cosh ( θ k ) (cid:104) α † k α k ( t ) (cid:105) + sinh ( θ k ) (cid:104) β † k β k ( t ) (cid:105) + sinh ( θ k )]= 1 N A (cid:88) k e − i k · ( r − r ) [cosh(2 θ k )¯ n k (cid:0) − e − γ k t (cid:1) + sinh ( θ k )] , (B3) (cid:104) a r b r ( t ) (cid:105) = (cid:104) a † r b † r ( t ) (cid:105) ∗ = − N A (cid:88) k e − i k · ( r − r ) cosh( θ k ) sinh( θ k )[ (cid:104) α † k α k ( t ) (cid:105) + (cid:104) β † k β k ( t ) (cid:105) + 1]= − N A (cid:88) k e − i k · ( r − r ) sinh(2 θ k )2 [2¯ n k (cid:0) − e − γ k t (cid:1) + 1] , (B4)Thus, on the one hand, the correlation function for r ∈ A is S ⊥ ( r , r , t ) = 2 SN A (cid:88) k cos[ k · ( r − r )] cosh(2 θ k )[2¯ n k (cid:0) − e − γ k t (cid:1) + 1]+ i 2 SN A (cid:88) k sin[ k · ( r − r )] . (B5)Since 1 N A (cid:88) k sin[ k · ( r − r )] = Im (cid:40) N A (cid:88) k e i k · ( r − r ) (cid:41) = Im ( δ r , r ) = 0 , (B6)and using (13), S ⊥ ( r , r , t ) = 2 SdN A (cid:88) k cos[ k · ( r − r )] (cid:34) n k (1 − e − γ k t ) + 1 (cid:112) d − ξ k (cid:35) . (B7)And on the other hand, for r ∈ B , S ⊥ ( r , r , t ) = − SN A (cid:88) k cos[ k · ( r − r )] sinh(2 θ k )[2¯ n k (cid:0) − e − γ k t (cid:1) + 1]= − SN A (cid:88) k cos[ k · ( r − r )] ξ k [2¯ n k (1 − e − γ k t ) + 1] (cid:112) d − ξ k . (B8) [1] A. Auerbach, Interacting Electrons and Quantum Mag-netism (Springer-Verlag, New York, 1994). [2] J. Gonzalez, M. A. Martin-Delgado, G. Sierra and A. H. Vozmediano,
Quantum Electron Liquids and High-Tc Su-perconductivity (Lecture Notes in Physics Monographs;Springer-Verlag, New York, 1995).[3] W. F. Brown Jr., Phys. Rev. , 16771686 (1963).[4] A. Rivas and S.F. Huelga,
Open Quantum Systems. AnIntroduction (Springer, Heidelberg, 2011).[5] U. Weiss,
Quantum Dissipative Systems (World Scien-tific, Singapore, 2008).[6] H. -P. Breuer and F. Petruccione,
The Theory of OpenQuantum Systems (Oxford University Press, Oxford,2002).[7] C. W. Gardiner and P. Zoller,
Quantum Noise (Springer,Berlin, 2004).[8] M. A. Nielsen and I. L. Chuang,
Quantum Computationand Quantum Information (Cambridge University Press,Cambridge, 2000).[9] A. Galindo and M. A. Martin-Delgado, Rev. Mod. Phys. , 347-423, (2002); arXiv:quant-ph/0112105.[10] J. Hick, T. Kloss and P. Kopietz, arXiv:1206.6689.[11] G. Shirane, Y. Endoh, R. J. Birgeneau, M. A. Kastner,Y. Hidaka, M. Oda, M. Suzuki and T. Murakami, Phys.Rev. Lett. , 1613–1616 (1987).[12] G. Aeppli, S. M. Hayden, H. A. Mook, Z. Fisk, S.-W.Cheong, D. Rytz, J. P. Remeika, G. P. Espinosa and A.S. Cooper, Phys. Rev. Lett. , 2052–2055 (1989).[13] T. G. Perring, G. Aeppli, S. M. Hayden, S. A. Carter,J. P. Remeika and S-W Cheong, Phys. Rev. Lett. ,711–714 (1996).[14] K. B. Lyons, P. A. Fleury, J. P. Remeika, A. S. Cooperand T. J. Negran, Phys. Rev. B , 2353–2356 (1988).[15] G. Cottamm and B. Stinchcombre, J . Phys. C: SolidState Phys. , 2283, 2305, 2326 (1970).[16] M. Greiner, O. Mandel, T. Esslinger, T. W. Hnsch andI. Bloch, Nature , 39–44 (2002).[17] T. Oguchi, Phys. Rev. , 117 (1960). [18] W. C. Koehler, H. R. Child, R. M. Nicklow, H. G. Smith,R. M. Moon and J. W. Cable, Phys. Rev. Lett. , 16–18(1970).[19] K. Clausen, J. J. Rhyne, B. Lebech and N. C. Koon, J.Phys. C: Solid State Phys. , 3587–3596 (1982).[20] A. K. Bera, S. M. Yusuf, N. S. Kini, I. Mirebeau and S.Petit, AIP Conf. Proc. , 1135–1136 (2011).[21] J. Li, V. O. Garlea, J. L. Zarestky and D. Vaknin, Phys.Rev. B , 024410 (2006).[22] T. G. Perring, D. T. Adroja, G. Chaboussant, G. Aeppli,T. Kimura and Y. Tokura, Phys. Rev. Lett. , 217201(2001).[23] E. Manousakis, Rev. Mod. Phys. , 1 (1991).[24] T. Holstein and H. Primakoff, Phys. Rev. , 1098–1113(1940).[25] P. W. Anderson, Phys. Rev. , 694–701 (1952).[26] K. Kubo, Phys. Rev. Lett. , 110 (1988).[27] A. J. Leggett, S. Chakravarty, A. T. Dorsey, M. P. A.Fisher, A. Garg and W. Zwerger, Rev. Mod. Phys. , 1(1987).[28] P. Pincus and J. Winter, Phys. Rev. Lett. , 269–270(1961).[29] E. B. Davies, Commun. Math. Phys. , 91 (1974).[30] By pressuposition we assume that the coupling is alwayssmall enough to safely neglect the counter-rotating termsfor every k , ω ( k ) t (cid:29) g [ ω ( k )]. Note that there is not aproblem with cases where ω ( k ) = 0, because the cou-pling, i.e. the spectral density, also approaches to zero,see Eqs. (26) and (27).[31] A. Rivas, A. D. K. Plato, S. F. Huelga and M. B. Plenio,New J. Phys. , 113032 (2010).[32] N. D. Mermin, H. Wagner, Phys. Rev. Lett.17