Magnon dressing by orbital excitations in ferromagnetic planes of K$_2$CuF$_4$ and LaMnO$_3$

We show that even when spins and orbitals disentangle in the ground state, spin excitations are renormalized by the local tuning of $e_g$ orbitals in ferromagnetic planes of K$_2$CuF$_4$ and LaMnO$_3$. As a result, dressed spin excitations (magnons) obtained within the electronic model propagate as quasiparticles and their energy renormalization depends on momentum ${\vec k}$. Therefore magnons in spin-orbital systems go beyond the paradigm of the effective Heisenberg model with nearest neighbor spin exchange derived from the ground state --- spin-orbital entanglement in excited states predicts large magnon softening at the Brillouin zone boundary, and in case of LaMnO$_3$ the magnon energy at the $M=(\pi,\pi)$ point may be reduced by $\sim 45$\%. In contrast, simultaneously the stiffness constant near the Goldstone mode is almost unaffected. We elucidate physics behind magnon renormalization in spin-orbital systems and explain why long wavelength magnons are unrenormalized while simultaneously energies of short wavelength magnons are reduced by orbital fluctuations. In fact, the ${\vec k}$-dependence of the magnon energy is modified mainly by dispersion which originates from spin exchange between second neighbors along the cubic axes $a$ and $b$.

Orbital degeneracy opens the route towards complex types of spin-orbital order with coexisting antiferromagnetic (AF) and ferromagnetic (FM) exchange bonds. Frequently such systems are analyzed using the classical Goodenough-Kanamori rules [26] which emphasize the complementarity of spin and orbital order, i.e. alternating orbital (AO) order supports FM spin exchange and ferro-orbital order supports AF exchange. They follow from the assumption that spin and orbital excitations are independent of each other and spin exchange interactions may be derived from the spin-orbital superexchange by averaging over the orbital state. Indeed, when joint spin-orbital fluctuations are quenched, e.g. by lattice distortions, these rules apply and the disentangled superexchange helps to understand experimental observations [13]. A good example of this approach is the parent compound of colossal magnetoresistance manganites LaMnO 3 [27], with small spinorbital entanglement [28]. Therefore, spin waves measured in inelastic neutron scattering [29][30][31] and optical spectral weights [32] could be successfully interpreted using the effective anisotropic Heisenberg model.
In doped manganites double exchange provides a large FM exchange interaction [27]. It is responsible for the onset of FM order and modifies occupied e g orbitals involved in the hopping process as demonstrated in the onedimensional (1D) spin-orbital model [33]. Hole-orbital and orbital-lattice fluctuations were identified as the main origin of the observed unusual softening of the magnon spectrum at the zone boundary [34][35][36]. It has been shown that orbitons depend thereby on magnons in Mott insulators with orbital degrees of freedom [37][38][39][40], and both contribute to spectral properties [41,42].
In 1D cuprates [43] (2D iridates [44]) orbitons (excitons) are dressed by magnons, while the opposite effect of orbital excitations on magnons was considered only in the context of the strong zone boundary magnon softening observed experimentally in manganites close to half doping [35]. In this paper we demonstrate that spin excitations in a FM plane with AO order, as in K 2 CuF 4 [45] or LaMnO 3 [46], are indeed renormalized by the changes of occupied e g orbitals, leading to magnons dressed by orbital fluctuations and propagating together as a quasiparticle in a Mott insulator. This phenomenon is similar to the local changes of AF spin order by an added hole in superconducting cuprates [47].
The remaining of the paper is organized as follows. In section 2 we introduce a general form of spin-orbital Hamiltonian with e g degrees of freedom and present the magnon excitations starting from orbital order in the ground state. Next we release the constraint of frozen orbitals and present the variational way of finding magnon excitations for optimized orbitals in section 3. A simplified version of this approach and a numerical ansatz (NA) which serves to verify the predictions of the variational approach are presented in section 4. The results for magnons in K 2 CuF 4 are given in section 5. We consider a spin excitation in the FM planes of LaMnO 3 and analyze the optimal orbital angles near the excitation in section 6. There we also show that the effective spin model will include nearest neighbor J 1 , next-nearest neighbor J 2 , and third next neighbor J 3 spin exchange, although the spin-orbital superexchange couples only nearest neighbors. Analytic estimation of the renormalized interaction J 1 which determines the magnon bandwidth is presented in the appendix. The paper is concluded with a short summary in section 7.

Spin-orbital model and magnons for frozen orbitals
We begin with the e g orbital basis (labeled in analogy to ñ | and ñ | spin = S 1 2 states): When tetragonal distortion is ignored (E z =0), the occupied orbital states on two sublattices A and B are symmetric/antisymmtric combinations of e g orbital basis z The orbital operators t ( ) In K 2 CuF 4 one finds FM order at J H /U∼0.2 [49], coexisting with AO order [50] of hole orbital states with angles ±θ opt on the two sublattices, A and B. Averaging the orbital operators over this AO order in ab planes gives j , which served to interpret the experimental data [51,52]. To investigate magnons (spin waves) we create a spin excitation at site i=0 by decreasing the value of In the simplest approach we disentangle [24] spin-orbital superexchange both in the ground and in excited states and use the same frozen AO order shown in figure 1(a) to determine spin exchange J ◊ .
A spin excitation (a magnon) itself is best described by the transformation to Holstein-Primakoff (HP) bosons [53]. In the linear spin-wave theory magnon energy consists of two contributions and we introduce: The latter originates from quantum fluctuations µ- , . Here d  stands for one of four nearest neighbors of the central site i=0 shown in figure 1(a). The above two terms determine the magnon dispersion in a 2D ferromagnet, which serves as a reference below. The breaking of SU(2) symmetry is reflected by a Goldstone mode (at =  k 0), and w = à  J Sk k 2 for   k 0-we find that this result is insensitive to spin-orbital coupling. In general however orbitals are not frozen in a spin-orbital system and will respond locally to a spin excitation. One might expect that this reduces spin exchange,  àJ J and the magnon dispersion would soften. Indeed, we have found that the magnon energy w  ( ) k 0 is reduced but this effect is rather subtle and the renormalization of exchange interaction J ♦ depends on momentum  k .
Mwhere orbital states in the shaded cluster are radically different from those in the ground state, see frozen orbitals in (a).

Variational approximation
To capture the response of orbital background to a spin excitation we invoke the following variational approximation (VA): significant changes of occupied orbitals with respect to the reference AO order are expected at the nearest neighbors of excited spin and at the site of spin excitation itself, see figure 1 , occurs at the site of spin excitation itself. For the neighboring sites we use the lattice symmetry and search for the same optimal orbitals given by angles 3 at equivalent neighbors along each cubic axis, a or b. It is crucial that the VA is performed for each value of momentum  k independently. We have evaluated the matrix elements of the Hamiltonian  (3) for a single spin excitation in the thermodynamic limit, and determined six variational parameters q { } iL (i=1, 2, 3; L=A, B). In this way we obtained the renormalized magnon dispersion which replaces equation (8), ; . 9 By construction the angles {θ iL } are real as in equation (4); we have verified that complex coefficients do not lead to further significant energy lowering.
The AO order has a unit cell consisting of two atoms which defines the reduced Brillouin zone (RBZ). The magnon dispersion consists then of two branches in the RBZ, the lower one for , is the nesting vector. The two magnon branches give a gapless dispersion and are determined in two steps to take the full advantage of variational parameters. First we find the magnon energies from the lower magnon band-they depend on stand for the Ising HP boson parts, while is obtained from the HP boson propagation along the bonds parallel to the a and b axis, from sublattice A to B (or vice versa). The eigenstates of the second magnon band are determined in the second step-magnon states which momenta do not belong to the RBZ. As a magnon states with momentum Ï  k RBZ is orthogonal to its partner magnon state with momentum -  ( ) k Q , then, at this stage, the variational principle has to be applied together with rigorous orthogonality condition.

Simplified variational approximation (SVA) and NA
Assuming that orbital optimization for both sublattices is equivalent, we use the constraint q q q º = i . The SVA is equivalent to the VA when the magnon happens to be a symmetric linear combination of the two waves propagating over the two sublattices-this concerns the G -M direction; otherwise one may expect that the amplitude of the spin wave is larger in one sublattice and the magnon wave function differs qualitatively from that obtained for a Heisenberg ferromagnet. Below we show that the VA gives indeed better results than the SVA, and the magnon dressing occurs differently on both sublattices.
Finally, we verified the predictions of the VA by exact diagonalization employing a NA with six states per sublattice: a spin defect with or without orbital excitation and four spin-orbital states with spin excitation at the central site together with an orbital excitation at one of nearest neighbors. The state with excitations within a shaded cluster depicted in figure 1(b) may be thus expressed in terms of these six states. Here the constraint for equal orbital angles at two neighbors along the same axis is released. The eigenstates and the spin excitation energy w  k are found separately by exact diagonalization of a 12×12 matrix obtained for each momentum  k .

Magnons for K 2 CuF 4
Taking as an example the K 2 CuF 4 state at = -E J 0.8 z shown in figure 1(b), one finds that the orbital renormalization is large-at the central site with spin excitation it is largely modified to~-( ) x y 2 2 and the orbitals at the four neighboring sites are also changed. The latter orbitals found within the VA are only weakly changed as these latter sites have three neighbors belonging to the neighbors with undisturbed AO order [50], but the one at the spin excitation itself is radically different. For this reason we introduce a cutoff and assume that the orbitals at further neighbors of the excited spin are unchanged. One expects then large dressing of the magnon, with the corresponding reduction of the effective FM interaction to J ♦ , particularly in the neighborhood of the M point. This is confirmed by the results shown in figure 2(a)-the magnon energy ω M is reduced by ∼27% from w ( ) M 0 . Internal consistency of the theory is confirmed by this reduction being nearly the same in all three methods treating spin-orbital coupling: VA, SVA, and NA.
At the X point we recognize the importance of independent optimization of orbitals at the two sublatticesthe energy ω X is reduced by ∼25% from w ( ) path. While the VA may underestimate somewhat the magnon dressing effect, altogether we find indeed a comparison of the VA with the NA very encouraging. The renormalization of magnon energy increases fast when the orbital splitting | | E z is reduced, and one finds that the magnon energy reduction is large for = -E J 0.3 z , e.g. by ∼60% at the M point, see figure 2(b). The agreement between the VA and the NA is somewhat worse here but remains still in qualitative agreement. Altogether, we suggest that the magnon softening may be very large for spin-orbital systems with low spin = S 1 2 as in K 2 CuF 4 .

Magnons for FM planes of LaMnO 3
For LaMnO 3 we consider electrons in e g orbitals at E z >0 and use a representative value of the orbital spitting which gives θ opt ;120°. Spin and orbital excitations depend on θ opt for frozen orbitals [55]. The magnon dispersion is modified within the VA or the SVA, see figure 3(a). In agreement with our initial intuition, the magnon energies predicted for dynamical orbitals soften. The energy lowering from w  ( ) k 0 to w  k is substantial for this value of E z -up to about 45% at the M point. We emphasize that the VA and the NA agree almost perfectly and this agreement confirms a posteriori our initial choice of real orbital phases in the VA. One observes that the energy w  k is somewhat lower than in the NA in the neighborhood of the M point, indicating that the orbitally doubly excited states become important when at least two orbital deformations are large enough (such states are not included in the NA).
We remark that the reported experimental spin exchange constants are the final product of processing the experimental data concerning the magnons energies. A link between them and the measured energies is established by a parameterized form of the dispersion relation for some conceived pure-spin models defined by a specific interactions pattern. In case of LaMnO 3 , the simplest Heisenberg model with nearest neighbor interaction J 1 was successfully used to interpret the experimental data in the past [29]-it predicts the magnon dispersion given by equation (8).
We decided to follow the same strategy and studied our calculated magnon dispersion w  k in figure 3(a). We tested whether or not one may fit the calculated magnon energies with the effective Heisenberg model and how many exchange interactions are needed using the dispersion w  k discretized over a mesh of (k a , k b ) values. It turned out that to reproduce the magnon bandwidth J S 8 1 the fit requires nearest neighbor exchange interaction =´-J J 6.34 10 1 3 , see figure 3(b). Although the reduced value of the magnon bandwidth is then reproduced, the  k -dependence of w  k near the M point is not. It is clear that the Heisenberg model with nearest neighbor exchange is insufficient as the obtained dispersion w  k deviates then from the one derived from the VA, particularly near the M point. The fit may be refined by taking into account the next-nearest and third neighbor interactions, J 2 and J 3 , in the effective spin model. One finds that =´-J J 0.25 10 2 3 is rather small but´- J J 1.35 10 3 3 is significant and plays an important role. Both fits are shown in figure 3(b).
We emphasize that in this way an outstanding question in the theory how these two effects may occur simultaneously [34] is explained. Altogether, the lowering of magnon energy is quite similar for all the methods although some discrepancies occur. We observe that the SVA is here again insufficient when the magnon momentum has large imbalance between its components k a and k b . Indeed, for  k being close to the X point, the SVA is able to give only half of the magnon softening seen in the VA or in the NA, see figure 3. Good agreement between the VA and the NA found here and in K 2 CuF 4 at = -E J 0.8 z justifies a posteriori the idea of independent determination of orbital angles for the sublattices A and B. The optimal orbital angles q iL for a magnon dressed by orbital excitations are changed in LaMnO 3 by less than ±30°and remain quite similar to the ground state orbitals with θ opt =120°, see figure 4. In general orbital angles increase for the dressed HP bosons. This may be explained because the optimal values of orbital angles q iL follow from the interplay between superexchange interaction and tetragonal crystal fieldE z . The first one favors θ=90°, while the second one favors θ=180°. When a HP boson is created, the spin exchange effectively decreases while the value of E z is not affected.

Conclusions
Summarizing, we have shown that spin-orbital superexchange tunes the orbital angles near local spin excitations and is responsible for novel dressed magnon quasiparticles. The magnon-orbiton coupling is local and reduces nearest neighbor spin exchange J 1 responsible for magnon dispersion at the M point, while orbital fluctuations couple predominantly to spin excitations at neighboring sites and this generates third nearest neighbor J 3 exchange couplings. Thus spin-orbital entanglement has here similar consequences to the exotic phase suggested as a possible ground state of the 2D Kugel-Khomskii model [8]. There the effective spin model derived near a quantum phase transition includes next nearest (J 2 ) and third nearest (J 3 ) neighbor exchange interactions and the latter are of crucial importance to stabilize spin orientation. Here spin-orbital entanglement generates also J 3 interactions which couple spins distant by two lattice constants along the cubic axes a and b, and thus the magnon dispersion is different from that given by the Heisenberg model with nearest neighbor exchange constants derived from frozen orbitals when spin-orbital interactions are disentangled in the ground state. We suggest that such effects would be weaker but still measurable in 3D ordered phases with FM planes, as for instance in LaMnO 3 , and it is very challenging to detect them.
In the electronic model considered here spin-orbital degrees of freedom are entangled and thus respond jointly, giving renormalized spin excitations. However, strong coupling between orbitals and lattice distortions caused by the Jahn-Teller effect will reduce the magnon-orbiton entanglement and thus the renormalization of magnon dispersion reported here will decrease. We suggest that only future experiments could establish importance of spin-orbital entanglement in excited states.
We suggest that similar analysis of the magnon dispersion could be performed using the variational approach for doped FM manganites with statistically averaged interactions between initial S=2 spins at Mn 3+ ions and S=3/2 spins at Mn 4+ ions [57]. We expect that it would reproduce the reduction the magnon energies at the Brillouin zone boundary obtained in the diagrammatic approach [35]. Finally, we remark that the present variational method could also be used to investigate magnon dispersion in the charge, orbital, and spin ordered phase in La 1/2 Sr 3/2 MnO 4 [58].
Appendix. Analytic estimation of the nearest neighbor exchange J 1 Creation of magnons characterized by =  k 0 (being Goldstone modes) does not entail any changes in the orbital background for a spin-orbital system. When magnons characterized by finite   k 0 are created, the coupled orbitals may be slightly modified. As a result, I and  ( ) P k terms deviate from I (0) and To highlight the minute changes due to spin-orbital entanglement, we introduce a vector x consisting of differences in variational parameters with respect to their values in the ground state, and expand I and P in terms of x treated as a small parameter: In the above formulae: (i) in a considered spin-orbital model (3), the  ( ) P k term is factorized into the coefficient T 4 and a  k -dependent term g  k describing the dispersion; (ii) for the sake of clarity the following symbols were introduced: u and U for the gradient and the Hessian of I (x), and w and W for the gradient and the Hessian of T(x), all at x=0.
The above expansions were truncated at quadratic terms, so that the corresponding variational function for a magnon energy has a quadratic form Note that as long as I(x) and T(x) may be expressed analytically, so do u, U, w, W, and, finally, w  k as well. Owing to this, the above formula offers analytically the approximate results without involving any further numerical minimization (as opposed to the strategy used in the main part of the article).
If J 1 is perceived as a parameter in a generic form of a dispersion relation w g = -+¼ , where all other terms that may be introduced for better accounting of the functional dependence such as µ ( ) J k cos 2 a 3 are not written explicitly, then its approximate value may be extracted directly from equation (14) as a coefficient in front of g -S 4 k : where the second term captures the deviation form the frozen orbital description. In order to rationalize the reported value of J 1 for LaMnO 3 we used equation (15) together with the SVA parametrization. The obtained value of the correction is equal to = -´-J J 4.90 10 1 app 3 , in fairly good agreement with the value -´-J 5.22 10 3 resulting from the fit. To avoid lengthy formulae we do not present here a similar approach to determine J 2 and J 3 , and restrict this analytic consideration to the SVA.