Thermal boundary conductance across Co/Cu interfaces with spin–lattice interactions

The origin of spin caloritronics is the so-called spin Seebeck effect (SSE) related to the generation of spin voltage driven by a temperature gradient. Both conduction electrons and spin waves (magnons) carry spin currents, but the range of the former is only hundreds of nanometers while the latter can persist for millimeters. Experiments have observed SSE in magnetic insulators where conduction electrons are absent, and this phenomenon highlights the role of spin waves in producing SSE.^1,2 SSE has also been detected in a $Ni81Fe19$/Pt wire on an electrically and magnetically isolated sapphire substrate³ to demonstrate that SSE originates from a non-equilibrium distribution of magnons in the ferromagnet that interact with phonons in the substrate. The phonon-driven redistribution of magnons has also been studied theoretically with a phenomenological phonon–magnon drag model.⁴ 

In the area of heat transport, magnon–phonon interactions offer potential applications in engineering thermal devices with magnetic materials. In our previous work,⁵ the effect of magnon–magnon and magnon–phonon interactions on the thermal conductivity of BCC iron was modeled by combining molecular dynamics and spin dynamics. In this paper, we present a method of calculating thermal boundary conductance between heterogeneous materials by integrating first-principles exchange constants and spin–lattice dynamics into the non-equilibrium Green’s function (NEGF) method, and we apply the method to bulk Co and a Co/Cu interface.

Fundamental understanding of carrier scattering and thermal transport is crucial in the design of electronic devices and has long been of interest to researchers. However, theoretical work on interfacial heat transport involving spin–lattice interactions is rare. Zhang et al.⁶ conducted thermal conductance calculations of a hypothetical 1D atomic chain within the framework of the NEGF method and self-consistent Born approximation (SCBA). The authors found that magnon–phonon interactions cause thermal rectification and negative differential thermal conductance in the presence of an external magnetic field at a ferromagnetic–normal insulator interface; these effects could be tuned by adjusting the external magnetic field. Despite the fact that SCBA is computationally intensive and a 1D model was used, the system was not representation of a real material because arbitrary force constants, exchange constants, and magnon–phonon coupling strength constants were chosen.

The purpose of this paper is to develop a method to predict thermal boundary conductance at 3D heterogeneous interfaces with spin–lattice interactions. To ensure that the problem is computationally tractable, we calculate the magnetic properties of bulk Co and a Co/Cu (001) interface with small lattice mismatch from first-principles calculations. The obtained Heisenberg exchange constants are then incorporated into spin–lattice dynamics. We further compute interfacial conductances at different temperatures by calculating autocorrelations of heat flux across the interface sampled at different time steps after the spin and lattice subsystems have reached equilibrium. We also compare the results from the equilibrium molecular dynamics (EMD) and the NEGF methods. Inelastic scattering is considered by attaching Büttiker probes to atoms in the device, and the probe scattering rates are obtained in each contact by fitting their bulk thermal conductivities. The remainder of this paper details the various computational methods employed, their integration, and predictions for bulk and interfacial transport.

Spin–lattice dynamics calculations were conducted with the modified SPIN package⁷ incorporated into LAMMPS.⁸ The Hamiltonian for a spin–lattice coupled system is

Here, the second term is the inter-atomic potential; we use the Embedded-Atom Method (EAM) potential for the system. The third term considers the classical Heisenberg exchange interactions between atoms. The last term is the energy compensation for the spin system at 0 K.

The equation of motion for the spin degree of freedom is derived from a generalized Poisson bracket⁹ as

Equation (2) describes the spin precession about the direction of the local effective field $Hi$⁠. $Si$ is the spin angular momentum associated with atomic site $i$⁠, and the magnitude of $Si$ is proportional to the magnetic moment $μi$ calculated from first-principles. Similarly, the governing equations for updating the lattice degrees of freedom are

The heat flux generated by lattice vibrations is defined as

where $ei$⁠, $vix$⁠, and $τi$ are the sum of kinetic and potential energy, velocity in the x direction, and stress on atom $i$⁠, respectively. If only two-body interactions are considered, the second term in Eq. (5) becomes $−12∑i<j(fij⋅(vi+vj))rijx$⁠, where $fij$ is the force exerted on atomic site $i$ by atomic site $j$⁠. The generation of force is not only due to the variation of inter-atomic potentials but also due to the variation of exchange interactions that depend on inter-atomic distance, as indicated in Eq. (4).

Heat flux generated by spin fluctuations is calculated following the approach used in our previous work.⁵ With an imaginary interface separating the left and right sides of the system, interfacial heat flux is the rate of magnetic energy change on one side of the interface due to interactions with spins on the other side of the interface.

The rate of change in the magnetic potential energy is related to the change in the local magnetic field via

The heat flux induced by spin fluctuations across the interface is then derived from Eq. (6) as

where the $12$ prefactor indicates that the magnetic energy is pairwise and shared by two individual spins.

Thermal conductivity is proportional to the autocorrelation of heat flux at equilibrium and is calculated based on the fluctuation–dissipation theorem,¹⁰ 

The general Green’s function is written as

where $Hd$ is the dynamical matrix in the device region, $Σ1$ and $Σ2$ are self energies in the two contacts, and $ΣBP$ is the Büttiker probe self-energy that incorporates inelastic scattering and takes the form

Energy conservation in each Büttiker probe must be satisfied,

where $Σiin$ represents the in-scattering from the Büttiker probe $i$⁠, $A$ is the spectral function, $Γi$ is the imaginary part of the Büttiker probe self-energy, and $Gn$ is the lesser Green’s function. $Σiin$ and $Γi$ are both block-diagonal; therefore, only block-diagonal terms of $A$ and $Gn$ are needed, and the recursive Green’s function (RGF) method can be used. The avoidance of direct inversions of Green’s functions makes the computation feasible especially for large-scale devices. Because $Σiin$ and $Gn$ depend on Büttiker probe temperatures, Eq. (11) is solved iteratively using the Newton–Raphson technique. Details of derivation and implementation related to the RGF formulation can be found in Refs. 11 and 12. Similarly, phonon temperatures are computed by solving

where $fBE(ω,Tp)$ is the phonon distribution function at equilibrium phonon temperature $Tp$⁠.

Electronic structure calculations were performed using the spin-polarized scalar relativistic Korringa–Kohn–Rostoker Green’s function method (SPR-KKR)¹³ and linear muffin-tin orbital method (LMTO).¹⁴ In both methods, spherically symmetric potentials inside muffin-tin spheres are matched by flat potentials in the interstitial regions at the boundary of muffin-tin spheres, corresponding to a certain energy level. The search of the energy values is simplified by using multiple-scattering theory in the former method and by linearizing basis functions in the latter.

The exchange–correlation was evaluated within the local density approximation (LDA) using the Vosko–Wilk–Nusair functional.¹⁵ The angular momentum expansion cutoff $lmax=2$ and a regular k-mesh of $20×20×2$ were chosen in the self-consistent calculations. The energy convergence threshold was set to be $10−5$ Ry. The atomic sphere approximation (ASA) is used in our current calculations considering that the system under study is closely packed and muffin-tin potentials are close to full potentials.

Figure 1(a) shows the band structure of fcc cobalt with spin up and spin down contributions highlighted in different colors. Figure 1(b) shows a comparison of fcc cobalt band structures calculated from KKR-ASA and LMTO-ASA methods. The band structures produced by these two methods overlap within 1 eV of the Fermi level, confirming the validity of calculations from both methods.

The exchange constants are obtained with the Lichtenstein formula¹⁶ 

where $t$ is single site t-matrix and $τ$ is the scattering path matrix. The computed exchange constants are further fitted to a Bethe–Slater curve as

Lattice dynamics was conducted in a NVT ensemble and the Langevin thermostat was applied to reach the targeted lattice temperature. The spin–lattice dynamics simulation was then run for 400 ps in order for the spin and lattice subsystems to equilibrate. The exchange interaction between spins was included as an inter-atomic distance dependent function, obtained by fitting first-principles data to Eq. (14). Heat flux was calculated within a period of 800 ps based on a NVE ensemble. A 0.2 fs time step was chosen in all the molecular dynamics calculations. Depending on the system of interest and the temperature conditions, the autocorrelation time in Eq. (8) varies.

Figure 2 shows first-principles calculated exchange constant values for fcc cobalt at discrete lattice points, as well as the fitted continuous curve to Eq. (14). Here, $α=18.603$ meV, $δ=1.486$ Å, $γ=1.7×10−5$⁠, and $R=6$ Å. The exchange interaction between first nearest neighbors is the strongest, and the interaction decays exponentially with inter-atomic distances. The cutoff radius was chosen to be 6 Å, beyond which the exchange interaction is negligible.

To obtain a converged thermal conductivity, 30 independent runs for each individual temperature were performed and the results were averaged. Figure 3 shows that the thermal conductivity of bulk Co at 300 K reaches a plateau with an autocorrelation time longer than 10 ps. Magnon thermal conductivity is comparable to phonon thermal conductivity in bulk Co. Magnon thermal conductivity decreases with temperature, as shown in Fig. 4, due to the reduced effective total magnetization. The difference between phonon thermal conductivities calculated with and without spin–lattice interactions is insignificant, especially at high temperatures. The inter-atomic potential is dominant in determining the motion of cobalt atoms. At high temperatures, the effect of exchange interactions is less prominent because atoms are farther from each other than at low temperatures.

To calculate interfacial phonon heat flux, Eq. (5) is modified as

where $fij=ϕijdRj−dJdRij(Si⋅Sj−1)eij$ is the force acting on atom $i$ by atom $j$ and $ϕij$ is the second-order force constant between $i$ and $j$⁠.

Interfacial thermal conductance is calculated using the formula of Puech et al.,¹⁷ 

Thermal boundary conductance at the interface is then calculated by subtracting the ballistic contact resistances from Eq. (16),

Similar to bulk Co, we calculated magnetic properties of the Co–Cu interface from first-principles calculations. As shown in Fig. 5(a), the interface structure is a supercell that consists of eight layers of Co and eight layers of Cu. Periodic boundary conditions are applied in all directions.

In contact with Cu, the magnetic properties of Co atoms next to the interface are modified. As shown in Fig. 5(b), the magnetic moment in Co atoms adjacent to the interface is reduced compared to that in bulk Co (indicated by the horizontal magenta line). Furthermore, the Cu layer adjacent to the interface is magnetized and has a small magnetic moment of 0.001 15 $μB$⁠.

Figure 6(a) shows the exchange interaction constants between Co and Co atoms next to the interface. The nearest-neighbor exchange constants differ for inter-layer and intra-layer Co–Co interactions. The inter-layer Co–Co exchange constant is higher than that in bulk Co while the intra-layer Co–Co exchange constant remains unchanged. Figure 6(b) indicates that the exchange interaction between Co and Cu atoms next to the interface exists because the Cu layer connected to Co layer is magnetized, but the interaction is weak and two orders of magnitude smaller than the interaction between Co and Co atoms.

To clarify the effects of spin–lattice interactions, we calculated interfacial conductances with and without spins. Thermal boundary conductance decreases when magnon heat flux is considered, as indicated by Fig. 7. Although the Cu layer at the interface is magnetized by Co, the magnetization and exchange interaction with Co layers are small. Within the Cu contact away from the interface, magnetic energy is not significantly transmitted. Magnon waves mainly interact with phonon waves on the Co side. The scattering between magnons and phonons on the Co side of the interface is responsible for the reduction of thermal boundary conductance when magnon waves are considered.

The weak effects of magnons on phonons are also observed in the small change of phonon density of states after including magnon effects. In Fig. 8(a), we compare phonon density of states contributed by Co and Cu atoms at the interface with and without spin–lattice interactions. The densities of states were obtained from renormalized force constants at 300 K. In the presence of spin–lattice interactions, density of states slightly changes in the middle range of phonon frequencies, but the change is not significant. Figure 8(b) illustrates the heat flux oscillations with time for cases with and without spins, also showing that the effects of spin–lattice interactions on the motion of atoms are minimal.

The most essential part of solving the NEGF equations involves constructing Green’s function [Eq. (9)]. The dynamical matrix $Hd$ in the device region (harmonic force constants) were obtained from spin–lattice dynamics with EAM potentials and were renormalized to finite temperatures. Büttiker probe scattering rates $τ(ω)$ in Eq. (10) on both sides of the interface were assumed to take the form $τ(ω)=Aω2$⁠, where the temperature-dependent parameter $A$ was obtained by fitting bulk thermal conductivities of Co and Cu obtained from molecular dynamics at corresponding temperatures.

Figure 9 illustrates two fitting processes using the NEGF method. Figure 9(a) calculates the length-dependent resistance while Fig. 9(b) extrapolates the inverse of the thermal conductivity to infinite length. Both methods successfully fit the target phonon thermal conductivity of bulk Co at 300 K (15.1 W/mK) with the same Büttiker probe scattering rates.

By applying a small temperature gradient of 10 K across the device and iteratively solving Eqs. (11) and (12), lattice temperatures in different layers of the device are calculated. Figure 10(a) shows the lattice temperatures in different regions of the device with a temperature gradient of 10 K applied across the Co and Cu contacts and the Co contact being at 300 K. The region on the left of the dashed line (0 Å from the center of the device) represents the Co contact, and the region on the right of the dashed line represents the Cu contact. The difference in bulk thermal conductivities of Co and Cu is larger when magnon effects are considered than when magnon effects are neglected; therefore, the slope on the Cu side is steeper in the presence of magnons.

Spectral heat flux in both contacts has also increased in the middle range of phonon frequencies, as indicated in Fig. 10(b). Even though the total heat flux has increased when magnon effects are included, the accumulated heat flux normalized by the total heat flux changes little in the entire frequency range as shown in Fig. 10(c); therefore, we observe a very small change in the thermal boundary conductance with and without magnon effects as shown in Fig. 11.

Figure 11 also displays similar trends of thermal boundary conductances predicted by the NEGF method and by the EMD method. With spins taken into account, thermal boundary conductance decreases because the interface acts as a barrier for the transmission of magnetic waves incident on the Co side. The handling of inelastic scattering is a major factor that contributes to discrepancies between thermal boundary conductances calculated by the EMD and the NEGF methods. While the EMD method incorporates inelastic scattering through the anharmonicity of the EAM potentials, the NEGF method uses empirical Büttiker probe scattering rates that were obtained by fitting the thermal conductivities of the two contacts (Co and Cu contacts). The domain size effect is another factor that could result in the difference in the EMD and the NEGF calculations. This effect is more significant in the EMD method because the size of the interface in the EMD calculations cannot extend to infinity due to the infeasible computational cost, but the NEGF method inherently treats the two contacts as infinite leads in the calculation of self energies.

In this paper, we have developed a method that combines first-principles calculations, spin–lattice dynamics, and the NEGF method to predict thermal boundary conductance with spin–lattice interactions across a heterogeneous Co–Cu interface. The magnetic moments and exchange constants in the Co layers next to the interface differ from those in bulk Co. The Cu layer adjacent to the interface is also slightly magnetized. However, the exchange interaction between Co and Cu layers is weak such that the magnetic wave incident from the Co side does not efficiently transmit across the interface. The magnons mainly interact with phonons inside the Co contact. The effects of spins on the movement of atoms are small. Both EMD and NEGF predicted a reduction of thermal boundary conductance in the presence of spins because the interface is a barrier to the transfer of magnetic energy.

The effect of spins at the Co–Cu interface is not significant, but the effect could be different for interfaces with stronger exchange interactions. Future work can be performed to find such heterogeneous structures in which spin–lattice interactions play a more significant role.

This work used the Extreme Science and Engineering Discovery Environment (XSEDE), which is supported by the National Science Foundation under Grant No. ACI-1548562, and computational and storage services associated with the Hoffman2 Shared Cluster provided by UCLA Institute for Digital Research and Education’s Research Technology Group.

The data that support the findings of this study are available from the corresponding author upon reasonable request.

The authors have no conflicts to disclose.