Magnon-phonon interactions in magnetic insulators

TL;DR

This paper develops a theoretical framework for magnon-phonon interactions in magnetic insulators, specifically yttrium iron garnet, analyzing their impact on quasi-particle and transport lifetimes at intermediate temperatures.

Contribution

It introduces a detailed theoretical analysis of magnon-phonon interactions, including anisotropy- and exchange-mediated mechanisms, and evaluates their effects on magnon transport and damping.

Findings

Phonon scattering has weak effects on magnon transport.

Magnon transport lifetime differs from quasi-particle lifetime at short wavelengths.

The study provides insights into magnon damping and transport in yttrium iron garnet.

Abstract

We address the theory of magnon-phonon interactions and compute the corresponding quasi-particle and transport lifetimes in magnetic insulators with focus on yttrium iron garnet at intermediate temperatures from anisotropy- and exchange-mediated magnon-phonon interactions, the latter being derived from the volume dependence of the Curie temperature. We find in general weak effects of phonon scattering on magnon transport and the Gilbert damping of the macrospin Kittel mode. The magnon transport lifetime differs from the quasi-particle lifetime at shorter wavelengths.

Equations

$${\mathcal{H}}_m=-\frac{J}{2}\sum _{⟨i\ne j⟩}{\mathbf{S}}_i\cdot {\mathbf{S}}_j-\sum _{i}g{\mu }_{B}B{S}_i^z,$$

$$S_i^{+}$$

$$S_i^{-}$$

$$S_i^z$$

$${\mathcal{H}}_m\to \sum _{\mathbf{k}}{E}_{\mathbf{k}}{b}_{\mathbf{k}}^{†}{b}_{\mathbf{k}},$$

$$b_i$$

$$b_i^{†}$$

$$E_{\mathbf{k}}=g{\mu }_{B}B+4SJ\sum _{\alpha =x,y,z}{\sin }^2(k_{\alpha }a/2)$$

$$E_{\mathbf{k}}=g{\mu }_{B}B+E_{ex}{k}^2{a}^2,$$

$${\mathbf{X}}_i={\mathbf{r}}_i-{\mathbf{R}}_i,$$

$$X_i^{\alpha }=\frac{1}{\sqrt{N}}\sum _{\mathbf{q},\lambda }{e}_{\mathbf{q}\lambda }^{\alpha }{X}_{\mathbf{q}\lambda }{e}^{i\mathbf{q}\cdot {\mathbf{R}}_i},$$

$${\mathbf{e}}_{\mathbf{q}1}$$

$${\mathbf{e}}_{\mathbf{q}2}$$

$${\mathbf{e}}_{\mathbf{q}3}$$

$$\mathbf{q}=q\left(\sin {\theta }_{\mathbf{q}}\cos {\varphi }_{\mathbf{q}},\sin {\theta }_{\mathbf{q}}\sin {\varphi }_{\mathbf{q}},\cos {\theta }_{\mathbf{q}}\right),$$

$${\mathcal{H}}_p$$

$$=\sum _{\mathbf{q}\lambda }{\epsilon }_{\mathbf{q}\lambda }\left(a_{\mathbf{q}\lambda }^{†}{a}_{\mathbf{q}\lambda }+\frac{1}{2}\right),$$

$${\epsilon }_{\mathbf{q}\lambda }=\hslash c_{\lambda }|\mathbf{q}|,$$

$$X_{\mathbf{q}\lambda }=\frac{a_{\mathbf{q}\lambda }+a_{-\mathbf{q}\lambda }^{†}}{\sqrt{2m{\epsilon }_{\mathbf{q}\lambda }/{\hslash }^2}},P_{\mathbf{q}\lambda }=\frac{1}{i}\sqrt{\frac{m{\epsilon }_{\mathbf{q}\lambda }}{2}}\left(a_{\mathbf{q}\lambda }-a_{-\mathbf{q}\lambda }^{†}\right),$$

$$E_{me}=$$

$$+B_{\alpha \beta }^{\prime }\frac{\partial \mathbf{M}(\mathbf{r})}{\partial r_{\alpha }}\cdot \frac{\partial \mathbf{M}(\mathbf{r})}{\partial r_{\beta }}]X_{\alpha \beta }(\mathbf{r}),$$

$$X_{\alpha \beta }(\mathbf{r})=\frac{1}{2}\left[\frac{\partial X_{\alpha }(\mathbf{r})}{\partial r_{\beta }}+\frac{\partial X_{\beta }(\mathbf{r})}{\partial r_{\alpha }}\right],$$

$$B_{\alpha \beta }$$

$$B_{\alpha \beta }^{\prime }$$

$$B_{\parallel }$$

$$B_{\perp }$$

$$B_{\parallel }^{\prime }/a^2$$

$$B_{\perp }^{\prime }/a^2$$

$${\mathcal{H}}_{mp}^{an}$$

$$+\frac{1}{\sqrt{N}}\sum _{\mathbf{q},\mathbf{k},{\mathbf{k}}^{\prime }}{\delta }_{\mathbf{k}-{\mathbf{k}}^{\prime }-\mathbf{q},0}\sum _{\lambda }{\Gamma }_{{\mathbf{kk}}^{\prime },\lambda }^{an}{b}_{\mathbf{k}}^{†}{b}_{{\mathbf{k}}^{\prime }}{X}_{\mathbf{q}\lambda }$$

$$+\frac{1}{\sqrt{N}}\sum _{\mathbf{q},\mathbf{k},{\mathbf{k}}^{\prime }}{\delta }_{\mathbf{k}+{\mathbf{k}}^{\prime }+\mathbf{q},0}\sum _{\lambda }{\Gamma }_{{\mathbf{kk}}^{\prime },\lambda }^{bb}{b}_{\mathbf{k}}{b}_{{\mathbf{k}}^{\prime }}{X}_{\mathbf{q}\lambda }$$

$$+\frac{1}{\sqrt{N}}\sum _{\mathbf{q},\mathbf{k},{\mathbf{k}}^{\prime }}{\delta }_{\mathbf{k}+{\mathbf{k}}^{\prime }-\mathbf{q},0}\sum _{\lambda }{\Gamma }_{{\mathbf{kk}}^{\prime },\lambda }^{\bar{b}\bar{b}}{b}_{\mathbf{k}}^{†}{b}_{{\mathbf{k}}^{\prime }}^{†}{X}_{\mathbf{q}\lambda },$$

$${\Gamma }_{\mathbf{q}\lambda }$$

$$+\left(i{q}_x+q_y\right)e_{\mathbf{q}\lambda }^z],$$

$${\Gamma }_{{\mathbf{kk}}^{\prime },\lambda }^{an}$$

$${\Gamma }_{{\mathbf{kk}}^{\prime },\lambda }^{bb}$$

$${\Gamma }_{{\mathbf{kk}}^{\prime },\lambda }^{\bar{b}\bar{b}}$$

$$U_{\mathbf{q},\lambda }$$

$$V_{\mathbf{q},\lambda }$$

$$+\frac{B_{\perp }}{S}\left[q_y{e}_{\mathbf{q}\lambda }^x+q_x{e}_{\mathbf{q}\lambda }^y\right].$$

Full text

I Introduction

Magnons are the elementary excitations of magnetic order, i.e. the quanta of spin waves. They are bosonic and carry spin angular momentum. Of particular interest are the magnon transport properties in yttrium iron garnet (YIG) due to its very low damping ($\alpha<10^{-4}$), which makes it one of the best materials to study spin-wave or spin caloritronic phenomena [1, 2, 3, 4, 5, 6]. For instance, the spin Seebeck effect (SSE) in YIG has been intensely studied in the past decade [7, 8, 9, 10, 11, 12, 13]. Here, a temperature gradient in the magnetic insulator injects a spin current into attached Pt contacts that is converted into a transverse voltage by the inverse spin Hall effect. Most theories explain the effect by thermally induced magnons and their transport to and through the interface to Pt [7, 14, 15, 16, 17, 18, 19]. However, phonons also play an important role in the SSE through their interactions with magnons [20, 21, 22].

Magnetoelastic effects in magnetic insulators were addressed first by Abrahams and Kittel [23, 24, 25], and by Kaganov and Tsukernik [26]. In the long-wavelength regime, the strain-induced magnetic anisotropy is the most important contribution to the magnetoelastic energy, whereas for shorter wavelengths, the contribution from the strain-dependence of the exchange interaction becomes significant [27, 28, 29]. Rückriegel *et al. *[28] computed very small magnon decay rates in thin YIG films due to magnon-phonon interactions with quasi-particle lifetimes $\tau_{qp}\apprge 480\;\mathrm{ns,}$even at room temperature. However, these authors do not consider the exchange interaction and the difference between quasi-particle and transport lifetimes.

Recently, it has been suggested that magnon spin transport in YIG at room temperature is driven by the magnon chemical potential [3, 30]. Cornelissen et al.  [3] assume that at room temperature magnon-phonon scattering of short-wavelength thermal magnons is dominated by the exchange interaction with a scattering time of $\tau_{qp}\sim 1\;\mathrm{ps}$, which is much faster than the anisotropy-mediated magnon-phonon coupling considered in Ref. [28] and efficiently thermalizes magnons and phonons to equal temperatures without magnon decay. Recently, the exchange-mediated magnon-phonon interaction [31] has been taken into account in a Boltzmann approach to the SSE, but this work underestimates the coupling strength by an order of magnitude, as we will argue below.

In this paper we present an analytical and numerical study of magnon-phonon interactions in bulk ferromagnetic insulators, where we take both the anisotropy- and the exchange-mediated magnon-phonon interactions into account. By using diagrammatic perturbation theory to calculate the magnon self-energy, we arrive at a wave-vector dependent expression of the magnon scattering rate, which is the inverse of the magnon quasi-particle lifetime $\tau_{qp}$. The magnetic Grüneisen parameter $\Gamma_{m}=\partial\ln T_{C}/\partial\ln V$ [32, 33], where $T_{C}$ is the Curie temperature and $V$ the volume of the magnet, gives direct access to the exchange-mediated magnon-phonon interaction parameter. We predict an enhancement in the phonon scattering of the Kittel mode at the touching points of the two-magnon energy (of the Kittel mode and a finite momentum magnon) and the longitudinal and transverse phonon dispersions, for YIG at around $1.3\;\mathrm{T}$ and $4.6\;\mathrm{T}$. We also emphasize the difference in magnon lifetimes that broaden light and neutron scattering experiments, and the transport lifetimes that govern magnon heat and spin transport.

The paper is organized as follows: in Sec. II we briefly review the theory of acoustic magnons and phonons in ferro-/ferrimagnets, particularly in YIG. In Sec. III we derive the exchange- and anisotropy-mediated magnon-phonon interactions for a cubic Heisenberg ferromagnet with nearest neighbor exchange interactions in the long-wavelength limit. In Sec. IV we derive the magnon decay rate from the imaginary part of the magnon self-energy in a diagrammatic approach and in Sec. V we explain the differences between the magnon quasi-particle and transport lifetimes. Our numerical results for YIG are discussed in Sec. VI. Finally in Sec. VII we summarize and discuss the main results of the present work. The validity of our long-wavelength approximation is analyzed in Appendix A and in Appendix B we explain why second order magnetoelastic couplings may be disregarded. In Appendix C we briefly discuss the numerical methods used to evaluate the k-space integrals.

II Magnons and phonons in ferromagnetic insulators

Without loss of generality, we focus our treatment on yttrium iron garnet (YIG). The magnon band structure of YIG has been determined by inelastic neutron scattering [34, 35, 36] and by ab initio calculation of the exchange constants [37]. The complete magnon spectral function has been computed for all temperatures by atomistic spin simulations [38], taking all magnon-magnon interactions into account, but not the magnon-phonon scattering. The pure phonon dispersion is known as well [39, 29]. In the following, we consider the interactions of the acoustic magnons from the lowest magnon band with transverse and longitudinal acoustic phonons, which allows a semi-analytic treatment but limits the validity of our results to temperatures below $100\;\mathrm{K}$. Since the low-temperature values of the magnetoelastic constants, sound velocities, and magnetic Grüneisen parameter are not available for YIG, we use throughout the material parameters under ambient conditions.

II.1 Magnons

Spins interact with each other via dipolar and exchange interactions. We disregard the former since at the energy scale $E_{\mathrm{dip}}\approx 0.02\;\mathrm{meV}$ [28] it is only relevant for long-wavelength magnons with wave vectors $k\lesssim 6\times 10^{7}\;\mathrm{m}^{-1}$ and energies $E_{\mathbf{k}}/k_{B}\lesssim 0.2\;\mathrm{K}$, which are negligible for the thermal magnon transport in the temperature regime we are interested in. The lowest magnon band can then be described by a simple Heisenberg model on a course-grained simple cubic ferromagnet with exchange interaction $J$

[TABLE]

where the sum is over all nearest neighbors and $\hbar\mathbf{S}_{i}$ is the spin operator at lattice site $\mathbf{R}_{i}$. The lattice constant of the cubic lattice or YIG is $a=12.376\;\mathrm{\mathring{A}}$ and the effective spin per unit cell $\hbar S=\hbar M_{s}a^{3}/(g\mu_{B})\approx 14.2\hbar$ at room temperature [28] ($S\approx 20$ for $T\lesssim 50\;\mathrm{K}$ [40]), where the g-factor $g\approx 2$, $\mu_{B}$ is the Bohr magneton and $M_{s}$ the saturation magnetization. The parameter $J$ is an adjustable parameter that can be fitted to experiments or computed from first principles. $B$ is an effective magnetic field that orients the ground-state magnetization vector to the $z$ axis and includes the (for YIG small) magnetocrystalline anisotropy field. The $1/S$ expansion of the spin operators in terms of Holstein-Primakoff bosons reads [41],

[TABLE]

where $b_{i}^{\dagger}$ and $b_{i}$ are the magnon creation and annihilation operators with boson commutation rule $\left[b_{i},b_{j}^{\dagger}\right]=\delta_{i,j}$. Then

[TABLE]

where the magnon operators $b_{\mathbf{k}}^{\dagger}$ and $b_{\mathbf{k}}$ are defined by

[TABLE]

and $N$ the number of unit cells. The dispersion relation

[TABLE]

becomes quadratic in the long-wavelength limit $ka\ll 1$:

[TABLE]

where $E_{ex}=SJ$. With $E_{ex}=k_{B}\times 40\;\mathrm{K}=3.45\;\mathrm{meV}$ the latter is a good approximation up to $k_{0}=1/a\approx 8\times 10^{8}\;\mathrm{m}^{-1}$ [34]. The effective exchange coupling is then $J\approx 0.24\;\mathrm{meV}$. The lowest magnon band does not depend significantly on temperature [38], which implies that $E_{ex}=SJ$ does not depend strongly on temperature. The temperature dependence of the saturation magnetization and effective spin $S$ should therefore not affect the low-energy exchange magnons significantly. By using Eq. (9) in the following, our theory is valid for $k\lesssim k_{0}$ (see Fig. 1) or temperatures $T\lesssim 100\;\mathrm{K}$. In this regime the cut-off of an ultraviolet divergence does not affect results significantly (see Appendix A). We disregard magnetostatic interactions that affect the magnon spectrum only for very small wave vectors since at low temperatures the phonon scattering is not significant.

II.2 Phonons

We expand the displacement $\mathbf{X}_{i}$ of the position $\mathbf{r}_{i}$ of unit cell $i$ from the equilibrium position $\mathbf{R}_{i}$

[TABLE]

into the phonon eigenmodes $X_{\mathbf{q}\lambda}$,

[TABLE]

where $\alpha\in\left\{x,y,z\right\}$ and $\mathbf{q}$ a wave vector. We define polarizations $\lambda\in\left\{1,2,3\right\}$ for the elastic continuum [42]

[TABLE]

where the angles $\theta_{\mathbf{q}}$ and $\phi_{\mathbf{q}}$ are the spherical coordinates of

[TABLE]

which is valid for YIG up to $3\;\mathrm{THz}$ ($12\;\mathrm{meV}$) [39, 29]. The phonon Hamiltonian then reads

[TABLE]

where the canonical momenta $P_{\mathbf{q}\lambda}$ obey the commutation relations $\left[X_{\mathbf{q}\lambda},P_{\mathbf{q}^{\prime}\lambda^{\prime}}\right]=i\hbar\delta_{\mathbf{q},-\mathbf{q}^{\prime}}\delta_{\lambda\lambda^{\prime}}$ and the mass of the YIG unit cell $m=\rho a^{3}=9.8\times 10^{-24}\;\mathrm{kg}$ [27]. The phonon dispersions for YIG then read

[TABLE]

where $c_{1,2}=c_{t}=3843\;\mathrm{m/s}$ is the transverse sound velocity and $c_{3}=c_{l}=7209\;\mathrm{m/s}$ the longitudinal velocity at room temperature [27]. In terms of phonon creation and annihilation operators

[TABLE]

and $\left[a_{\mathbf{q}\lambda},a_{\mathbf{q}^{\prime}\lambda^{\prime}}^{\dagger}\right]=\delta_{\mathbf{q},\mathbf{q}^{\prime}}\delta_{\lambda,\lambda^{\prime}}$.

In Fig. 1 we plot the longitudinal and transverse phonon and the acoustic magnon dispersion relations for YIG at zero magnetic field. The magnon-phonon interaction leads to an avoided level crossing at points where magnon and phonon dispersion cross, as discussed in Refs. [27] and [28].

III Magnon-phonon interactions

We derive in this section the magnon-phonon interactions due to the anisotropy and exchange interactions for a cubic lattice ferromagnet.

III.1 Phenomenological magnon-phonon interaction

In the long-wavelength/continuum limit ($k\lesssim k_{0}$) the magnetoelastic energy to lowest order in the deviations of magnetization and lattice from equilibrium reads [28, 23, 24, 25, 26]

[TABLE]

where $n=1/a^{3}$. The strain tensor $X_{\alpha\beta}$ is defined in terms of the lattice displacements $X_{\alpha}$,

[TABLE]

with, for a cubic lattice [28],

[TABLE]

$B_{\alpha\beta}$ is caused by magnetic anisotropies and $B_{\alpha\beta}^{\prime}$ by the exchange interaction under lattice deformations. For YIG at room temperature [27, 33]

[TABLE]

We discuss the values for $B_{\parallel}^{\prime}$ and $B^{\prime}_{\perp}$ in Sec. III.3.

III.2 Anisotropy-mediated magnon-phonon interaction

The magnetoelastic anisotropy (19) is described by the Hamiltonian [28],

[TABLE]

with interaction vertices

[TABLE]

and

[TABLE]

The one magnon-two phonon process is of the same order in the total number of magnons and phonons as the two magnon-one phonon processes, but its effect on magnon transport is small, as shown in Appendix B.

III.3 Exchange-mediated magnon-phonon interaction

The exchange-mediated magnon-phonon interaction is obtained under the assumption that the exchange interaction $J_{ij}$ between two neighboring spins at lattice sites $\mathbf{r}_{i}$ and $\mathbf{r}_{j}$ depends only on their distance, which leads to the expansion to leading order in the small parameter $(|\mathbf{r}_{i}-\mathbf{r}_{j}|-a)$

[TABLE]

where $a$ is the equilibrium distance and $J^{\prime}=\partial J/\partial a$. With $\mathbf{r}_{i}=\mathbf{R}_{i}+\mathbf{X}_{\mathbf{R}_{i}},$ the Heisenberg Hamiltonian (1) is modulated by

[TABLE]

where $\mathbf{e}_{\alpha}$ is a unit vectors in the $\alpha$ direction. Expanding the displacements in terms of the phonon and magnon modes

[TABLE]

with interaction

[TABLE]

where the last line is the long-wavelength expansion. The magnon-phonon interaction

[TABLE]

conserves the magnon number, while (30) and (31) do not. Phonon numbers are not conserved in either case.

The value of $J^{\prime}$ for YIG is determined by the magnetic Grüneisen parameter [32, 33]

[TABLE]

where $V=Na^{3}$ is the volume of the magnet. The only assumption here is that the Curie temperature $T_{C}$ scales linearly with the exchange constant $J$ [43]. $\Gamma_{m}$ has been measured for YIG via the compressibility to be $\Gamma_{m}=-3.26$ [32], and via thermal expansion, $\Gamma_{m}=-3.13$ [33], so we set $\Gamma_{m}=-3.2$. For other materials the magnetic Grüneisen parameter is also of the order of unity and in many cases $\Gamma_{m}\approx-10/3$ [32, 44, 33]. A recent ab initio study of YIG finds $\Gamma_{m}=-3.1$ [45].

Comparing the continuum limit of Eq. (35) with the classical magnetoelastic energy (19)

[TABLE]

where for YIG $B_{\parallel}^{\prime}/a^{2}\approx 235\;\mathrm{meV}$. We disregard $B_{\perp}^{\prime}$ since it vanishes for nearest neighbor interactions by cubic lattice symmetry.

The coupling strength of the exchange-mediated magnon-phonon interaction can be estimated from the exchange energy $SJ^{\prime}a\approx E_{ex}=SJ$ [46, 31] following Akhiezer et al. [47, 48]. Our estimate of $SJ^{\prime}a=3\Gamma_{m}SJ$ is larger by $3\Gamma_{m}$, i.e. one order of magnitude. Since the scattering rate is proportional to the square of the interaction strength, our estimate of the scattering rate is a factor $100$ larger than previous ones. The assumption $J^{\prime}a\approx J$ is too small to be consistent with the experimental Grüneisen constant [32, 33]. Ref. [3] educatedly guessed $J^{\prime}a\approx 100J,$ which we now judge to be too large.

III.4 Interaction vertices

The magnon-phonon interactions in the Hamiltonian (27) are shown in Fig. 2 as Feynman diagrams. Fig. 2(a) illustrates magnon and phonon interconversion, which is responsible for the magnon-phonon hybridization and level splitting at the crossing of magnon and phonon dispersions [27, 28]. The divergence of this diagram at the magnon-phonon crossing points is avoided by either direct diagonalization of the magnon-phonon Hamiltonian [42] or by cutting-off the divergence by a lifetime parameter [31]. This process still generates enhanced magnon transport that is observable as magnon polaron anomalies in the spin Seebeck effect [22] or spin-wave excitation thresholds [49, 50], but these are strongly localized in phase space and disregarded in the following, where we focus on the magnon scattering rates to leading order in $1/S$ of the scattering processes in Fig. 2(b)-(d).

IV Magnon scattering rate

Here we derive the magnon reciprocal quasi-particle lifetime $\tau_{qp}^{-1}=\gamma$ as the imaginary part of the wave vector dependent self-energy, caused by acoustic phonon scattering [28],

[TABLE]

This quantity is in principle observable by inelastic neutron scattering. The total decay rate

[TABLE]

is the sum of the magnon number conserving decay rate $\gamma^{c}$ and the magnon number non-conserving decay rate $\gamma^{nc}$, which are related to the magnon-phonon scattering time $\tau_{mp}$ and the magnon-phonon dissipation time $\tau_{mr}$ by

[TABLE]

$\gamma^{\mathrm{other}}$ is caused by magnon-magnon and magnon disorder scattering, thereby beyond the scope of this work.

The self-energy to leading order in the $1/S$ expansion is of second order in the magnon-phonon interaction [28],

[TABLE]

where the magnon number conserving magnon-phonon scattering vertex $\Gamma_{\mathbf{k},\mathbf{k}^{\prime},\lambda}^{\bar{b}b}=\Gamma_{\mathbf{k},\mathbf{k}^{\prime},\lambda}^{ex}+\Gamma_{\mathbf{k},\mathbf{k}^{\prime},\lambda}^{an}$ and the Planck (Bose) distribution function $n_{B}(\varepsilon)=(e^{\beta\varepsilon}-1)^{-1}$ with inverse temperature $\beta=1/(k_{B}T)$. The Feynman diagrams representing the magnon number conserving and non-conserving contributions to the self-energy are shown in Fig. 3.

We write the decay rate in terms of four contributions

[TABLE]

where $out$ and $in$ denote the out-scattering and in-scattering parts. The contributions to the decay rate read [28]

[TABLE]

where the sum is over all momenta $\mathbf{q}$ in the Brillouin zone. Here the magnon/phonon annihilation rate is proportional to the Boson number $n_{B}$, while the creation rate scales with $1+n_{B}$. For example, in the out-scattering rate $\gamma_{out}^{c}(\mathbf{k})$ the incoming magnon with momentum $\mathbf{k}$ gets scattered into the state $\mathbf{k}-\mathbf{q}$ and a phonon is either absorbed with probability $\sim n_{B}$ or emitted with probability $\sim(1+n_{B})$. The out- and in-scattering rates are related by the detailed balance

[TABLE]

For high temperatures $k_{B}T\gg E_{\mathbf{k}}$, we may expand the Bose functions $n_{B}(E_{\mathbf{k}})\sim k_{B}T/E_{\mathbf{k}}$ and we find $\gamma_{in}\sim\gamma_{out}\sim T^{2}$ and $\gamma=\gamma_{out}-\gamma_{in}\sim T$. For low temperatures $k_{B}T\ll E_{\mathbf{k}}$, the out-scattering rate $\gamma_{out}\to\mathrm{const}.$ and the in-scattering rate $\gamma_{in}\sim e^{-\beta E_{\mathbf{k}}}\to 0$. The scattering processes (c) and (d) in Fig. 2 conserve energy and linear momentum, but not angular momentum. A loss of angular momentum after integration over all wave vectors corresponds to a mechanical torque on the total lattice that contributes to the Einstein-de Haas effect [51].

V Magnon transport lifetime

In this section we compare the transport lifetime $\tau_{t}$ and the magnon quasi-particle lifetime $\tau_{qp}$ that can be very different [52, 53, 54], but, to the best of our knowledge, has not yet been addressed for magnons. The magnon decay rate is proportional to the imaginary part of self energy, as shown in Eq. (41). On the other hand, the transport is governed by transport lifetime $\tau_{t}$ in the Boltzmann equation that agrees with $\tau_{qp}$ only in the relaxation time approximation. The stationary Boltzmann equation for the magnon distribution can be written as [3, 42]

[TABLE]

where $f_{\mathbf{k}}(\mathbf{r})$ is the magnon distribution function. The $in$ and $out$ contributions to the collision integral are related to the previously defined in- and out-scattering rates by

[TABLE]

where the equilibrium magnon distribution $n_{B}(E_{\mathbf{k}})$ is replaced by the non-equilibrium distribution function $f_{\mathbf{k}}$. The factor $(1+f_{\mathbf{k}})$ corresponds to the creation of a magnon with momentum $\mathbf{k}$ in the in-scattering process and the factor $f_{\mathbf{k}}$ to the annihilation in the out-scattering process. The phonons are assumed to remain at thermal equilibrium, so we disregard the phonon drift contribution that is expected in the presence of a phononic heat current.

Magnon transport is governed by three linear response functions, i.e. spin and heat conductivity and spin Seebeck coefficient [42]. These can be obtained from the expansion of the distribution function in terms of temperature and chemical potential gradients and correspond to two-particle Green functions with vertex corrections, that reflect the non-equilibrium in-scattering processes, captured by a transport lifetime $\tau_{t}$ that can be different from the quasi-particle (dephasing) lifetime $\tau_{qp}$ defined by the self-energy. We define the transport life time of a magnon with momentum $\mathbf{k}$ in terms of the collision integral

[TABLE]

with $f_{0,\mathbf{k}}=n_{B}(E_{\mathbf{k}})$ and we assume a thermalized quasi-equilibrium distribution function

[TABLE]

where $\mu$ is the magnon chemical potential. We linearize the function $f_{\mathbf{k}}$ in terms of small deviations $\delta f_{\mathbf{k}}$ from equilibrium $f_{0,\mathbf{k}}$,

[TABLE]

leading to [3]

[TABLE]

where the gradients of chemical potential $\boldsymbol{\nabla}\mu$ and temperature $\boldsymbol{\nabla}T$ drive the magnon current. In the relaxation time approximation we disregard the dependence of $\tau_{\mathbf{k},t}[f]$ on $\delta f$ and recover the quasi-particle lifetime $\tau_{\mathbf{k},t}\rightarrow\tau_{\mathbf{k},qp}$.

To first order in the phonon operators and second order in the magnon operators the collision integral for magnon number non-conserving processes,

[TABLE]

where the interaction vertex $\Gamma_{\mathbf{k},\mathbf{k^{\prime}},\lambda}^{bb}$ is given by Eq. (30) and $n_{\mathbf{\mathbf{q}\lambda}}=n_{B}(\varepsilon_{\mathbf{\mathbf{q}\lambda}})$. By using the expansion (56) in the collision integral that vanishes at equilibrium,

[TABLE]

we arrive at

[TABLE]

For the magnon number conserving process the derivation is similar and we find

[TABLE]

with interaction vertex $\Gamma_{\mathbf{k},\mathbf{k^{\prime}},\lambda}^{\bar{b}b}$ given by Eq. (38). Due to the $\delta f_{\mathbf{\mathbf{k}-\mathbf{q}}}/\delta f_{\mathbf{k}}$ term this is an integral equation. It can be solved iteratively to generate a geometric series referred to as vertex correction in diagrammatic theories. By simply disregarding the in-scattering with terms $\delta f_{\mathbf{\mathbf{k}-\mathbf{q}}}/\delta f_{\mathbf{k}}$ the transport lifetime reduces to the the quasi-particle lifetime of the self-energy. We leave the general solution of this integral equation for future work, but argue in Sec. VI.4 that the vertex corrections are not important in our regime of interest.

VI Numerical results

VI.1 Magnon decay rate

In the following we present and analyze our results for the magnon decay rates in YIG. We first consider the case of vanishing effective magnetic field ($B=0$) and discuss the magnetic field dependence in Sec. VI.3. Since our model is only valid in the long-wavelength ($k<8\text{\texttimes}10^{8}\;\mathrm{m}^{-1}$) and low-temperature ($T\lesssim 100\;\mathrm{K}$) regime, we focus first on $T=50\;\mathrm{K}$ and discuss the temperature dependence in Sec. VI.2.

In Fig. 4 we show the magnon number conserving decay rate $\gamma^{c}(\mathbf{k})$, which is on the displayed scale dominated by the exchange-mediated magnon-phonon interaction and is isotropic for long-wavelength magnons.

In Fig. 5 we compare the contribution from the exchange-mediated magnon-phonon interaction ($\gamma^{c}\sim k^{4}$) and from the anisotropy-mediated magnon-phonon interaction ($\gamma^{c}\sim k^{2}$). We observe a cross-over at $k\approx 4\times 10^{7}\;\mathrm{m}^{-1}$: for much smaller wave numbers, the exchange contribution can be disregarded and for larger wave numbers the exchange contribution becomes dominant.

The magnon number non-conserving decay rate $\gamma^{nc}$ in Fig. 6 is much smaller than the magnon-conserving one. This is consistent with the low magnetization damping of YIG, i.e. the magnetization is long-lived. We observe divergent peaks at the crossing points (shown in Fig. 1) with the exception of the (001) direction. These divergences occur when magnons and phonons are degenerate at $k=0.48\times 10^{9}\;\mathrm{m}^{-1}$ ($1.2\;\mathrm{meV}$) and $k=0.9\times 10^{9}\;\mathrm{m}^{-1}$ ($4.3\;\mathrm{meV}$), respectively, at which the Boltzmann formalism does not hold; a treatment in the magnon-polaron basis [42] or a broadening parameter [31] would get rid of the singular behavior. The divergences are also suppressed by arbitrarily small effective magnetic fields (see Sec. VI.3). There are no peaks along the (001) direction because in the (001) direction the vertex function $V_{\mathbf{q},\lambda}$ (see Eq. (33)) vanishes for $\mathbf{q}=(0,0,k_{z})$. For $k>\hbar c_{l}/(D(\sqrt{8}-2))=1.085\times 10^{9}\;\mathrm{m}^{-1}$ the decay rate $\gamma^{nc}$ vanishes because the decay process does not conserve energy ($\delta(E_{\mathbf{k}}+E_{\mathbf{q}-\mathbf{k}}-\varepsilon_{\mathbf{q}\lambda})=0$).

VI.2 Temperature dependence

Above we focused on $T=50\;\mathrm{K}$ and explained that we expect a linear temperature dependence of the magnon decay rates at high, but not low temperatures. Fig. 7 shows our results for the temperature dependence at $k_{x}=10^{8}\;\mathrm{m}^{-1}$. Deviations from the linear dependence at low temperatures occurs when quantum effects set in, i.e. the Rayleigh-Jeans distribution does not hold anymore,

[TABLE]

VI.3 Magnetic field dependence

The numerical results presented above are for a mono-domain magnet in the limit of small applied magnetic fields. A finite magnetic field $B$ along the magnetization direction induces an energy gap $g\mu_{B}B$ in the magnon dispersion, which shifts the positions of the magnon-phonon crossing points to longer wavelengths. The magnetic field suppresses the (unphysical) sharp peaks at the crossing points (see Fig. 8) that are caused by the divergence of the Planck distribution function for a vanishing spin wave gap.

In the magnon number conserving magnon-phonon interactions, the magnetic field dependence cancels in the delta function and enters only in the Bose function via $n_{B}$ (magnetic freeze-out). Fig. 9 shows that the magnetic field mainly affects magnons with energies $\lesssim 2g\mu_{B}B=0.23(B/\mathrm{T})\;\mathrm{meV}$.

As shown in Fig. 10 the magnon decay by phonons does not vanish for the $\mathbf{k}=0$ Kittel mode, but only in the presence of a spin wave gap $E_{0}=g\mu_{B}B$. Both magnon conserving and non-conserving scattering processes contribute. The divergent peaks at $B\approx 1.3\;\mathrm{T}$ and $B\approx 4.6\;\mathrm{T}$ in $\gamma^{nc}$ are caused by energy and momentum conservation in the two-magnon-one-phonon scattering process,

[TABLE]

when the gradient of the argument of the delta function vanishes,

[TABLE]

i.e., the two-magnon energy $E_{\mathbf{k}=0}+E_{\mathbf{q}}$ touches either the transverse or longitudinal phonon dispersion $\varepsilon_{\mathbf{q}\lambda}$. The total energy of the two magnons is equivalent to the energy of a single magnon with momentum $q$ but in a field $2B$, resulting in the divergence at fields that are half of those for the magnon-polaron observed in the spin Seebeck effect [42, 31]. The two-magnon touching condition can be satisfied in all directions of the phonon momentum $\mathbf{q}$, which therefore contributes to the magnon decay rate when integrating over the phonon momentum $\mathbf{q}$. For $\mathbf{k}\neq 0$ this two-magnon touching condition can only be fulfilled for phonons along a particular direction and the divergence is suppressed.

The magnon decay rate is related to the Gilbert damping $\alpha_{\mathbf{k}}$ as $\hbar\gamma_{\mathbf{k}}=2\alpha_{\mathbf{k}}E_{\mathbf{k}}$ [55]. We find that phonons contribute only weakly to the Gilbert damping, $\alpha_{0}^{nc}=\hbar\gamma_{0}^{nc}/(2E_{0})\sim 10^{-8}$ at $T=50\;\mathrm{K}$, which is much smaller than the total Gilbert damping $\alpha\sim 10^{-5}$ in YIG, but the peaks at $1.3\;\mathrm{T}$ and $4.6\;\mathrm{T}$ might be observable. The phonon contribution to the Gilbert damping scales linearly with temperature, so is twice as large at 100 K. At low temperatures ($T\lesssim 100\;\mathrm{K}$) Gilbert damping in YIG has been found to be caused by two-level systems [56] and impurity scattering [40], while for higher temperatures magnon-phonon [57] and magnon-magnon scattering involving optical magnons [34] have been proposed to explain the observed damping. Enhanced damping as a function of magnetic field at higher temperatures might reveal other van Hove singularities in the joint magnon-phonon density of states.

VI.4 Magnon transport lifetime

We do not attempt a full solution of the integral equations (60) and (LABEL:eq:tau_t^c) for the transport lifetime. However, we can still estimate its effect by the observation that the ansatz $\tau_{\mathbf{k},t}^{-1}\sim k^{n}$ can be an approximate solution of the Boltzmann equation with in-scattering.

Our results for the magnon number conserving interaction are shown in Fig. 11 (for $\boldsymbol{\nabla}T=0$ and finite $\boldsymbol{\nabla}\mu||\mathbf{e}_{x}$), where $\gamma_{t}=\tau_{t}^{-1}$. We consider the cases $n=0,2,4$, where $n=0$ or $\tau_{\mathbf{k},t}=\mathrm{const}.$ would be the solution for a short-range scattering potential. For very long wavelengths ($k\lesssim 4\times 10^{7}\;\mathrm{m}^{-1}$) the inverse quasi-particle lifetime $\tau_{\mathbf{k},qp}^{-1}\sim k^{2}$ and for shorter wavelengths $\tau_{\mathbf{k},qp}^{-1}\sim k^{4}$. $n=2$ is a self-consistent solution only for very small $k\lesssim 4\times 10^{7}\;\mathrm{m}^{-1}$, while $\tau_{\mathbf{k},qp}^{-1}\sim k^{4}$ is a good ansatz up to $k\lesssim 0.3\times 10^{9}\;\mathrm{m}^{-1}$. We see that the transport lifetime approximately equals the quasi-particle lifetime in the regime of the validity of the $n=4$ power law.

For the magnon number non-conserving processes in Fig. 12 the quasi-particle lifetime behaves as $\tau_{\mathbf{k},qp}^{-1}\sim k^{2}$. The ansatz $n=2$ turns out to be self-consistent and we see deviations of the transport lifetime from the quasi-particle lifetime for $k\gtrsim 5\times 10^{7}\;\mathrm{m}^{-1}$. The plot only shows our results for $k<1\times 10^{8}\;\mathrm{m}^{-1}$ because our assumption of an isotropic lifetime is not valid for higher momenta in this case.

We conclude that for YIG in the long-wavelength regime the magnon transport lifetime (due to magnon-phonon interactions) should be approximately the same as the quasi-particle lifetime, but deviations at shorter wavelengths require more attention.

VII Summary and conclusion

We calculated the decay rate of magnons in YIG induced by magnon-phonon interactions in the long-wavelength regime ($k\lesssim 1\times 10^{9}\;\mathrm{m}^{-1}$). Our model takes only the acoustic magnon and phonon branches into account and is therefore valid at low to intermediate temperatures ($T\lesssim 100\;\mathrm{K}$). The exchange-mediated magnon-phonon interaction has been recently identified as a crucial contribution to the overall magnon-phonon interaction in YIG at high temperatures [3, 29, 45]. We emphasize that its coupling strength can be derived from experimental values of the magnetic Grüneisen parameter $\Gamma_{m}=\partial\ln T_{C}/\partial\ln V$ [32, 33]. In previous works this interaction has been either disregarded [28], underestimated [29, 46], or overestimated [3].

In the ultra-long-wavelength regime the wave vector dependent magnon decay rate $\gamma(\mathbf{k})$ is determined by the anisotropy-mediated magnon-phonon interaction with $\gamma(\mathbf{k})\sim k^{2}$, while for shorter wavelengths $k\gtrsim 4\times 10^{7}\;\mathrm{m}^{-1}$ the exchange-mediated magnon-phonon interaction becomes dominant, which scales as $\gamma(\mathbf{k})\sim k^{4}$. The magnon number non-conserving processes are caused by spin-orbit interaction, i.e., the anisotropy-mediated magnon-phonon interaction, and are correspondingly weak.

In a finite magnetic field the average phonon scattering contribution, from the mechanism under study, to the Gilbert damping of the $k=0$ macrospin Kittel mode is about three orders of magnitude smaller than the best values for the Gilbert damping $\alpha\sim 10^{-5}$. However, we predict peaks at $1.3\;\mathrm{T}$ and $4.6\;\mathrm{T}$, that may be experimentally observable in high-quality samples.

The magnon transport lifetime, which is given by the balance between in- and out-scattering in the Boltzmann equation, is in the long-wavelength regime approximately the same as the quasi-particle lifetime. However, the magnon quasi-particle and transport lifetime differ more significantly at shorter wavelengths. A theory for magnon transport at room temperature should therefore include the “vertex corrections”.

A full theory of magnon transport at high temperature requires a method that takes the full dispersion relations of acoustic and optical phonons and magnons into account. This would also require a full microscopic description of the magnon-phonon interaction, since the magnetoelastic energy used here only holds in the continuum limit.

Acknowledgements.

N. V-S thanks F. Mendez for useful discussions. This work is part of the research program of the Stichting voor Fundamenteel Onderzoek der Materie (FOM), which is financially supported by the Nederlandse Organisatie voor Wetenschappelijk Onderzoek (NWO) as well as a Grant-in-Aid for Scientific Research on Innovative Area, ”Nano Spin Conversion Science” (Grant No. 26103006), CONICYT-PCHA/Doctorado Nacional/2014-21140141, Fondecyt Postdoctorado No. 3190264, and Fundamental Research Funds for the Central Universities.

Appendix A Long-wavelength approximation

The theory is designed for magnons with momentum $k<0.8\times 10^{9}\;\mathrm{m}^{-1}$ and phonons with momentum $q<2.5\times 10^{9}\;\mathrm{m}^{-1}$ (corresponding to phonon energies/frequencies \leq$$12\;\mathrm{meV}/$3\;\mathrm{THz}$), but relies on high-momentum cut-off parameters $k_{c}$ because of the assumption of quadratic/linear dispersion of magnon/phonons. We see in Fig. 13 that the scattering rates only weakly depend on $k_{c}$.

The dependence of the scattering rate on the phonon momentum cut-off $q_{c}$ is shown in Fig. 14. $q_{c}=3.15\times 10^{9}\;\mathrm{m}^{-1}$ corresponds to an integration over the whole Brillouin zone, approximated by a sphere. From these considerations we estimate that the long-wavelength approximation is reliable for $k\lesssim 8\times 10^{8}\;\mathrm{m}^{-1}$. Optical phonons (magnons) that are thermally excited for $T\apprge 100\;\mathrm{K}\,(300\;\mathrm{K})$ are not considered here.

Appendix B Second order magnetoelastic coupling

The magnetoelastic energy is usually expanded only to first order in the displacement fields. Second order terms can become important e.g. when the first order terms vanish. This is the case for one-magnon two-phonon scattering processes. The first order term

[TABLE]

only contributes when phonon and magnon momenta and energies cross, giving rise to magnon polaron modes [42]. In other areas of reciprocal space the second order term should therefore be considered. Eastman [58, 59] derived the second-order magnetoelastic energy and determined the corresponding coupling constants for YIG. In momentum space, the relevant contribution to the Hamiltonian is of the form

[TABLE]

where the interaction vertices are symmetrized,

[TABLE]

and obey

[TABLE]

The non-symmetrized vertex function is

[TABLE]

with

[TABLE]

and $x\leftrightarrow y$ denotes an exchange of $x$ and $y$. The relevant coupling constants in YIG are [59, 58]

[TABLE]

The magnon self-energy (see Fig. 15) reads

[TABLE]

with phonon propagator

[TABLE]

and leads to a magnon decay rate

[TABLE]

where

[TABLE]

The first term in curly brackets on the right-hand-side of Eq. (78) describes annihilation and creation of a phonon as a sum of out-scattering minus in-scattering contributions,

[TABLE]

while the second term can be understood in terms of out-scattering by the creation of two phonons and the in-scattering by annihilation of two phonons,

[TABLE]

For this one-magnon-two-phonon process the quasi-particle and the transport lifetimes are the same,

[TABLE]

since this process involves only a single magnon that is either annihilated or created. The collision integral is then independent of the magnon distribution of other magnons and the transport lifetime reduces to the quasi-particle lifetime.

The two-phonon contribution to the magnon scattering rate in YIG at $T=50\;\mathrm{K}$ and along (100) direction as shown in Fig. 16 is more than two orders of magnitude smaller than that from one-phonon processes and therefore disregarded in the main text. The numerical results depend strongly on the phonon momentum cutoff $q_{c}$, even in the long-wavelength regime, which implies that the magnons in this process dominantly interact with short-wavelength, thermally excited phonons. Indeed, the second order magnetoelastic interaction (69) is quadratic in the phonon momenta, which favors scattering with short-wavelength phonons. Our long-wavelength approximation therefore becomes questionable and the results may be not accurate at $T=50\;\mathrm{K}$, but this should not change the main conclusion that we can disregard these diagrams.

Our finding that the two-phonon contributions are so small can be understood in terms of the dimensionful prefactors of the decay rates (Eqs. (48-49) and (78)): The one-phonon decay rate is proportional to $\hbar/(ma^{2})\approx 7\times 10^{6}\;\mathrm{s^{-1}}$, while the two-phonon decay rate is proportional to $\hbar^{3}/(m^{2}a^{4}\varepsilon)\approx 33\;\mathrm{s}^{-1}$, where $\varepsilon\approx 1\;\mathrm{meV}$ is a typical phonon energy. The coupling constants for the magnon number non-conserving processes are $B_{\parallel,\perp}\sim 5\;\mathrm{meV}$ while the strongest two phonon coupling which enhances the two-phonon process by about a factor 100, but does not nearly compensate the prefactor. The two phonon process is therefore three orders of magnitudes smaller than the contribution of the one phonon process. The physical reason appears to be the large mass density of YIG, i.e. the heavy yttrium atoms.

Appendix C Numerical integration

The magnon decay rate is given be the weighted density of states

[TABLE]

that contain the Dirac delta function $\delta(\varepsilon)$ that can be eliminated to yield

[TABLE]

where the $\mathbf{q}_{i}$ are the zeros of $\varepsilon(\mathbf{q})$ and $A_{i}$ the surfaces inside the Brillouin zone with $\varepsilon(\mathbf{q})=\varepsilon(\mathbf{q}_{i})$. The calculation these integrals is a standard numerical problem in condensed matter physics.

For a spherical Brillouin zone of radius $q_{c}$ and spherical coordinates $(q,\theta,\phi)$,

[TABLE]

When $\varepsilon(q_{i},\theta,\phi)=0$

[TABLE]

where $\varepsilon^{\prime}=\partial\varepsilon/\partial q$ and

[TABLE]

which is particularly useful when the zeros of $\varepsilon(q,\theta,\phi)$ can be calculated analytically for linear and quadratic dispersion relations.

We can also evaluate the integral $I$ fully numerically by broadening the delta function [60] e.g. replacing it by a Gaussian [60],

[TABLE]

where $\sigma$ is the broadening parameter. An alternative is the Lorentzian (Cauchy-Lorentz distribution),

[TABLE]

which has fat tails that are helpful in finding the zeros of the delta function for an adaptive integration grid. Here we use the cubature package by Steven G. Johnson [61], which implements an adaptive multidimensional integration algorithm over hyperrectangular regions [62, 63].