Finite temperature magnon spectra in yttrium iron garnet from mean field approach in tight-binding model

We study magnon spectra at finite temperature in yttrium iron garnet from tight-binding model with nearest neighboring exchange interaction. The spin reduction due to thermal magnon excitation are taken into account via the mean field approximation to the local spin and found to be different at two sets of iron atoms. The resulting temperature dependence of the spin wave gap shows good agreement with experiment. We find only two magnon modes are relevant to ferromagnetic resonance.


I. INTRODUCTION
Since its discovery decades ago [1], yttrium iron garnet (YIG) has been regarded as one of the most important magnetic materials due to its extremely low magnetic damping and other intriguing properties [2]. As an insulator, YIG is free from Joule heating, revealing its promising applications in future low energy consumption devices, which makes it a popular material for recent study on magnonics [3,4] and spin caloritronics [5]. While the applications in magnonics mainly rely on the coherent transport properties of spin waves with relatively long wavelength (∼ µm) and low frequency (∼GHz), in spin caloritronic devices, the short wavelength spin waves, namely, magnons, dominate because of their significant population from thermal excitation. In the latter case, recent experiments show that the diffusion length of thermal magnons can reach tens of microns [6,7] and present very interesting behaviors with respect to variation of external magnetic field and temperature [8][9][10] which have not been fully understood so far. Moreover, the magnon diffusion length is also found to be sensitive to the method of measurement [6,7]. The full understanding on those observation obviously requires a comprehensive study of dissipation process beyond the long wavelength limit into a wide range of magnon frequency up to THz.
In the literature, great efforts [11][12][13] have been devoted to explain the origin of magnetic damping in YIG measured by ferromagnetic resonant [14][15][16][17][18][19] or parametric pumping [20,21]. Several mechanisms are believed to be relevant, such as magnon-impurity scattering, magnon-phonon scattering, and magnon-magnon interaction as well as spin-lattice relaxation. However, both the uniform mode in ferromagnetic resonant experiment and finite wave vector magnon from parametric pumping lie in GHz regime, therefore, the applicability of the conclusions from GHz magnons to thermal magnons of THz is questionable. In some recent works, the lifetime of the thermal magnons in YIG were estimated with a single parabolic dispersion [8], which however is valid only for GHz magnons and could depart far from the real spectra. Moreover, at room temperature, several magnon bands could be excited in YIG [22,23], which again reveals the limitation of the single band model.
The full magnon band structure of YIG were modeled by Harris [24] within nearest neighboring exchange interaction and measured from neutron scattering by Plant [22], from which the exchange interaction coefficients were fitted out [2]. The next nearest neighboring exchanges were recently also taken into account with the exchange integrals calculated from first principle calculation [25] or fitted to neutron scattering data [26,27]. Atomistic dynamic simulation shows that the nearest neighboring exchange model with stochastic thermal fluctuation is sufficient to quantitatively reproduce the decrease of the spin wave gap, minimal frequency of the antiferromagnetic-like mode, with increasing temperature [23]. The reason for such a red shift behavior is that the effective magnetic field a local spin felt is suppressed by the reduction of the expectation value of neighboring spins due to magnon excitation. In the low temperature regime, the reduction of average spin is too small to be relevant, resulting in saturate spin wave gap and magnetization, observed in experiment. The gap and magnetization from the atomistic dynamic simulation of classical spin however keep linearly increasing and fails to saturate at very low temperature [23], in contrast to the results from experiment and Harris's quantum model with the same parameters.
In the present work, we develop a mean field approach based on nearest neighboring exchange model, where the magnon excitation at finite temperature is taken into account. In our approach, the temperature dependence of the entire magnon spectra can be easily obtained by selfconsistently solving the magnon population and its influence on dispersion. With temperature independent exchange constants, the resulting spin wave gap and magnetization show good agreement with experiment from zero up to room temperature. In addition to the proper description of the entire spectra beyond the widely used long wavelength approximation, another apparent advantage of our approach is that the output tight-binding type wave functions contain only twenty components, which is therefore convenient for quantitative computation of various spin wave properties, e.g., spin wave dissipation due to different mechanisms, as well as the responses to dif-ferent sources, e.g., light or microwave. Moreover, the information of the polarization, responsible to spin Seebeck effect [23], is also naturally included in the wave functions and easy to read out. Therefore our approach should be useful for study various physics in YIG at finite temperature. As an application, we discuss ferromagnetic resonance and find only two of twenty modes are relevant. Our approach can also be extended to nano structure.

II. EXCHANGE MODEL
In one unit cell in YIG, the collinearly aligning iron atoms at ground state are separated into two sets, eight a atoms and twelve d atoms, according to the polarization or the configuration of neighboring oxygen atoms [2]. The closest distances are r aa = ( √ 3/4)a 0 between two asite atoms, r dd = ( √ 6/8)a 0 between two d-atoms, and r ad = ( √ 5/8)a 0 between a and d atoms with a 0 being the lattice constant of face-centered cubic lattice. Following Harris [24], we start from the nearest neighboring exchange Hamiltonian where J aa , J dd , and J ad represent the corresponding exchange coupling constants. N is the number of unit cell. The last two terms are Zeeman energy. In YIG, all these coupling constants are negative [2]. The coordinates of atoms can be found in Ref. [24].
Without loss of generality, we assume that the equilibrium magnetization follows the small external magnetic field along z axis. By defining transverse components S ± = S x ± iS y , the Holstein-Primakoff transformation of two antiferromagnetically coupled sublattices can be written as [28,29] and S − a,d = (S + a,d ) † where we have taken into account the antiferromagnetic alignment of two sublattices. Thus, a † and d † are the creation operators of clockwise rotation in a-site and anticlockwise rotation in d-site, respectively. Applying Eq. (2) into Eq. (1) and Fourier transform a in = 1 √ N k a i (k)e ik · Rin , one obtains the tight-binding type Bogoliubov-de Gennes Hamiltonian for an arbitrary wave vector k (the quadratic order) [2,24] where with γ ηη ij (k) = |rij |=r ηη e ik · rij . The last term in Eq. (3) leads to angular momentum exchange between the two sublattices preserving the chirality.

III. RESULTS
In this section, we present our results from the diagonalization of the above Bogoliubov-de Gennes Hamiltonian. The exchange constants are specified to be Both a-and d-sites are approximately of S a = S d = 5/2. Note that, instead of unitary transform in fermionic system, a paraunitary transform is needed in the diagonalization procedure of bosonic Hamiltonian matrix under basis (a k , d k , a † −k , d † −k ) [30,31]. The resulting eigenstates contains two sets which satisfy the Boson commutation rule guaranteed by the paraunitary transform. Einstein summation convention has been applied with i, i = 1, ..., 8 and j, j = 1, ..., 12. Since the creation of an anticlockwise rotation of d-site (d † ) is equivalent to the annihilation of a clockwise rotation, α ik are annihilation operators of pure clockwise rotation. In the same sense, β jk are about purely anticlockwise rotation. The dynamics of the transverse components of local spin for two types of modes read where the time evolution of the coefficients are c ii 1, The inverse transform gives the original site-operators in form of a combination of magnon operators of eigenstates The fluctuation of local spins is then given by which shows that at very low temperature, where the thermal excitation is negligible, i.e., α † ik α ik ≈ β † jk β jk ≈ 0, there is still non-zero fluctuation due to the coupling between a-and d-site rotations via B ij terms in Eq. (3) [32,33].
A. Zero temperature spectra Fig. 1(a) shows all twenty magnon branches in [110] and [100] directions carried out from zero temperature in the absence of external magnetic field. The characteristics of polarization α ik and β ik are distinct from the color red and blue, respectively. Typically, the eight clockwise modes have higher energy than the twelve anticlockwise modes, since the former correspond to the a-site spins rotating under the exchange field of six nearby d-site spins, which is stronger (∼ 12|J ad |S d ) than the one felt by the d-site spin from the four nearby a-site spins (∼ 8|J ad |S a ). This is confirmed by Fig. 1(b) from the calculation with the same parameters but setting all B ij = 0, therefore, no coupling between clockwise and anticlockwise rotations. From comparison between the results with and without B ij terms, we can see that the angular momentum exchange between the sublattices partially lifts the degeneracy at high symmetry points, and more importantly leads to the correct dispersion of acoustic mode and the antiferromagnetic-like mode, labeled as β 1 and α 1 , separately. The β 1 mode is almost isotropic up to 5 THz and starts to become anisotropic approaching the boundary of the Brillouin zone.
In Fig. 1(c) we plot the instantaneous spin configuration of each mode at the center of Brillouin zone (k = 0). The transverse spins in both α 1 and β 1 are parallel within the sublattices and antiparallel between sublattices. The magnitudes of spins in β 1 mode are all the same, as a result, this tilted configuration is a ground state of the system and stable, explaining its zero frequency. The a spins in α 1 contain larger transverse spin components and therefore dominate the rotation direction of this mode. Without B ij , the two sublattices are decoupled, as a result, the parallel configuration of each sublattice (α 1 and β 1 ) is of highest energy among the separately clockwise and anticlockwise modes due to negative exchange coupling constant between the same type of atoms [see Fig. 1(b)].

B. Finite temperature spectra
At low temperature, the density of the thermally excited magnon is low enough so that the magnon dispersion calculated at zero temperature still works well. However as more magnons are excited at high temperatures, the enhanced magnon-magnon interaction can strongly affect the magnon dispersion. Physically, for example, the large magnon density reduces the magnitude of the expectation value of all local spins, hence the exchange field felt by the neighboring spins. The effect of magnonmagnon interaction can be taken into account by includ- ing the higher order terms beyond quadratic approximation, e.g., those in form of The modification in magnon spectra then can be carried out by diagrammatic calculation on the real part of the magnon self-energy due to these interaction terms.
Here, alternatively, we perform a mean field approximation with a self-consistent procedure, whose advantage has been discussed in the Introduction.
Specifically, we keep all magnon-magnon interactions with up to three magnon operators at the same site and apply for interactions between a-site spins to reduce the interaction terms back to quadratic order. The resulting effective quadratic Hamiltonian is exactly of the same expression as Eq. (3) after a replacement of S a(d) by where ∆S a(d) stands for the average of thermal excitation over all a(d) site spins Here, n B corresponds to Plank distribution at thermal equilibrium. By diagonalizing Eq. (3) with S eff a(d) , we obtain modified spin wave spectra, from which we compute updated S eff a(d) and repeat such process until reaching a self-consistent solution.
In Figs. 2(a) and (b), we plot the convergent magnon spectra at 100 K and 300 K, respectively. As we can see, the shape of each band roughly remain the same at different temperatures. However, as the temperature increases, the clockwise rotating mode (in red) move towards to lower frequency regime, which is consistent with experimental observation [22] and atomistic simulation [23]. The anticlockwise rotating modes are relatively insensitive to the temperature change.
In Fig. 2(c), the effective magnitudes of two sets spins, i.e., Eq. (19), are plotted as function of temperature, which shows weak dependence in the low temperature regime, and an approximately linear decrease in the temperature regime. Interestingly, the reduction at d-site is much larger than that at a-site. The reason for such behavior is that the d site has larger components than a site in the low energy magnon bands, who are efficiently excited. At room temperature, the effective spin reduces to 2.21 and 2.03 respectively. These distinct reductions are actually the origin of the relative shift between anticlockwise and clockwise sets in Figs. 2(a) and (b), because the frequencies of the clockwise modes rely on the magnitude of d spins. If the two effective magnitudes are equal S eff a = S eff d , all frequencies would vary simultaneously by a global prefactor according to Eq. (3), which is proportional to S eff a = S eff d at vanishing magnetic field. We summarize in Fig. 2(d) the frequency at the button of the α 1 mode (i.e., the spin wave gap mentioned in the Introduction) at different temperatures as black solid squares. For comparison, the experimental data [22] and results from atomistic simulation [23] are plotted in the same figure as red bullet and open squares, respectively. As we can see that our results show good agreement with experiment especially at low temperature, where limited magnon are excited hence approaching to the zero temperature case. Our approach therefore avoids the low temperature deviation in atomistic simulation based on classical spin. In the inset, we use the magnetization from our calculation at finite temperature and plot the frequency as function of the corresponding magnetization of the same temperature, which presents good linear dependence. The magnetization at room temperature is around µ 0 M s (300 K) = 0.163 T , being reasonable agreement with experimetal data 0.175 T [34]. The deviation from experimental value is because of the fact that the temperature dependence of the magnon spectra was not taken into account in the procedure to determine the exchange constants from temperature dependence of magnetization [2]. The determination of the correction in exchange constants due to the change in magnon spectra is beyond the scope of the present work.

IV. FERROMAGNETIC RESONANCE
As discussed in the Introduction, with the knowledge of wave functions, one can calculate various properties. Here, we take ferromagnetic resonance for an example. Assume spatially uniform ac magnetic field transverse to the magnetization, we write out the time dependent perturbation Hamiltonian to describe its coupling with magnetic atoms For linear polarization magnetic field along x-direction, with G α i = i C i i 1,0 + j C j i * 3,0 , G β j = i C i j * 2,0 + j C j j 4,0 , and B x (t) = B 0 cos(ωt). Here, only k = 0 survives in the space integration. The absorption rates are then given by Fermi Golden rule which is proportional to the number of unit cells in magnetic sample and the input power (∝ B 2 0 ). By using the coefficients in wave function C 1,2,3,4 , we find only two of the prefactors G α 1 and G β 1 are non-zero. This actually can be easily understood from the spin configuration in Fig. 1(c), where only α 1 and β 1 modes have net spin. The vanishing of net spin in other modes makes them transparent to driving field.

V. CONCLUSION AND DISCUSSION
In conclusion, we calculated the entire magnon spectra at finite temperature within mean field approximation on the magnon-magnon interaction. The temperature dependence of magnon spectra shows good agreement with experiment. We find the reduction in magnitude of d spin is larger than that of a atoms. By using the wave functions of eigenstates, we analyze the ferromagnetic resonance and find that only two modes can be drived by spatially uniform ac magnetic field. We suppose our approach which provides proper magnon spectra and corresponding wave functions can be usefully to study various properties in yttrium iron garnet at finite temperature as discussed in the Introduction.
One may notice that in Fig. 2 the residual magnetization (effective local spin) near Curie temperature (560 K) is still finite, which reveals the limitation of the spin wave approximation in the high temperature regime [2]. In the literature, the mean field theory based on sublattices' magnetizations was found to perform well up to Curie temperature [35]. How to link this theory to our tightbinding type mean field approach and how to extend the present spin-wave-based model up to Curie temperature in ferrimagnetic systems are still under study.