Collective excitations in two-dimensional antiferromagnet in strong magnetic field
aa r X i v : . [ c ond - m a t . s t r- e l ] F e b Collective excitations in two-dimensional antiferromagnet in strong magnetic field
A. V. Syromyatnikov ∗ Petersburg Nuclear Physics Institute, Gatchina, St. Petersburg 188300, Russia (Dated: November 29, 2018)We discuss spin- two-dimensional (2D) Heisenberg antiferromagnet (AF) on a square latticeat T = 0 in strong magnetic field H near its saturation value H c . A perturbation approach isproposed to obtain spectrum of magnons with momenta not very close to AF vector in the leadingorder in small parameter ( H c − H ) /H c . We find that magnons are well-defined quasi-particles at H > . H c although the damping is quite large near the zone boundary. A characteristic rotonlikelocal minimum in the spectrum is observed at k = ( π,
0) accompanied by decrease of the dampingnear ( π, PACS numbers: 75.10.Jm, 75.50.Ee, 75.40.Gb
I. INTRODUCTION
Spin- two-dimensional (2D) Heisenberg antiferromagnet (AF) on a square lattice has been one of the most attrac-tive theoretical objects in the last two decades because this model describes parent compounds of high- T c supercon-ducting cuprates . Whereas the theory of long-wavelength magnons in quantum square 2D AF is well developed anddescribes well existing experimental data , there are some surprising recent experimental findings indicating thatthe standard theoretical approaches do not work for short-wavelength magnons. Thus, a rotonlike local minimum wasobserved at small T in the spin-wave dispersion at k = ( π,
0) in a number of recent experiments on square spin- . This minimum is not reproduced quantitatively within second order of 1 /S expansion and phase flux RVBtechniques, while numerical computations using Quantum Monte-Carlo and series expansion describe it satisfactorily(see Ref. and references therein). The origin of the rotonlike minimum has not been clarified yet. It is attributed tothe entanglement of spins on neighboring sites . Rotonlike minimum in the spectrum of short-wavelength magnons(at k = ( π, π/ √ triangular AF at H = 0 obtained using 1 /S expansion and series expansion .Even more peculiar finding concerning short-wavelength magnons in square spin-
2D AF was obtained in Ref. .The authors investigated renormalization of the spin-wave spectrum within first order in 1 /S in strong magneticfield H smaller than the saturation field H c . Due to noncollinearity of sublattices and the field, three-magnon termsappear in the Hamiltonian after transformation of spin operators to bosons. Such terms lead to the loop diagramfor the self-energy part Σ( k ) in the first order in 1 /S corresponding to processes of magnon decay into two magnons.The spin-wave damping obtained within the first order in 1 /S was not very large. To discuss decaying processesself-consistently authors went beyond the first order in 1 /S carrying out a self-consistent calculation of the Green’sfunction (GF) G ( ω, k ) using the dressed GF considering the loop diagram. As a result they obtained that Im G ( ω, k )as a function of ω at fixed k does not resemble a peak at all in most of the Brillouin zone (BZ) when 0 . H c < H < H c (the worse situation was at H ≈ . H c ). Thus the authors concluded that magnons are unstable according to decayinto two magnons at H ∼ H c . Magnons reappeared only at H > . H c due to decreasing of the three-magnonvertex. Meantime within the approach suggested in Ref. only some higher-order 1 /S corrections are taken intoaccount and it is not clear what is the role of the ignored 1 /S terms. This question looks very serious because there isno small parameter in 1 /S expansion for S = 1 /
2. On the other hand the unsuccessful attempts to describe rotonlikeminimum in the spectrum within second order in 1 /S at H = 0 and results of Ref. themselves signify that higherorder 1 /S terms can renormalize the spectrum of short-wavelength magnons noticeably if S ∼ H c − H ) /H c ≪ k = ( π, π ). The magnon spectrum is calculated below in the leading order in ( H c − H ) /H c . Our results can betested experimentally on the growing family of AF compounds with small or moderate exchange coupling constant J (i.e., with experimentally achievable H c ∝ J ) which can be modeled by square spin- Heisenberg AF . Theapproach suggested below can be used in discussion of short-wavelength excitations in other 2D Bose gases of particlesor quasi-particles.The rest of the paper is organized as follows. We discuss the approach in Sec. II. Results of the spectrum calculationare presented in Sec. III. Sec. IV contains our conclusion. One appendix is added devoted to diagrams estimations. (a) –q+k q+k –p+k k k G( k ) G( p,k ) G( q,k ) (b) (c) S k k ª + = + p+k p+k k k –p+k –p+k k k k k k k q (d) p+k FIG. 1: (Color online.) (a) Equation for the vertex Γ( p, k ), where k = ( ω, k ). (b) Expression for the self-energy part Σ( k ) inthe leading order in µ/H c for k not very close to k and ( π, k ) which aremuch smaller than those presented on slide (a) if µ ≪ H c and k is not very close to k and ( π, k ) at k ∼ ( π, II. MODEL AND TECHNIQUE
We discuss spin- Heisenberg AF on the square lattice in magnetic field which Hamiltonian has the form H = 12 X i,j J ij S i S j + H X i S zi . (1)We express spin projections via Pauli operators a † i and a i as S zi = − + a † i a i , S † i = a † i , S − i = a i . To treat them asBose operators one should introduce to Hamiltonian (1) the constraint term U P i a † i a † i a i a i , where U → ∞ , describinginfinite repulsion of particles on the same site . As a result Hamiltonian (1) reads as H = X k ( ǫ k − µ ) a † k a k + 12 N X k + k = k + k ( J k − k + U ) a † k a † k a k a k , (2)where ǫ k = ( J + J k ), J k = P i J ij cos( kR ij ), R ij connects sites i and j , µ = H c − H , H c = J is the saturationfield, and N is the number of spins in the lattice. The bare spectrum ǫ k is quadratic near k .Formally Eq. (2) describes a dilute Bose gas with pare interaction between particles and with the chemical potential µ . In 3D AF this equivalence leads to the solution at µ ≪ H c by borrowing the results from the theory of the diluteBose gas . Within this approach the long-range AF ordering in the plane perpendicular to the field appearing at H < H c corresponds to condensation of magnons with momentum equal to k . The density of particles is proportionalto µ = H c − H . Therefore one can use results for dilute Bose gas when µ ≪ H c . The ladder approximation for thevertex is valid and one has for the vertex the equation shown in Fig. 1(a). Normal Σ( k ) and anomalous Π( k ) self-energy parts are expressed via the vertex, where k = ( ω, k ) and the same notation is used below. In particular onehas at H < H c for Σ( k ) the expression shown in Fig. 1(b) in the leading order in µ/H c and corrections to it are smallfrom other diagrams some of which are shown in Figs. 1(c) and 1(d).It is tempting to apply the same approach to 2D AF at T = 0 but in contrast to 3D dilute Bose gas the theory of2D one at T = 0 can be developed only in the limit of exponentially small density of particles . It means in the caseof 2D AF that there is a good perturbation theory only for fields exponentially close to H c , ln( H c /µ ) ≫
1, while itis desirable to extend the theory to a more acceptable smallness of the parameter µ/H c . We show now that one canovercome this obstacle and construct a perturbation theory on µ/H c to obtain the magnon spectrum for momenta | k − k | ∼ (0,0) ( p , (0, p ) ( p, p ) FIG. 2: (Color online.) A quarter of Brillouin zone is shown and an area around k is drawn in red in which the suggestedperturbation approach does not work when ln( H c /µ ) ∼ We remind first that self-energy parts are small at µ ≪ H c and they vanish at µ = 0 . Then, one notices thatfor such k that the inequality ǫ k − µ ≫ | Σ( ǫ k , k ) | (3)satisfies, one can neglect | Π( k ) | in the GF denominator D ( k ) = ( ω − ǫ k + µ − Σ( k ))( ω + ǫ k − µ + Σ( − k )) + | Π( k ) | and the normal GF G ( k ) = ( ω + ǫ k − µ + Σ( − k )) / D ( k ) (see Refs. ) has the form G ( k ) = 1 ω − ( ǫ k − µ + Σ( k )) . (4)It is implied in Eq. (3) that | Σ( ǫ k , k ) | ∼ | Π( ǫ k , k ) | that is really the case in the entire BZ as particular estimationsshow. Inequality (3) is not satisfied only in an area around k that grows as the field reduces (see Fig. 2). Thus onehas to calculate only normal self-energy part to obtain the magnon spectrum at | k − k | ∼ for dilute 3D Bose gas, to show usingthe smallness of Σ and Π that we have expression for Σ( k ) presented in Fig. 1(b) in the leading order in µ/H c if k is not very close to k and ( π,
0) (see Appendix for some details). Integrals in other diagrams, some of which arepresented in Figs. 1(c) and 1(d), become less singular at small µ at such external momenta as compared with cases of k ∼ k and ( π, µ/H c . Diagrams, the more simple of which are shown in Fig. 1(d), should be also taken into account at k ∼ ( π, ρ andthe spectrum in the entire BZ which can be found only using self-energy parts near k . The trick here bases onthe fact that momenta of integration near k are essential in the second diagram in Fig. 1(b). Thus, one can put themomentum q = ( ω q , q ) in the vertex equal to k = (0 , k ) and factor out the vertex. As a result we haveΣ( k ) = 4 n Γ(0 , k ) , (5) n = ρ + iN X q = k Z dω π e iωt G ( ω, q ) = 1 N X i h a † i a i i , (6)where t → +0 and n is the density of particles which is expressed in spin- AF via the magnetization M = − P i h S zi i as n = − MN . The problem remains to calculate M but we can make use of results of previous numerical calculations.Magnetization in spin- square 2D AF in magnetic field was investigated before in the second order in 1 /S andby numerical diagonalization of finite clusters . It was obtained in Ref. that expression for the magnetizationderived within the first order in 1 /S fits the numerical data very good except for close vicinity of H c . Thus we canapproximate n to high accuracy by its value obtained in the first order in 1 /S for which one has after calculationsfollowing those of Ref. n = ˜ µ − ˜ µN X q J q J s J + J q J + J q − J q ˜ µ (2 − ˜ µ ) , (7)where ˜ µ = µ/H c . In particular, one obtains from Eq. (7) n ≈ .
083 and 0.048 for H = 0 . H c and 0 . H c , respectively.We conclude from Eq. (7) that n ∼ µ/H c if H is not exponentially close to H c .It remains only to obtain the vertex Γ(0 , k ). We imply for simplicity that there is exchange coupling betweenneighboring spins only with the exchange coupling constant J . It is convenient to try solution of the equation shownin Fig. 1(a) in the form Γ( p, k ) = α ( k ) + β ( k ) sin (cid:18) p x + 12 k x (cid:19) + γ ( k ) sin (cid:18) p z + 12 k z (cid:19) . (8)Putting this expression in the equation one leads to a set of three linear algebraic equations on α ( k ), β ( k ) and γ ( k )that can be readily solved. The result is quit cumbersome and we do not present it here. One leads to very compactexpressions for the vertex at p = 0 and k = ( ω = ǫ k + iδ, k ) in the particular cases of | k x | = | k z | and k z = 0 havingthe form Γ(0 , k ) = J J k + k ) / T J k + k ) / − J for | k x | = | k z | , (9)Γ(0 , k ) = J (1 − cos k x )4 T (1 − cos k x ) − k z = 0 , (10)where T = 1 N X q J ǫ q + ǫ q + k + k − ǫ k − iδ ) . (11)We use these results below to calculate Σ( k ) for k directed along vectors k and ( π, µ/H c , we discard µ in the expressions for the vertex.It is difficult to calculate the vertex analytically in the entire BZ and we perform in the next section numericalintegration to find it along some particular directions in BZ and for some particular H values. III. RESULTS OF THE SPECTRUM CALCULATION
Results for the momentum scan along k obtained using Eqs. (9) and (11) are plotted in Fig. 3 for H = 0 . H c and 0 . H c . Solid (dashed) line in Fig. 3 is for the area in which inequality (3) does (does not) hold. Momenta k c at which the solid line turns into the dashed one have been found from the assumption | Σ( ǫ k c , k c ) | = 0 . ǫ k c − µ )that gives k c ≈ . k and 0 . k for H = 0 . H c and 0 . H c , respectively. Insets in Fig. 3 show that the ratio ofthe imaginary and the real part of the spectrum does not exceed 0.17 upon increasing | k | up to | k c | . The spin-wave spectrum is also drawn in Fig. 3 obtained within the linear spin-wave theory (LSWT) and having the form p ( J + J k )( J + (2( H/H c ) − J k ). Notice that at k = the spectrum obtained within our approach is real andcoincides with that derived in LSWT. This finding is in agreement with quite general arguments presented in Ref. according which quantum fluctuations do not change the classical spectrum at k = in quantum 2D AF at H = 0.Quantum fluctuations shift down the spectrum at finite k as the first 1 /S corrections do .Let us discuss the spectral function defined as − π Im G ( ω + iδ, k ). As is explained above, the vertex Γ(0 , k ), where k = ( ω, k ), can be found exactly using representation (8) with the following result for | k x | = | k z | Γ(0 , k ) = − ( ω − J ) J (2 T ω + (1 − T ) J ) + J ( − T ω + ( − T ) J ) J k + k ) / − T J k + k ) / J (2 T ωJ + (1 − T ) J − T J k + k ) / ) , (12)where T = 1 /N P q J / ǫ q + ǫ q + k + k − ω − iδ ) . These expressions can be simplified further at ω = ǫ k and one leadsto Eqs. (9) and (11). We plot in Fig. 4 the spectral function for some momenta k k k obtained using Eq. (12). Onecan see that curves look like peaks while peaks widths are quite large near the zone boundary. It should be notedthat Fig. 4 is in contrast to the corresponding figure of Ref. for H = 0 . H c obtained using 1 /S expansion. FIG. 3: (Color online.) Real and imaginary parts of magnon spectrum (in units of the exchange coupling constant J betweennearest spins) for k = ξ k calculated in the leading order in µ/H c for (a) H = 0 . H c and (b) H = 0 . H c . Insets show ratiosof imaginary and real parts of the spectrum. Dashed lines are for area in which the suggested perturbation technique does notwork (inequality (3) does not hold). Results of the scan along another path in BZ are shown in Fig. 5. The most remarkable feature of this scan is therotonlike minimum at ( π,
0) similar to that obtained recently at H = 0 . It should be noted that imaginary partof T diverges at k = ( π,
0) and the vertex vanishes because one has in Eq. (11) ǫ q + ǫ q + k + k − ǫ k = 2 J (1 + cos q x ).Notice also that for k k ( π,
0) the damping is at least of one order smaller than the real part of the spectrum even for H = 0 . H c . As is noticed above, extra diagrams should be taken into account at k ∼ ( π,
0) the more simple of whichare presented in Fig. 1(d) and which sum cannot be calculated analytically. It is shown in the Appendix that theygive zero contribution to ImΣ( k ) at k = ( π,
0) in the leading order and one concludes that there is a decrease of thedamping near ( π, k ) at k ∼ ( π, µ . Then the real part of the spectrum can deviate near ( π,
0) from curves drawn in Fig. 5. Thesediagrams are of different signs at k = ( π,
0) (e.g., the first two diagrams in Fig. 1(d) are negative and the last one ispositive) and we cannot prove rigorously that there is the rotonlike minimum in the real part of the spectrum if H is not very close to H c . Meantime we demonstrate now using the theory of 2D dilute Bose gas that the rotonlikeminimum does exist if ln( H c /µ ) ≫
1. One estimates in the leading order ReΣ( k ) ∼ − µ/ ln( H c /µ ) + Jκ ln( H c /µ ) (thefirst and the second term here stem from the first diagrams shown in Fig. 1(d) and 1(b), respectively) at k = ( π, κ )and κ ≪ µ/H c while ReΣ( k ∼ k ) ∼ µ ln( H c /µ ) and ǫ k ≈ J − Jκ /
2. It seems to us likely that the minimumremains also at smaller H .It is difficult to go beyond the first order in µ/H c because, in particular, one needs to know self-energy parts near k to calculate the next order diagrams. We can estimate them to be O ( µ p µ/H c ) from the following considerationthat works near ( π,
0) only for ImΣ( k ). Using relation Σ( k ) − Π( k ) = µ one obtains for the spectrum near AF vector p ǫ k ( ǫ k + 2Π( k )) . As a result diagrams, some of which are shown in Figs. 1(c) and 1(d), are of the order of ρ p J Π( k ) or smaller (see Appendix for details). On the other hand one concludes from Eqs. (6) and (7), bearing inmind that the second term in Eq. (6) is positive and of the order of Π( k ), that ρ , Π( k ) ∼ µ if H is not exponentially I m ( G ( i , k )) ( / J ) J H=0.9H C H=0.95H C FIG. 4: (Color online.) Spectral function − π Im( G ( ω + iδ, k )) at fixed momenta k = ξ k (the value of ξ is indicated near eachcouple of curves). Only the curve for H = 0 . H c is shown for ξ = 0 . FIG. 5: (Color online.) Same as in Fig. 3 but along the path in BZ shown in insets. Rotonlike local minima are seen at ( π, π, close to H c .Thus, one can expect from estimations by the order of magnitude that the approach discussed in the present paperworks for H > . H c . It should be noted that numerical evidences have appeared recently that short-wavelengthmagnons are unstable in some regions of BZ at H . . H c . IV. CONCLUSION
In conclusion, we discuss spin-
2D Heisenberg AF on a square lattice at T = 0 in strong magnetic field H near itssaturation value H c . A perturbation approach is proposed to obtain spectrum of magnons with momenta not veryclose to AF vector k = ( π, π ) in the leading order in small parameter µ/H c = ( H c − H ) /H c . It is shown that onlynormal self-energy part Σ( k ) contributes to the spectrum renormalization for which we have two diagrams shownin Fig. 1(b). The sum of these diagrams is proportional to the number of particles which is expressed via uniformmagnetization investigated before numerically (see Eqs. (5) and (6)). These diagrams are of the order of µ . Higher-order diagrams some of which are shown in Fig. 1(c) and 1(d) are of the order of O ( µ p µ/H c ) at | k − k | ∼ k ∼ ( π, k ). We find that magnons are well-defined quasi-particles at H > . H c (see Figs. 3 and 4).A characteristic rotonlike local minimum in the spectrum is observed at k = ( π,
0) accompanied by decrease of thedamping (see Fig. 5).The approach suggested in the present paper can be used in discussion of short-wavelength excitations in 2D Bosegases of particles or quasi-particles. In the latter case one needs to calculate numerically the uniform magnetizationto find the particle density n in Eq. (5). We note also that similar approach basing on the hard-core bosons formalismwas proposed in Ref. for AFs with singlet ground state. Acknowledgments
This work was supported by Russian Science Support Foundation, President of Russian Federation (grant MK-1056.2008.2), RFBR grant 07-02-01318, and Russian Programs ”Quantum Macrophysics”, ”Strongly correlated elec-trons in semiconductors, metals, superconductors and magnetic materials” and ”Neutron Research of Solids”.
APPENDIX A: DIAGRAMS ESTIMATIONS
This Appendix is devoted to estimation of diagrams. We demonstrate first that when k is not very close to k and( π,
0) the leading order corrections to the normal self-energy part Σ( k ) come from two diagrams shown in Fig. 1(b)and other diagrams are small some of which are shown in Figs. 1(c) and 1(d). As it is explained in the main text,diagrams presented in Fig. 1(b) are of the order of Jρ ∼ µ .Let us discuss the first diagram in Fig. 1(d). Representing the normal GF as G ( k ) = ( E ( k ) + ω ) / ( ω − ε k ), where E ( k ) = ǫ k − µ + Σ( − k ) and ε k is the renormalized spectrum, one estimates for the first diagram in Fig. 1(d) afterintegration over the frequency and putting k = ( ε k , k ) − ρ N X q E ( q ) − ε q ε q ( E ( q − k + k ) + ε q − k + k )( ε q + ε q − k + k ) ε q − k + k [( ε q + ε q − k + k ) − ε k ] Γ( k − q, q ) . (A1)One has in Eq. (A1) when q ∼ k E ( q ) − ε q ε q ≈ | Π( k ) | ( E ( q ) + ε q ) ε q ≈ | Π( k ) | (Π( k ) + ǫ k + ε q ) ε q , (A2)where we use the relations Σ( k ) − Π( k ) = µ and ε q = p ǫ q ( ǫ q + 2Π( k )). Then, one has when q ∼ k and κ ≫ Π( k ) /H c : ε q ≈ ǫ q + Π( k ) and ε q + ε q − k + k − ε k ≈ J (cid:0) κ − ( κ x cos k x + κ z cos k z ) − κ x sin k x + κ z sin k z ) (cid:1) + Π( k ) , (A3)where κ = q − k . One concludes from the above consideration bearing in mind that Σ , Π . µ that momenta ofsummation q near k are essential in Eq. (A1).Let us discuss the point k = ( π , π ). One has from Eq. (A3) in this case ε q + ε q − k + k − ε k ≈ − J ( κ x + κ z ) + Π( k ).As a result one obtains from Eqs. (A1), (A2), and (A3) that the first diagram shown in Fig. 1(d) is proportional to iJρ p µ/H c ∼ iJ ( µ/H c ) / , where we use also that Π( k ) ∼ µ , as it is explained in the main text.One leads to the same estimations of Eq. (A1) for other k except for k ∼ ( π,
0) and k ∼ k . One has at k = ( π, ε q + ε q − k + k − ε k ≈ Jκ x +Π( k ). As a result the first diagram shown in Fig. 1(d) at k = ( π,
0) contributesto the main correction to ReΣ( k ) being of the order of Jρ ∼ µ . Its contribution to ImΣ( k ) is of the higher order in theparameter µ/H c . The same is true also for k ∼ k because we have at k = k : ε q + ε q − k + k − ε k = 2 ε q ≈ Jκ +Π( k ).We assume in the above estimations that Γ( k − k , k ) ∼ J in Eq. (A1). It is really the case as the resultsshow of particular calculations of Γ( k − k , k ) carried out as it is described in the main text. Meantime the vertexΓ(0 , k ) ≪ J in the vicinity of the point k = ( π,
0) because T diverges (see the main text). Thus only diagramscontaining Γ( k − k , k ) contribute to the leading correction to ReΣ( k ) at k ∼ ( π, ∗ Electronic address: [email protected] E. Manousakis, Rev. Mod. Phys. , 1 (1991). S. Chakravarty, B. I. Halperin, and D. R. Nelson, Phys. Rev. B , 2344 (1989). N. B. Christensen, H. M. Ronnow, D. F. McMorrow, A. Harrison, T. G. Perring, M. Enderle, R. Coldea, L. P. Regnault,and G. Aeppli, Proc. Natl. Acad. Sci. U.S.A. , 15264 (2007). Y. J. Kim, A. Aharony, R. J. Birgeneau, F. C. Chou, O. Entin-Wohlman, R. W. Erwin, M. Greven, A. B. Harris, M. A.Kastner, I. Y. Korenblit, et al., Phys. Rev. Lett. , 852 (1999). H. M. Rønnow, D. F. McMorrow, R. Coldea, A. Harrison, I. D. Youngson, T. G. Perring, G. Aeppli, O. Sylju˚asen, K. Lefmann,and C. Rischel, Phys. Rev. Lett. , 037202 (2001). M. D. Lumsden, S. E. Nagler, B. C. Sales, D. A. Tennant, D. F. McMorrow, S.-H. Lee, and S. Park, Phys. Rev. B ,214424 (2006). O. A. Starykh, A. V. Chubukov, and A. G. Abanov, Phys. Rev. B , 180403 (2006). A. L. Chernyshev and M. E. Zhitomirsky, Phys. Rev. Lett. , 207202 (2006). W. Zheng, J. O. Fjaerestad, R. R. P. Singh, R. H. McKenzie, and R. Coldea, Phys. Rev. Lett. , 057201 (2006). W. Zheng, J. O. Fjaerestad, R. R. P. Singh, R. H. McKenzie, and R. Coldea, Phys. Rev. B , 224420 (2006). M. E. Zhitomirsky and A. L. Chernyshev, Phys. Rev. Lett. , 4536 (1999). F. M. Woodward, A. S. Albrecht, C. M. Wynn, C. P. Landee, and M. M. Turnbull, Phys. Rev. B , 144412 (2002). T. Lancaster, S. J. Blundell, M. L. Brooks, P. J. Baker, F. L. Pratt, J. L. Manson, M. M. Conner, F. Xiao, C. P. Landee,F. A. Chaves, et al., Phys. Rev. B , 094421 (2007). F. C. Coomer, V. Bondah-Jagalu, K. J. Grant, A. Harrison, G. J. McIntyre, H. M. Ronnow, R. Feyerherm, T. Wand,M. Meissner, D. Visser, et al., Phys. Rev. B , 094424 (2007). E. G. Batyev and L. S. Braginskii, Sov. Phys. JETP , 781 (1984). V. N. Popov,
Functional Integrals and Collective Excitations (Cambridge University Press, Cambridge, 1987). S. T. Belyaev, Sov. Phys. JETP , 299 (1958). M. Schick, Phys. Rev. A , 1067 (1971). M. E. Zhitomirsky and T. Nikuni, Phys. Rev. B , 5013 (1998). M. S. Yang and K. H. M¨utter, Z. Phys. B , 117 (1997). D. I. Golosov and A. V. Chubukov, Sov. Phys. Solid State , 893 (1988). O. F. Sylju˚asen, Phys. Rev. B , 180413(R) (2008). A. Luscher and A. Laeuchli, arXiv:0812.3420. V. N. Kotov, O. Sushkov, Z. Weihong, and J. Oitmaa, Phys. Rev. Lett.80