Finite temperature magnon spectra in yttrium iron garnet from mean field approach in tight-binding model
Ka Shen
cond-mat.mes-hall








[],n















↑
↓

    ξ∇ rẑ p V v j k q A E P R S T u vϵξ𝒜 A[2]
        The Center for Advanced Quantum Studies and Department of Physics, Beijing Normal University, Beijing 100875, China
        
        
        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.
        
        
        72.25.Dc, 75.70.Tj, 85.75.-dFinite temperature magnon spectra in yttrium iron garnet from mean field approach in tight-binding model
            Ka Shen
            December 30, 2023
        ========================================================================================================
        
        
        
        § INTRODUCTION
        
        Since its discovery decades ago <cit.>, 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 <cit.>. 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 <cit.> and spin caloritronics <cit.>. 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 <cit.> and present very interesting behaviors with respect to variation of external magnetic field and temperature <cit.> which have not been fully understood so far. Moreover, the magnon diffusion length is also found to be sensitive to the method of measurement <cit.>. 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 <cit.> have been devoted to explain the origin of magnetic damping in YIG measured by ferromagnetic resonant <cit.> or parametric pumping <cit.>. 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 <cit.>, 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 <cit.>, which again reveals the limitation of the single band model.
        
        The full magnon band structure of YIG were modeled by Harris <cit.> within nearest neighboring exchange interaction and measured from neutron scattering by Plant <cit.>, from which the exchange interaction coefficients were fitted out <cit.>. The next nearest neighboring exchanges were recently also taken into account with the exchange integrals calculated from first principle calculation <cit.> or fitted to neutron scattering data <cit.>. 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 <cit.>. 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 <cit.>, 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 different sources, e.g., light or microwave. Moreover, the information of the polarization, responsible to spin Seebeck effect <cit.>, 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.
        
        
        
        § 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 <cit.>. The closest distances are r_aa=(√(3)/4)a_0 between two a-site 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 <cit.>, we start from the nearest neighboring exchange Hamiltonian
        
        H   =   -∑_n=1^N[J_aa∑_i^8∑_|𝐫_ij|=r_aa𝐒_a(𝐑_in)·𝐒_a(𝐑_in+𝐫_ij)
              +2J_ad∑_i=1^8∑_|𝐫_ij|=r_ad𝐒_a(𝐑_in)·𝐒_d(𝐑_in+𝐫_ij)
              +J_dd∑_i=9^20∑_|𝐫_ij|=r_dd𝐒_d(𝐑_in)·𝐒_d(𝐑_in+𝐫_ij)
              +∑_i=1^8gμ_B𝐁·𝐒_a(𝐑_in)+∑_i=9^20gμ_B𝐁·𝐒_d(𝐑_in)], 
        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 <cit.>. The coordinates of atoms can be found in Ref. <cit.>.
        
        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 <cit.>
        S_a^z=S_a-a^†a,    S_a^+=√(2S_a-a^†a)a,
        
        S_d^z=-S_d+d^†d,    S_d^+=d^†√(2S_d-d^†d),
        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. (<ref>) into Eq. (<ref>) and Fourier transform a_in=1/√(N)∑_𝐤a_i(𝐤)e^i𝐤·𝐑_in, one obtains the tight-binding type Bogoliubov-de Gennes Hamiltonian for an arbitrary wave vector 𝐤 (the quadratic order) <cit.>
        H_𝐤   =   ∑_i,j=1^8A_ij(𝐤)a_i^†(𝐤)a_j(𝐤)+∑_i,j=9^20D_ij(𝐤)d_i^†(𝐤)d_j(𝐤)
              +∑_i=1^8∑_j=9^20[B_ij(𝐤)a_i^†(𝐤)d_j^†(-𝐤)+h.c.], 
        where
        
        A_ij(𝐤)   =   (-12J_adS_d+16J_aaS_a+gμ_BB)δ_ij
              -2J_aaS_aγ_ij^aa(𝐤),
        
        D_ij(𝐤)   =   (-8J_adS_a+8J_ddS_d-gμ_BB)δ_ij
              -2J_ddS_dγ_ij^dd(𝐤),
        
        B_ij(𝐤)   =   -2J_ad√(S_aS_d)γ_ij^ad(𝐤),
        
        with γ_ij^ηη'(𝐤)=∑_|r_ij|=r_ηη'e^i𝐤·𝐫_ij. The last term in Eq. (<ref>) leads to angular momentum exchange between the two sublattices preserving the chirality.
        
        
        
        § 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 J_ad=-39.8 K, J_dd=-13.4 K, and J_aa=-3.8 K <cit.>. 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_𝐤,d_𝐤,a^†_-𝐤,d^†_-𝐤) <cit.>. The resulting eigenstates contains two sets
        α_ik   =   c_1,𝐤^ii'a_i'𝐤+c_2,𝐤^ij'd_j'-𝐤^†,
        β_jk   =   c_3,𝐤^ji'a_i'-𝐤^†+c_4,𝐤^jj'd_j'𝐤,
        
        which satisfy the Boson commutation rule 
           [α_i𝐤,α^†_i'𝐤']=δ_ii'δ_𝐤𝐤',
           [β_j𝐤,β^†_j'𝐤']=δ_jj'δ_𝐤𝐤',
        
        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, α_i𝐤 are annihilation operators of pure clockwise rotation. In the same sense, β_j𝐤 are about purely anticlockwise rotation. The dynamics of the transverse components of local spin for two types of modes read
        
        (S_i'a^x,S_i'a^y)_α_ik(t)   =   ( c_1,𝐤^ii'(t), c_1,𝐤^ii'(t))
        
        (S_j'd^x,S_j'd^y)_α_ik(t)   =   ( c_2,𝐤^ij'(t), c_2,𝐤^ij'(t))
        
        (S_i'a^x,S_i'a^y)_β_jk(t)   =   ( c_3,𝐤^ji'(t),- c_3,𝐤^ji'(t))
        
        (S_j'd^x,S_j'd^y)_β_jk(t)   =   ( c_4,𝐤^ij'(t),- c_4,𝐤^ij'(t))
        
        where the time evolution of the coefficients are
        c_1,𝐤^ii'(t)=e^ik·r_i'-iω_α_iktc_1,𝐤^ii',
        c_2,𝐤^ij'(t)=e^-ik·r_j'-iω_α_iktc_2,𝐤^ij',
        c_3,𝐤^ji'(t)=e^-ik·r_i'-iω_β_jktc_3,𝐤^ji',
        c_4,𝐤^ij'(t)=e^ik·r_j'-iω_β_jktc_4,𝐤^ij'.
        The inverse transform gives the original site-operators in form of a combination of magnon operators of eigenstates
        
        a_i'k   =   C_1,𝐤^i'iα_i𝐤+C_2,-𝐤^i'jβ_j-𝐤^†,
        
        d_j'k   =   C_3,-𝐤^j'iα_i-𝐤^†+C_4,𝐤^j'jβ_j𝐤.
        
        The fluctuation of local spins is then given by 
        ⟨  a_i'^† a_i'⟩   =   1/N∑_𝐤|C_1,𝐤^i'i|^2 ⟨α_i𝐤^†α_i𝐤⟩ + |C_2,𝐤^i'j|^2 (⟨β_j𝐤^†β_j𝐤⟩+1),
        
        ⟨  d_j'^† d_j'⟩   =   1/N∑_𝐤|C_3,𝐤^j'i|^2 (⟨α_i𝐤^†α_i𝐤⟩ +1)
        + |C_4,𝐤^j'j|^2 ⟨β_j𝐤^†β_j𝐤⟩,
        
        which shows that at very low temperature, where the thermal excitation is negligible, i.e., ⟨α_i𝐤^†α_i𝐤⟩≈⟨β_j𝐤^†β_j𝐤⟩≈0, there is still non-zero fluctuation due to the coupling between a- and d-site rotations via B_ij terms in Eq. (<ref>) <cit.>.
        
        
        
        
         §.§ Zero temperature spectra
        
        Fig. <ref>(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 α_i𝐤 and β_i𝐤 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. <ref>(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. <ref>(c) we plot the instantaneous spin configuration of each mode at the center of Brillouin zone (𝐤=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. <ref>(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 magnon-magnon interaction can be taken into account by including the higher order terms beyond quadratic approximation, e.g., those in form of a_i^† a_i^† a_i a_j,  a_i^† a^†_i a_i d_j, a_i^† a_i d^†_j d_j,  etc. 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
         
        a^†_i a_i→⟨  a^†_i a_i⟩,      
        a^†_i a_i^† a_i→ 2⟨  a^†_i a_i ⟩ a^†_i
        
        d^†_j d_j→⟨  d^†_j d_j⟩ ,      
        d^†_j d_j^† d_j→ 2⟨  d^†_j d_j ⟩ d^†_j
        
        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. (<ref>) after a replacement of S_a(d) by
        
        S_a(d)^ eff=S_a(d)-Δ S_a(d). 
        where Δ S_a(d) stands for the average of thermal excitation over all a(d) site spins
        Δ S_i'a   =   1/N∑_𝐤|C_1,𝐤^i'i|^2n_B(ħω_α_i𝐤)+|C_2,𝐤^i'j|^2n_B(ħω_β_i-𝐤),
        Δ S_j'd   =   1/N∑_𝐤|C_3,𝐤^j'i|^2n_B(ħω_α_i-𝐤)+|C_4,𝐤^j'j|^2n_B(ħω_β_i𝐤).
        
        Here, n_B corresponds to Plank distribution at thermal equilibrium. By diagonalizing Eq. (<ref>) 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. <ref>(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 <cit.> and atomistic simulation <cit.>. The anticlockwise rotating modes are relatively insensitive to the temperature change.
        
        In Fig. <ref>(c), the effective magnitudes of two sets spins, i.e., Eq. (<ref>), 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. <ref>(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. (<ref>), which is proportional to S_a^ eff=S_d^ eff at vanishing magnetic field.
        
        
        We summarize in Fig. <ref>(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 <cit.> and results from atomistic simulation <cit.> 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 μ_0M_s(300 K)=0.163 T, being reasonable agreement with experimetal data 0.175 T <cit.>. 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 <cit.>. The determination of the correction in exchange constants due to the change in magnon spectra is beyond the scope of the present work.
        
        
        
        § 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
        
        H'(t)=gμ_B 𝐁(t)· [∑_i,n𝐒_a(𝐑_in)+∑_j,n𝐒_d(𝐑_jn)].
        
        For linear polarization magnetic field along x-direction, 
        
        H'(t)=√(N)gμ_BB_x(t)( G^α_iα_i0+ G^β_jβ_j0)+h.c.
        
        with G_i^α=∑_i'C_1,0^i'i+∑_j'C_3,0^j'i∗, G_j^β=∑_i'C_2,0^i'j∗+∑_j'C_4,0^j'j, and B_x(t)=B_0cos(ω t). Here, only 𝐤=0 survives in the space integration. The absorption rates are then given by Fermi Golden rule
        Γ(ω)   =   π/2ħ(gμ_BB_0)^2N {| G_i^α|^2[N_B(ω_α_i0)+1]δ(ω-ω_α_i0)
              +| G_j^β|^2[N_B(ω_β_j0)+1]δ(ω-ω_β_j0)},
        
        which is proportional to the number of unit cells in magnetic sample and the input power (∝ B_0^2). 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. <ref>(c), where only α_1 and β_1 modes have net spin. The vanishing of net spin in other modes makes them transparent to driving field.
        
        
        
        § 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. <ref> 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 <cit.>. In the literature, the mean field theory based on sublattices' magnetizations was found to perform well up to Curie temperature <cit.>. How to link this theory to our tight-binding 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.
        
        
        The author thanks G. E. W. Bauer and J. L. Cheng for valuable discussions.
        This work is supported by the Recruitment Program of Global Youth Experts. 
        
        prsty
