Atomistic Simulation of Phonon and Magnon Thermal Transport across the Ferro-Paramagnetic Transition

TL;DR

This paper introduces a temperature-dependent simulation method combining atomic and spin dynamics to accurately model phonon and magnon thermal transport, including their interactions, across magnetic phase transitions, validated on body-centered cubic iron.

Contribution

The study develops a novel Green-Kubo equilibrium approach that captures phonon-magnon interactions and thermal conductivities across the ferro-para magnetic transition, filling a key gap in thermal transport modeling.

Findings

Successfully reproduces temperature-dependent thermal conductivities of BCC iron.

High-frequency phonon-magnon scattering rates are significantly larger than low-frequency ones.

The method aligns well with experimental measurements of non-electronic thermal conductivity.

Abstract

A temperature-dependent approach involving Green-Kubo equilibrium atomic and spin dynamics (GKEASD) is reported to assess phonon and magnon thermal transport processes accounting for phonon-magnon interactions. Using body-center cubic (BCC) iron as a case study, GKEASD successfully reproduces its characteristic temperature-dependent spiral and lattice thermal conductivities. The non-electronic thermal conductivity, i.e., the sum of phonon and magnon thermal conductivities, calculated using GKEASD for BCC Fe agrees well with experimental measurements. Spectral energy analysis reveals that high-frequency phonon-magnon scattering rates are one order of magnitude larger than those at low frequencies due to energy scattering conservation rules and high densities of states. Higher temperatures further accentuate this phenomenon. This new framework fills existing gaps in simulating thermal transport across the ferro- to para-magnetic transition. Future application of this methodology to phonon- and magnon-dominant insulators and semiconductors will enhance understanding of emerging thermoelectric, spin caloritronic and superconducting materials.

Equations

$$\frac{d{\vec{r}}_i}{dt}=\frac{{\vec{p}}_i}{m_i}$$

$$\frac{d{{{\vec{p}}}_{i}}}{dt}=\sum\limits_{i\neq j}^{N}{\left[-\frac{dV(|{{{\vec{r}}}_{ij}}|)}{d|{{{\vec{r}}}_{ij}}|}+\frac{dJ(|{{{\vec{r}}}_{ij}}|)}{d|{{{\vec{r}}}_{ij}}|}{{{\vec{s}}}_{i}}\cdot{{{\vec{s}}}_{j}}\right]{{{\vec{e}}}_{ij}}}\$$

$$\frac{d{\vec{s}}_i}{dt}={\vec{f}}_i\times {\vec{s}}_i$$

$$J\left(|r_{ij}|\right)=4\alpha {\left(\frac{|r_{ij}|}{\delta }\right)}^2\left(1-\gamma {\left(\frac{|r_{ij}|}{\delta }\right)}^2\right)e^{-{\left(\frac{|r_{ij}|}{\delta }\right)}^2}\Theta (R-|r_{ij}|)$$

$$-c(\vec{q},\nu )\frac{\partial \bar{n}(\vec{q},\nu )}{\partial T}=$$

$$\frac{1}{2}{\Gamma }_{(\vec{q},\nu )}^{({\vec{q}}^{\prime },{\nu }^{\prime }),({\vec{q}}^{\prime \prime },{\nu }^{\prime \prime })}(n_{(\vec{q},\nu )}^1-n_{({\vec{q}}^{\prime },{\nu }^{\prime })}^1-n_{({\vec{q}}^{\prime \prime },{\nu }^{\prime \prime })}^1)]$$

$$E=E_{kinetic}+E_{potential}+E_{spin}$$

$$=\sum _{i=1}^N\frac{|\vec{p}{|}^2}{2m_i}+\sum _{i,j=1}^{N}V({\vec{r}}_{ij})-\sum _{i,j,i\ne j}^{N}J(|{\vec{r}}_{ij}|){\vec{s}}_i\cdot {\vec{s}}_j$$

$$Q_{lattice}^{{}^{\prime \prime }x}=\frac{1}{V}\sum _i\frac{d(E_i^{kinetic}+E_i^{potential})r_i^x}{dt}$$

$$=\frac{1}{V}\left\{\sum _i{e}_i\cdot v_i^x+\frac{1}{2}\sum _{i=1}\sum _{j=1;j\ne i}\left({\vec{F}}_{ij}\cdot {\vec{v}}_i\right)\cdot r_{ij}^x\right\}$$

$$Q_{spin,ij}^{{}^{\prime \prime }}=$$

$$=$$

$$\frac{1}{S}\frac{dJ(|{\vec{r}}_{ij}|)}{dt}{\vec{s}}_i\cdot {\vec{s}}_j$$

$${\vec{r}}_{ij}$$

$$=[{\vec{r}}_i^0+{\vec{u}}_i(t)]-[{\vec{r}}_j^0+{\vec{u}}_j(t)]$$

$$={\vec{r}}_{ij}^0+{\vec{u}}_{ij}(t)$$

$$Q_{spin,ij}^{{}^{\prime \prime }}$$

$$\approx \frac{J(|{\vec{r}}_{ij}|)}{S}\left({\vec{s}}_i\cdot \frac{d{\vec{s}}_j}{dt}+{\vec{s}}_j\cdot \frac{d{\vec{s}}_i}{dt}\right)$$

$$Q_{spin}^{{}^{\prime \prime }}$$

$$=\frac{1}{2S}\sum _{i\in L}\sum _{i\in R}J(|{\vec{r}}_{ij}|)\left({\vec{s}}_i\cdot \frac{d{\vec{s}}_j}{dt}+{\vec{s}}_j\cdot \frac{d{\vec{s}}_i}{dt}\right)$$

$$=\frac{1}{S}\sum _{i\in L}\sum _{i\in R}J(|{\vec{r}}_{ij}|)\frac{d{\vec{s}}_j}{dt}\cdot {\vec{s}}_i$$

$$\frac{d{\vec{s}}_i}{dt}$$

$$Q_{spin}^{{}^{\prime \prime }x}=-\frac{1}{S\hslash }\sum _{i\in L}\sum _{i\in R}J(|{\vec{r}}_{ij}|)\left(\frac{\partial E_{mag}}{\partial {\vec{s}}_j}\times {\vec{s}}_j\right)\cdot {\vec{s}}_i$$

$$\kappa ={\kappa }_{phonon}+{\kappa }_{magnon}+{\kappa }_{cross}$$

$$=\frac{V}{k_b{T}^2}\int [\underset{phonon}{\underbrace{⟨Q_{lattice}^{{}^{\prime \prime }x}(t)\cdot Q_{lattice}^{{}^{\prime \prime }x}(0)⟩}}$$

$$+\underset{magnon}{\underbrace{⟨Q_{spin}^{{}^{\prime \prime }x}(t)\cdot Q_{spin}^{{}^{\prime \prime }x}(0)⟩}}+\underset{cross}{\underbrace{2⟨Q_{spin}^{{}^{\prime \prime }x}(t)\cdot Q_{lattice}^{{}^{\prime \prime }x}(0)⟩}}]dt$$

$$\omega (\vec{q},\nu )\vec{e}(\vec{q},\nu )=\omega (\vec{q},\nu )\text{D}(\vec{q})$$

$$\omega (\vec{k})\vec{e}(\vec{k})=\text{D}(\vec{k})\vec{e}(\vec{k})$$

$$\text{D}(\vec{q})=\frac{1}{m}\sum _{l}K(0,l)\exp (i\vec{k}[\vec{r}(l)-\vec{r}(0)])$$

$$\text{D}(\vec{k})=\frac{1}{\hslash }\sum _{l}J(0,l)\cdot \left[1-\exp (i\vec{k}[\vec{r}(l)-\vec{r}(0)])\right]$$

$${\kappa }_{phormag}=\sum _{(\vec{q},\nu )or(\vec{k})}{C}_V{v}_g^2\tau$$

$${\left|\Phi \right|}^2=\frac{I_p}{{\left[(\omega -{\omega }_p)/\Delta \right]}^2+1}$$

$$\Phi (\vec{q},\nu )\sim$$

$${\vec{e}}^{\ast }(j,\vec{q},\nu )\vec{v}(jl,t)dt$$

$$\Phi (\vec{k},\mu )\sim$$

$${\vec{e}}^{\ast }(j,\vec{k},\mu )\frac{d\vec{s}(jl,t)}{dt}dt$$

$$P_{order}=\sqrt{\sum _{\alpha =x,y,z}{\left(\sum _i{s}_{i,\alpha }^T\right)}^2/\sum _{\alpha =x,y,z}{\left(\sum _i{s}_{i,\alpha }^{T=0}\right)}^2}$$

Full text

I I. INTRODUCTION

A better understanding of heat transfer considering interactions among different temperature-induced excitations in crystals, e.g., phonons, electrons and spins, is of great importance in many disciplines, including thermoelectric Uchida2010 ; Flipse2012 ; Uchida2011 , spin caloritronic Bauer2012 and superconducting Chen1988 materials. Unlike phonons and electrons, whose thermal transport properties have been well studied, the heat transport behavior of magnons - collective excitations of magnetic spins - is poorly understood at the fundamental level, and little is known about the influence of phonon-magnon and magnon-magnon scattering on heat transfer. For example, the experimental lattice thermal conductivities of magnetic materials such as body-centered cubic (BCC) iron Fulkerson1966 ; Backlund1961 , face-centered cubic (FCC) nickel Powell1965 , YMnO${}_{\text{3}}$, LuMnO${}_{\text{3}}$ and ScMnO Sharma2004 , as well as CrN Tomes2011 ; Jankovsky2014 , show significantly different temperature dependences near the Curie temperature ($T_{\text{c}}$), compared to the typical $1/T$ relation at high temperatures when only phonons are considered Holland1963 ; Callaway1959 ; Zhou2018 . The magnetic configuration changes from a ferromagnetic state at low temperature, in which spins are aligned to a paramagnetic state above the Curie temperature, in which the spin configuration is disordered Kormann2014 .

To model the thermal transport properties of magnetic materials, temperature-dependent lattice and magnetic excitations must be taken into consideration when calculating the inputs, e.g., heat current or force constants, for thermal conductivity. To date, theoretical studies have sought to develop explanations for mutual interactions between phonons and magnons Kormann2014 ; Stockem2018 ; Sabiryanov1999 ; Liao2014 ; Liu2017 . However, most such studies provide only a partial treatment based on thermodynamic properties Kormann2014 ; Dick2012 , average phonon-magnon relaxation times Stockem2018 , and phonon-magnon temperature gradients Liao2014 . A robust framework to compute thermal transport properties such as thermal conductivity and modal scattering rates in magnetic materials is still lacking.

In this paper, a temperature-dependent method - Green-Kubo equilibrium atomic and spin dynamics (GKEASD) - based on linear response theory Kubo1957 and spin-lattice dynamics Tranchida2018 ; Ma2016 is reported to calculate the thermal transport properties of phonons and magnons in magnetic materials. In GKEASD, phonon-phonon scattering, phonon-magnon scattering and magnon-magnon scattering are inherently included. As a proof of concept, our study focuses on a model system consisting of a simple transition-metal ferromagnet, BCC iron, with a Curie temperature of 1043 K. Using this methodology, we successfully reproduce the characteristic temperature-dependent non-electronic thermal conductivity observed in experiments for magnetic materials Fulkerson1966 ; Backlund1961 ; Powell1965 ; Sharma2004 ; Tomes2011 ; Quintela2009 . Modal level phonon-phonon, magnon-magnon and phonon-magnon scattering rates are then quantified using spectral energy density analysis. The agreement between theoretical predictions and experimental measurements establishes the reliability of the methodology.

II II. Computational methods

II.1 A. Spin Lattice Dynamics

Here, a symplectic and scalable algorithm for spin lattice dynamics that was recently developed Tranchida2018 ; Ma2016 and embedded in the Large-scale Atomic/Molecular Massively Parallel Simulator (LAMMPS) Plimpton1995 is applied to describe the atomic spins in magnetic crystals from the ferromagnetic limit to the paramagnetic limit (Figure  1). Compared to classical molecular dynamics, this new algorithm augments the phase space by adding a classical spin vector $\vec{s}$ to each magnetic atom $i$, in addition to its position $\vec{r}$ and momentum $\vec{p}$. The motion equations of atoms and spins can be written as

[TABLE]

[TABLE]

[TABLE]

where ${{\vec{e}}_{ij}}$ is the unit vector along ${{\vec{r}}_{ij}}$, $\vec{f_{i}}$ is the analog of a spin force applied on the spin. $J(|{{{\vec{r}}}_{ij}}|)$ is the magnetic coupling exchange constant, which originates from two main contributions: (i) direct ferromagnetic exchange between the orbitals localized on centers of ions $i$ and $j$, and (ii) spin and charge polarizations effects carried by nonmagnetic orbitals. Spin dynamics simulations without lattice vibrations, i.e., fixed atomic position, are implemented in a $NVT$ ensemble (constant number of particles, volume and temperature) for 200 ps to reach the target temperature, and then a $NVE$ ensemble (constant number of particles, volume and energy) is run for 200 ps to compute the heat flux associated with the spins. In the $NVT$ simulations, the thermal fluctuations in the magnetic system are described following the Langevin approach Tranchida2018 ; Antropov1995 ; Carcia1998 .

For a system with lattice vibrations, a $NVT$ simulation for 200 ps is first run to reach the target temperature. Then, a $NVE$ ensemble for 200 ps is used to obtain the heat flux fluctuation of the magnons. A time step of 0.2 fs is used in all simulations. The spin potential used in our study is fitted using first-principles calculations (details can be found in Tranchida2018 ). The magnon dispersion generated from this potential fits experiments well (see details below). Mechanical interactions among atoms are computed using the embedded-atom method potential Chamati2006 , and the interactions between spins is modeled with the Bethe-Slater curve Yosida1996

[TABLE]

where the fitting parameters $\alpha$, $\gamma$ and $\delta$ are 25.498 meV, 0.281 and 1.999 Å, respectively based on the $J$ values for BCC iron Pajda2001 . $\Theta(R-{{r}_{ij}})$ is the Heaviside step function, and $R$ is the cutoff radius, which is 4 Å in our study.

II.2 B. First Principles Calculations

In this paper, all first-principles simulations are implemented via Vienna Ab-initio Simulation Package (VASP) based on density functional theory. The pseudopotential with a generalized gradient approximation parameterized by Perdew-Burke-Ernzerhof theory Perdew1986 for the exchange-correlation functional is used to depict the system. Periodic boundary conditions are applied in the three directions. A plane wave basis with a cutoff energy of 520 eV is used in all simulations. The Monkhorst-Pack scheme is used to generate 2$\times$2$\times$2 k-point mesh for the primitive cell and 6$\times$6$\times$6 supercell. Before performing electrostatic potential or interatomic force constant calculations, the atomic structure and cell size are fully relaxed until the energy difference and the Hellman-Feynman force converge to within 1$\times 10^{-6}$ eV and 1$\times 10^{-5}$ eV/Å, respectively. The 6$\times$6$\times$6 supercell with an energy criterion of 1$\times 10^{-6}$ eV in the self consistent calculation is used to extract 2nd and 3rd order force constants.

II.3 C. Green-Kubo Equilibrium Atomic and Spin Dynamic Simulations

During molecular dynamics simulations for lattice thermal conductivity calculations, we first run the $NVT$ ensemble for 200 ps with a time step of 0.2 fs to allow the lattice and spin fields attain the target temperature. Then a canonical ensemble is run for the next 400 ps to generate heat current, which is the input for the calculation of lattice thermal conductivity using Green-Kubo equilibrium atomic and spin dynamics. For each case, 30 independent runs are performed to obtain a stable average thermal conductivity. The correlation time considered in our simulations is 16 ps, which is long enough to obtain a converged and steady thermal conductivity (Figure  2a). An 8$\times$8$\times$8 unit cell box is used in all equilibrium molecular dynamics simulations, for which size effects may be ignored (Figure  2b). Lattice expansion due to temperature and spin is considered in our results (Figure  2c).

For the magnon thermal conductivity computations, 30 independent runs were used to obtain a converged thermal conductivity. A domain of 8$\times$8$\times$8 unit cells was used, and the correlation time was 16 ps, which is long enough to obtain a converged and size-independent magnon thermal conductivity (Figure  2b).

II.4 D. Validation of Classical Mechanical Potential

To assess the accuracy of the potential, we compare the phonon dispersion (Figure  3a) as well as lattice thermal conductivity (Figure  3b) from molecular dynamics and $ab$ $initio$ calculations. Both the phonon dispersion and the lattice thermal conductivity computed using a classical embedded-atom method potential agree with those calculated using a first-principles approach somehow. The $ab$ $initio$ thermal conductivity is calculated via the Boltzmann transport equation which can be solved by expanding the scattering term into its first-order perturbation $n^{1}$. Without considering impurities and boundaries, the linearized Boltzmann transport equation may be recast as Li2014

[TABLE]

where $(\vec{q},\nu)$ is the phonon mode with wave vector $\vec{q}$ and branch $\nu$. $c$, $\bar{n}$ and $T$ are the specific heat capacity, equilibrium phonon population and system temperature, respectively. $\Gamma$ is the scattering rate at equilibrium of a process where two phonons combine to generate a third phonon or when one phonon splits into two phonons. The scattering rate matrix can be obtained using Fermi’s Golden Rule. The lattice thermal conductivity may be found by computing the heat flux from the computed phonon distribution and using Fourier’s law. The wave vector mesh density is varied from 14$\times$14$\times$14 to 17$\times$17$\times$17, and the neighbor cutoffs range from 1st to 5th. The latter are used to compute the 3rd order force constant which is key parameter to obtain the lattice thermal conductivity in Boltzmann transport equation. Based on these variations, a wave vector mesh of 17$\times$17$\times$17 with a 5th neighbour cutoff was chosen in our all Boltzmann transport equation calculations to produce a converged lattice thermal conductivity (Figure  3c).

III III. Theory

III.1 A. Thermal Conductivity based on Heat Flux Fluctuation

The energy $E$ of magnetic crystals must account for terms coupling the magnetic spins ($E_{spin}$) to the lattice ($E_{kinetic}$ and $E_{potential}$) through the following expression Tranchida2018 :

[TABLE]

in which $\vec{r}$, $\vec{p}$ and $\vec{s}$ are the position, momentum and spin vectors of the atoms, respectively; the negative sign on the last term indicates that the ground state energy of the system considering spin becomes lower. Phonons and magnons are coupled in Eq. (6) via atomistic positions, i.e., $\vec{r}_{ij}$. The heat flux due to lattice vibrations takes the form Hardy1963 :

[TABLE]

where ${{e}_{i}}$ and ${{\vec{v}}_{i}}$ are the energy and the velocity of atom $i$, respectively. ${{\vec{F}}_{ij}}$ represents the force between two atoms, $V$ and $S$ are the volume and the cross section area of the system, respectively. and the heat flux associated with spin is expressed as

[TABLE]

where $|{{\vec{r}}_{ij}}|$ is the distance between atoms $i$ and $j$:

[TABLE]

in which, $\vec{r}^{0}$ and $\vec{u}$ are the equilibrium position and displacement of atoms, respectively. Since $|\vec{r}_{ij}^{0}|\gg|{{\vec{u}}_{ij}}|$ (Figure  4), we assume that ${{\vec{r}}_{ij}}$ is time independent. Therefore, Eq. (9) can be rewritten in the form

[TABLE]

At the same time, if we assume that atoms $i$ and $j$ are separated by an imaginary interface, the heat flux across the interface can be expressed as

[TABLE]

where $L$ and $R$ indicate left and right sides of the imaginary interface, respectively. The factor 1/2 addresses the ergodicity of the system. In addition, for spin lattice dynamics systems,

[TABLE]

in which ${{\vec{f}}_{i}}$ is the analog of a force applied on the spin. Finally, the heat flux across the imaginary interface can be written as

[TABLE]

Based on linear response theory Kubo1957 , Eq. (7) and Eq. (13), the thermal conductivity of a magnetic system can be divided into contributions from lattice vibrations, ${{\kappa}_{phonon}}$, spin-related fluctuations, ${{\kappa}_{magnon}}$, and a term resulting from lattice-spin interactions, ${{\kappa}_{cross}}$:

[TABLE]

where ${{k}_{b}}$ is the Boltzmann constant, and $T$ is the temperature of the system. We note that ${{\kappa}_{cross}}$ is not the result of phonon-magnon scattering, but of the cross-correlation between the phonon and magnon heat fluxes. Physically, the cross term represents the interaction between the heat carried by the phonons and the heat carried by the magnons that can alter pure phonon and spin heat flow. The effect of magnon scattering on phonon transport is elucidated by calculating the phonon thermal conductivity separately using ($E={{E}_{kinetic}}+{{E}_{potential}}+{{E}_{spin}}$) and ($E={{E}_{kinetic}}+{{E}_{potential}}$) and comparing the two values. Similarly, the influence of phonon scattering on magnon thermal conductivity is assessed by computing the magnon thermal conductivity using ($E={{E}_{kinetic}}+{{E}_{potential}}+{{E}_{spin}}$) and ($E={{E}_{spin}}$) and comparing the two values.

III.2 B. Spectral Energy Density

The atomistic velocity and spin change frequency used for the spectral energy density calculations were output every 4 fs over a total sampling time of 80 ps. All the reduced results were averaged across 3 runs and using two time intervals of 40 ps each. A Gauss window of 0.1 THz was used to filter noise in the original data.

The eigenvectors and frequencies of the phonon mode $(\vec{q},\nu)$ and magnon mode $(\vec{k})$ used in our simulations are computed by solving the phonon dynamical matrix equation

[TABLE]

and the magnon dynamical matrix equation

[TABLE]

where $\text{D}(\vec{q})$ and $\text{D}(\vec{k})$ are the lattice and spin dynamical matrices

[TABLE]

and

[TABLE]

respectively. Here, ${K}$ is the force constant matrix and $J(0,\ l)$ is magnetic exchange constant matrix. We note that Eq. (18) is derived assuming that the material is ferromagnetic, in which all spins have the same magnitude in one direction (only under this assumption can the spin motion equation be reduced to Eq. (16) Kittel1987 ). The volumetric heat capacity $C_{V}$ and group velocity $v_{g}$ of phonons and magnons are shown in Figure 5. Together with relaxation time $\tau$, the phonon or magnon thermal conductivity ${\kappa}_{ph\ or\ mag}$ can be obtained via

[TABLE]

where $\tau=1/2\Delta$. The linewidth $\Delta$, which is half the scattering rate $\Gamma$, can be calculated via spectral energy density analysis as:

[TABLE]

where ${{I}_{p}}$ and ${{\omega}_{p}}$ are the magnitude and frequency at the peak center, respectively. $\Phi$ is the spectral energy, which takes the following form for phonons Thomas2010 ; Zhou2015 :

[TABLE]

For magnons, the expression becomes Wu2018 :

[TABLE]

where $\vec{e}$ is the mode eigenvector of a phonon or magnon, and ${{t}_{0}}$ is the integration limit.

IV IV. Results

IV.1 A. Non-electronic Thermal Conductivity

From Eqs. (7) and (13), the net heat flux of lattice vibrations ($Q_{lattice}^{{}^{\prime\prime}x}$) and spin fluctuations ($Q_{spin}^{{}^{\prime\prime}x}$) in EMD simulations should be zero, which is validated by our numerical results (Figure  6a and 6****b). Furthermore, following linear response theory, the thermal conductivity of the two heat carriers should converge with increasing correlation time [Eq. (14)], and this is also reflected in our simulations (Figure  2b).

To evaluate thermal conductivity behavior at different temperatures, the order of the spin configuration in our system is calculated via

[TABLE]

where ${{P}_{order}}=1$ in the ferromagnetic limit and ${{P}_{order}}=0$ in the paramagnetic limit. Small values of ${{P}_{order}}$ indicate that spins in the system are predominantly disordered, which leads to broken periodicity.

By considering the contributions of both phonons and magnons, the predicted total thermal conductivity is in broad agreement with experimental measurements Fulkerson1966 ; Backlund1961 over the temperature range from 300 K to 1200 K (Figure  7a) (the electrical thermal conductivity was subtracted from the experimental measurements). We note that the two experimental results Fulkerson1966 ; Backlund1961 in Figure  7a are themselves quite different from each other, and the disagreement between them is larger than that caused by phonon-magnon scattering in the simulations. Thus, they cannot definitively establish the accuracy of our phonon-magnon scattering models. A number of reasons may account for these differences: (i) loss of heat from the specimen due to conduction through leads, and (ii) impurities in the sample, which in Ref. Fulkerson1966 may be as high as 1%.

An important feature in Figure  7a is the sharp fall in thermal conductivity reported in the experimental data of Ref. Backlund1961 near the Curie temperature (around 1043 K) due to the ferromagnetic-to-paramagnetic transition. Previously published atomic spin dynamics or spin lattice dynamics simulations have been known to smooth the ferromagnetic-to-paramagnetic phase transition near the Curie temperature Bergqvist2018 ; Szilva2013 (see ${{P}_{order}}$ with the original $J$ in Figure  7a). One explanation is that the simulation of atomic spins is performed within a classical framework and ignores quantum effects Turney2009 ; Szilva2013 . However, as discussed below and well established in Ref. Turney2009 , quantum effects can be ignored for phonon and magnon thermal transport properties for BCC iron in the temperature range considered here. Another explanation for the transition is that the value of the magnetic exchange parameter $J$ varies as a function of temperature Bergqvist2018 ; Ruban2016 .

Our EMD simulations (black and blue symbols) in Figure  7a use a single, temperature-invariant value of $J$ and treat the spins as classical, and therefore the results do not exhibit the experimentally observed sharp fall in thermal conductivity. To explore this issue further, an additional set of computations was performed (Figure  7b) wherein the value of $J$ varies with temperature. Here, we adjust the $J$ to make sure the spin in the system in the paramagnetic state is fully disordered because $J$ is only parameter that determines the spin configuration in the system. We assume that $J$ is constant when $T<{{T}_{c}}$, while $J$ is one tenth of the low temperature value when $T>{{T}_{c}}$. Thus, ${{P}_{order}}$ becomes zero when $T>{{T}_{c}}$ (Figure  7b), indicating that the spin configuration in the system is fully disordered (i.e., the paramagnetic limit). Consequently, the non-electronic thermal conductivity drops sharply around the Curie temperature due to the decrease in magnon thermal conductivity (Figure  7c), consistent with the experimental observations Fulkerson1966 ; Backlund1961 .

At the same time, phonon-magnon scattering does not exhibit a strong effect on thermal transport in BCC iron because the values of both the phonon and the magnon thermal conductivities are relatively small, i.e., phonon-phonon scattering and magnon-magnon scattering are the dominant scattering mechanisms in such materials.

IV.2 B. Phonon-magnon Scattering Process

To further understand the temperature behavior of thermal conductivity, mode-level phonon and magnon scattering rates have been calculated at three temperatures, 300 K, 700 K and 1100 K, by spectral energy density analysis. Using Matthiessen’s rule, phonon-magnon ($\Gamma_{phonon-magnon}$, in units of $1/ps$) and magnon-phonon ($\Gamma_{magnon-phonon}$) scattering rates can be calculated via:

[TABLE]

[TABLE]

where the superscripts $magnon$ and $no\leavevmode\nobreak\ magnon$ indicate that the lattice vibrates with and without spins in the system, respectively, and vice versa for $phonon$ and $no\leavevmode\nobreak\ phonon$.

First, the phonon (magnon) scattering rate of the same system with (without) considering the effects of magnons (phonons) was calculated based on spectral energy density. For phonons, our results indicate that (Figure  8a) the scattering rate in the low-frequency region changes little with the introduction of spin. Low-frequency phonons are known to be the main contributors to thermal conductivity at low temperatures whereas high-frequency phonons (i.e., short mean free path phonons) are important for lattice thermal conductivity at high temperatures Cuffe2015 . We explain why the lattice thermal conductivity considering spin is similar to that without it at low temperatures (see Figure  9a). For magnons (Figure  8b), the scattering rates for the systems with and without lattice vibrations are similar in this regime (below 50 THz), which indicates that the influence of phonon-magnon scattering on magnon transport in the low-frequency region is not important. On the other hand, in the high-frequency regime (above 50 THz), magnons can be strongly scattered by phonons.

From the phonon-magnon scattering results (Figure  9a), we observe a general tendency that high-frequency phonons ($\omega$ larger than 4.3 THz) are scattered by magnons more strongly than low-frequency phonons because the magnon energy is much larger than the phonon energy (Figure  9c). For phonon-dominant magnetic materials in which the thermal conductivity is mainly due to phonons, magnons may scatter phonons Chernyshev2015 ; Stamokostas2017 through phonon emission or absorption, $\hbar\omega({\vec{k}}^{\prime},\ {\mu}^{\prime})=\hbar\omega(\vec{q},\ \nu)+\hbar\omega(\vec{k},\ \mu)$. Referring to Figure  9c, and considering a magnon of energy $\hbar{{\omega}_{A}}$ scattered by a magnon of energy $\hbar{{\omega}_{B}}$, a high-frequency phonon would have a greater probability of involvement in a phonon-magnon scattering process than a low-frequency phonon. For instance, only two channels exist for the magnon $\hbar{{\omega}_{A}}$ to be scattered to the magnon $\hbar{{\omega}_{B}}$ when the phonon energy is $\hbar\omega=12.41$ meV ($\omega=3$ THz), whereas the number of channels for the magnon $\hbar{{\omega}_{A}}$ to be scattered to the magnon $\hbar{{\omega}_{C}}$ is four when the phonon energy is 37.6 meV ($\omega=9.1$ THz, the highest phonon frequency). From the magnon dispersion curve, the frequency changes for the small (P1-P2) and large peaks (P2-P3) are 4.3 and 3.6 THz, respectively; therefore, phonons (magnons) above 4.3 (68) THz have a much greater chance to be involved in phonon-magnon scattering processes.

Another reason that high-frequency phonons are preferentially involved in phonon-magnon scattering processes is the high density of states of phonons in the high-frequency region (Figure  9c). First-principles calculations Kormann2014 and experimental measurements Mauger2014 also confirm that high-frequency phonons are more strongly scattered by magnons. For materials in which the magnon thermal conductivity is dominant, magnons are scattered by phonons Chernyshev2015 ; Stamokostas2017 via phonon absorption or emission processes $\hbar\omega({\vec{k}}^{\prime},\ {\mu}^{\prime})\pm\hbar\omega(\vec{q},\ \nu)=\hbar\omega(\vec{k},\ \mu)$. The magnon-phonon scattering rate of high-frequency magnons is somewhat higher than that of low-frequency magnons (Figure  9b) because high-frequency magnons have a larger density of states (Figure  9c). The phenomena discussed above become more apparent at increased temperatures, which strengthen the scattering among the heat carriers.

IV.3 C. Accumulated Thermal Conductivity

To facilitate the analysis of thermal conductivity contributions, the thermal conductivity accumulation function has been computed with respect to the mean free path $\Lambda$ using

[TABLE]

where ${{c}_{v}}$, ${{v}_{g}}$ (Figure  5) and $\tau$ are the volumetric heat capacity, group velocity, and relaxation times of phonons and magnons, respectively. For phonons (Figure  10a), at $T$ = 300 K, the accumulated thermal conductivities with and without spin are similar, indicating that magnons do not have a strong effect on phonon thermal conductivity near room temperature. This is due to the fact that long mean free path phonons are the main contributors to thermal conductivity, and these phonons do not easily scatter with magnons, as discussed above. However, at 700 K, the phonon mean free path falls to 0.3 - 4 nm, which is much smaller than that at 300 K due to strong phonon-phonon scattering (black triangle symbols in Figure  10a). Furthermore, high-frequency phonons with short mean free paths (0.4 - 0.8 nm) are more easily scattered by magnons. Consequently, the total lattice thermal conductivity of magnetic BCC iron becomes slightly smaller than that of non-magnetic BCC iron (Figure  8). Finally, at 1100 K phonons are scattered strongly by both phonons and magnons because high temperatures increase the vibration magnitude of lattice and spin excitations.

For magnons (Figure  10b) at room temperature, only the long-mean-free-path, or equivalently, low-frequency magnons are scattered by phonons because only these magnons transport significant energy, i.e., the heat capacity of high frequency magnons is very small (Figure  5a). When the system temperature increases to 700 K and then to 1100 K, the mean free path of magnons becomes much shorter than that at 300K, and the short mean free path (high frequency) magnons are scattered strongly by phonons.

V V. Discussion

We also note that the magnon group velocities are calculated here from the magnon dispersion without considering temperature effects, i.e., assuming a ferromagnetic system. Such an assumption may introduce inaccuracies in the magnon mean free path and modal thermal conductivity computations. At 300 K, the total magnon thermal conductivity calculated using the Boltzmann transport equation, i.e., using $\kappa=\sum{{{c}_{v}}v_{g}^{2}\tau}$ with the values of relaxation time computed using spectral energy density, is 10.1 W/mK. This value is lower than the Boltzmann transport equation results of Wu $et.$ $al.$ (15.2 W/mK ) due to their overestimation of the magnon dispersion Wu2018 , whereas it is 7.8 W/mK computed by GKEASD. The closer correspondence between Boltzmann transport equation and GKEASD at room temperature indicates that the system may be treated as ferromagnetic at 300 K. However, the total magnon thermal conductivities calculated by the Boltzmann transport equation (GKEASD) are 11.1 (1.9) and 56.9 (1.2) W/mK at 700 K and 1100 K, respectively. The large differences between spectral energy density-Boltzmann transport equation and GKEASD indicate that the magnon group velocity is overestimated, and the calculation should consider spin disorder at high temperatures.

Before concluding our study, we briefly highlight two limitations that our methodology inherited from the spin lattice dynamics formalism. In our calculations, the exchange integral $J_{ij}$ is assumed to remain constant, and its dependence on temperature is neglected. Recent studies have investigated this dependence and proposed methodologies to account for it Szilva2013 ; Ruban2016 . Encapsulating them within our framework could improve the accuracy of predictions. Furthermore, the simulation of classical spins (instead of quantum spins) is known to make the ferromagnetic to paramagnetic phase transition smoother than in the experimental observations. It has been shown in Bergqvist2018 ; Szilva2013 that implementing quantum baths and statistics can reproduce more accurately the sharp transitin at $T_{c}$. However, as the GKESAD framework is based on equilibrium molecular dynamics, the associated simulations do not involve a connection to a random bath. Overall, an empirical parametrization of the exchange integral $J_{ij}$ could account for both its temperature dependence and the sharp transition occuring at $T_{c}$.

VI VI. Conclusions

In conclusion, we have developed a temperature-dependent method, the Green-Kubo equilibrium atomic and spin dynamics method, to calculate coupled phonon and magnon transport in magnetic materials. Reasonably good agreement is obtained between our simulation results and experimental measurements in computing dispersion curves and temperature-dependent thermal conductivity, and these results suggest that the approach captures the overall heat transfer behavior of phonons and magnons in magnetic crystals. However, the coupling between is governed by the magnetic exchange correlation constant, whose origins and fidelity deserve further scrutiny, as do other magnetic phenomena not considered here, such as anisotropy Finizio2014 . Analysis of scattering processes between phonons and magnons indicates that high-frequency phonon scattering rates due to phonon-magnon scattering are much larger than those at low frequencies because of energy-conserving rules for scattering and the high density of states. The application of this new methodology will yield deeper insights into the thermal transport properties of other phonon- or magnon-dominant materials.

VII Acknowledgments

This work was supported by the National Science Foundation (NSF) (Project number: NSF1758004) and used computational and storage services associated with the Hoffman2 Shared Cluster provided by UCLA Institute for Digital Research and Education’s Research Technology Group. This work used the Extreme Science and Engineering Discovery Environment (XSEDE), which is supported by National Science Foundation grant number DMR180057. Y. Z. gratefully acknowledges Dr. Zheyong Fan (Aalto University) for valuable comments and proof reading the manuscript. J. T. acknowledges that Sandia National Laboratories is a multi-mission laboratory managed and operated by National Technology and Engineering Solutions of Sandia, LLC, a wholly owned subsidiary of Honeywell International Inc., for the U.S. Department of Energy’s National Nuclear Security Administration under contract DE-NA0003525. This paper describes objective technical results and analysis. Any subjective views or opinions that might be expressed in the paper do not necessarily represent the views of the U.S. Department of Energy or the United States Government.