Dissociating the phononic, magnetic and electronic contributions to thermal conductivity: a computational study in alpha-iron

Computational tools to study thermodynamic properties of magnetic materials have, until recently, been limited to phenomenological modeling or to small domain sizes limiting our mechanistic understanding of thermal transport in ferromagnets. Herein, we study the interplay of phonon and magnetic spin contributions to the thermal conductivity in α\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\alpha$$\end{document}-iron utilizing non-equilibrium molecular dynamics simulations. It was observed that the magnetic spin contribution to the total thermal conductivity exceeds lattice transport for temperatures up to two-thirds of the Curie temperature after which only strongly coupled magnon-phonon modes become active heat carriers. Characterizations of the phonon and magnon spectra give a detailed insight into the coupling between these heat carriers, and the temperature sensitivity of these coupled systems. Comparisons to both experiments and ab initio data support our inferred electronic thermal conductivity, supporting the coupled molecular dynamics/spin dynamics framework as a viable method to extend the predictive capability for magnetic material properties.

Spin-lattice coupling effects can also be leveraged for dynamic heat flow control in ferromagnetic insulators, enabling the creation of novel nanoscale thermal management strategies [16,17]. Unlike typical insulators where thermal conductivity can be largely inferred solely from the phonons, in ferromagnetic insulators the heat carriers involved are the phonons, electrons, and magnons. Given that magnons are highly sensitive to the magnetic configuration as well as any external magnetic fields, changes in these properties can be exploited to vary the underlying thermal conductivity and regulate heat flow within the material. Spin caloritronic efforts like this can make both future and existing electronic devices more robust, due to a fine control of the temperature gradients within the material [18]. Unfortunately, the ability to directly gauge the coupling between spin and heat currents is limited by a lack of available computational tools.
In the past, methods based on density functional theory (DFT) [19,20] have been the overwhelmingly popular approach to probe the underlying spin structures in micromagnetic materials [21,22]. These methods, although highly accurate, do not scale well with increasing system size and temperature. As such many of these ab initio calculations have been performed at zero temperature [23,24]. Continuum micromagnetic codes can handle much larger temporal and spatial domains but their shortcomings stem from the underlying continuum assumptions which make it difficult to capture complex physics arising from material impurities, material defects, surface anisotropies, and ultrafast magnetization dynamics [25,26].
Alleviating both of the aforementioned shortcomings, atomistic spin-lattice models easily handle complex microstructures and their dynamics while still offering access to spatio-temporal domains that are orders of magnitude larger than what ab initio calculations can handle [25]. In addition to this, atomistic spin-lattice codes naturally couple to MD and DFT codes given each are particle-based computational methods. The ability to simultaneously capture the electronic, phononic and magnetic degrees of freedom within a material remains a significant modelling and simulation challenge.
In the last few years, a number of researchers have focused on incorporating spin-lattice effects into molecular dynamics simulations in this manner [27][28][29][30][31]. These coupled molecular dynamics/spin dynamics (MD-SD) schemes employ the stochastic Landau-Lifshitz-Gilbert equations of motion, which describe how the spin angular velocity evolves in time when the spins are coupled to a thermal bath obeying the fluctuation-dissipation theorem [32,33]. To achieve spin-lattice coupling, a lattice-dependent exchange interaction is incorporated into the magnetic Hamiltonian which describes the total energy of the magnetic spin subsystem. Parameterization of the exchange interaction is typically based on a Bethe-Slater function, which is fitted to ab initio data [31,34,35]. Similar parameterization procedures are employed to capture the phononic degrees of freedom (interatomic potential) wherein machine learned model forms can be trained to reproduce ab initio level accuracy, albeit at a drastically reduced computational cost [36]. The proper combination of quantum accurate interatomic potentials with a complex spin Hamiltonian (i.e., that avoids double counting) solves a grand challenge in magnetic materials modelling and opens the door for a myriad of research outlets [37].
In the present work, we probe the impact of spinlattice coupling on the magnon and phonon thermal conductivities of a-iron in the 300-1200K range. Within this effort, we also employ existing experimental data for the total thermal conductivity of iron to infer the corresponding changes in electronic conductivity with temperature. In addition to the thermal conductivity measurements, a complementary spectral analysis of the phonon and magnon heat carrying modes in iron is carried out. We hypothesize that the spin contribution to thermal conduction is non-negligible, even in electrically conductive metals like iron. Thus, the present novel approach for modelling magnetic materials using coupled MD-SD simulations gives us a first of its' kind opportunity to examine the impact that magnetic disorder has on the heat currents in (ferro-)magnetic materials. mechanical property predictions of the a-phase of iron. The interatomic potential used in the current work is a Spectral Neighborhood Analysis Potential (SNAP), which has been parameterized specifically for MD-SD simulations [37]. As noted in our earlier work, this potential was trained to reproduce the magnetoelastic behavior of a-phase iron (\20 GPa and \1200K) [37]. Additionally, the spin Hamiltonian was parameterized using a set of non-collinear magnetic DFT calculations, which were carried out in VASP [39][40][41][42]. Further details regarding the optimization and training of the SNAP potential and spin Hamiltonian can be found in Nikolov et al. [37]. A brief summary of the MD, SD, and MD-SD schemes is given in appendix A.
We measure the thermal conductivity directly in terms of the heat flux once a steady thermal gradient is established [43,44]. In order to gauge how the thermal conductivity varies throughout the a-phase of iron, we conduct non-equilibrium molecular dynamics (NEMD) [45,46] using triply periodic 20Â20Â240 unit cell replicas ( $ 5.75Â5.75Â68.9 nm with 192,000 atoms), with two 20Â20Â20 hot/cold sections ( $ 5.75Â5.75Â5.75 nm). The hot/cold sections where thermostats are applied are spaced 100 cells ( $ 28.8 nm) apart, an illustration of the computational setup is shown in Figure 1. Thermal analysis is conducted under the pressure-controlled conditions (PCC), as outlined in our earlier work [37]. The PCC scheme assumes that the spin and lattice target temperatures in their respective Langevin thermostats are equal: T l ¼ T s . Thus for brevity from this point forward we drop the spin/lattice notation in the temperature variable and just assume T ¼ T l ¼ T s . To establish a temperature gradient within the bar, the temperature of the hot section is set as T max ¼ 1:08T min where herein T min is varied over the range of 300-1200K.
In order to maintain the cold/hot sections at the target temperatures T min =T max , we apply Langevin thermostats to both the atoms and spins [33]. The region between the hot/cold section is not thermostatted and evolves in the microcanonical ensemble (NVE). The heat flux can be calculated by tracking the energies which are added/removed by either thermostat. For a Langevin thermostat, the net energy added is the sum of energy added by the random forces (f and gðtÞ) and the energy removed by the dampening forces (whose scale is set by c L and k), see appendix A for details. Once the heat fluxes are known, Fourier's law allows us to determine the corresponding thermal conductivities. For a given average temperature T ave ¼ 1:04T min , the lattice constant used to scale the computational box is based on our earlier results for the PCC scheme [37], which takes thermal expansion into account. Since the PCC scheme underpredicts the Curie transition [37], the reported temperatures are rescaled following the procedure outlined by Evans et al. [47].
To equilibrate each system, we run dynamics until the average temperature in the non-thermostated regions stabilizes to T ave ¼ 1:04T min . All simulation geometries start out initially at 0K, it was observed the total equilibration time varies from 200-1000 ps depending on T min , where higher temperatures take longer to equilibrate. Since the material is allowed to thermally expand, in all cases our samples equilibrate to pressures near 0 GPa (p eq \0:1 GPa). Once the temperature gradient in the non-thermostatted section has been equilibrated the simulation is continued for another 1000 ps in order to determine the heat fluxes for both the phonon and magnon sub-systems via their respective thermostats. The dampening constant for the MD thermostat is set to 1 picosecond; meanwhile, the dampening for the SD thermostat is set to 0.01 (dimensionless). The timestep in all simulations is set to 0.0005 ps and the magnitude of the magnetic spin vectors is set to 2.2 Bohr magneton.
The phonon density of states (DOS) and magnon norm/angular velocity DOS spectra calculations are performed using a 20Â20Â20 simulation cell which is thermostated (both atomic and spin thermostats are applied) to the average bar temperature, T ave ¼ 1:04T min for 150 ps. Subsequently, the atomic velocity, spin norms, and precessional velocities are outputted when the thermostats have been removed (NVE dynamics) for every atom/spin at 2.5 fs intervals for 100 ps, leading to a frequency space discretization of 0.01 THz when these time series data are fast Fourier transformed (FFT). The FFT is defined as where FFT Â XðtÞ Ã ðkÞ represents the discrete Fourier transform of a given time series vector X(t). Here, t and k represent the time and frequency domains, respectively, and W n ¼ e ðÀ2piÞ=n is one of n roots of unity. As a consistency check, we also tested a frequency resolution of 0.002857 THz, for which we found no discernible change in the phonon DOS and magnon spectra. We calculate the phonon DOS from the discrete Fourier transform of each atom's velocity components in the x, y, and z directions.
Meanwhile, for the magnons, we look at the frequency domain representation of both the spin orientation ( s ! ) and precessional velocity along the zdirection (x z ), which corresponds to the direction of initial magnetic alignment in the sample. Figure 2 illustrates graphically how the s ! and x ! vectors vary in time. The magnon power spectral density (PSD) is computed from the spin projections Finally, the angular velocity spectrum is defined as Here, M denotes the number of atoms/spins in the sample. Each spectrum is normalized to allow for ease of comparison between different temperatures. The ab initio predictions of the electronic thermal conductivity are carried out using the Kubo-Greenwood formalism [47][48][49][50][51] based on Kohn-Sham orbitals and eigenvalues obtained from DFT [20] calculations which are performed with VASP [39][40][41][42]. The unit cell contains 16 iron atoms with a grid consisting of 4 Â 4 Â 4 k-points on a total of 160 bands. A PAW pseudopotential [52] consistent with the PBE [53] exchange-correlation (XC) functional with a core radius of r c ¼ 1:9a B and containing 16 valence electrons is used in all the DFT calculations. The plane wave cutoff is set to 750 eV and the convergence in each self-consistency cycle is set to 10 À5 . A modified version of the KG4VASP package [54] is used for the post-processing of the electrical conductivity tensor terms. The electronic component of the thermal conductivity is evaluated from the dynamic Onsager coefficients in the DC limit (x ! 0), see Ref. [55] for further details. These calculations were performed on the basis of six random ionic configurations taken from the equilibrated MD configurations which were obtained from DFT-MD calculations performed with VASP [39][40][41][42].
T min T max Temperature Figure 1 Illustration of the computational setup for thermal conductivity measurements. The heat flux for both the lattice and spins is calculated between thermal reservoirs, which in conjunction with the temperature gradient defines the conductivity via Fourier's law. Figure 3a) shows the predicted thermal conductivity obtained from MD, SD, and coupled MD-SD simulations. Comparing phonon data between MD (blue dashed) and MD-SD models (blue and red markers), we infer that the inclusion of spin dynamics does alter the phonon thermal conductivity through the coupling between phonons and magnons, leading to a slight decrease in the phonon thermal conductivity. In the SD case (where atoms do not vibrate), the magnon thermal conductivity (red dashed) at low temperatures (near 380K) is approximately three times larger than the MD-SD magnon conductivity (red markers). As the temperature increases, the SD magnon thermal conductivity rapidly decreases, approaching values much closer to the MD-SD magnon conductivities. Based on this we can infer that the impact of lattice motion tends to be more significant at temperatures below 700K where only mild spin fluctuations are present (m [ 0:88), leaving magnon-phonon scattering to explain the sharp decrease in conductivity. At temperatures above 700K spin fluctuations become significantly higher as the net magnetization in the simulation cell approaches zero and the significance of lattice vibrations becomes weaker. For the MD-SD case, the magnon thermal conductivity is larger than the phonon thermal conductivity in the range of \700K. At higher temperatures, the magnon conductivity decreases below that of the phonons and near the Curie temperature it approaches values close to zero (j m;MDÀSD % 1 À 2 W/mK).    . 3b) shows the MD-SD thermal conductivity results along with a complementary experimental data set (black markers) from Fulkerson et al. [56]. In order to better compare MD-SD results with the data from experiment, we assume that the total thermal conductivity measured in experiment can be expressed in terms of an electronic and lattice thermal conductivity as j T ¼ j l þ j e , where the lattice thermal conductivity is composed of the phonon j p and the magnon j m thermal conductivities, j l ¼ j p þ j m . By using the experimental values for the total thermal conductivity j T;exp , we can use the phonon and magnon MD-SD thermal conductivities to infer an electronic thermal conductivity for our MD-SD simulations as: j e;MDÀSD ¼ j T;exp À j p;MDÀSD À j m;MDÀSD .

Results
The inferred j e;MDÀSD (green markers) is compared against the experimental data in Fig. 3b). Our inferred electronic thermal conductivity is in good agreement with the experimental data which was obtained using electrical resistivity measurements. Moreover, the inferred lattice thermal conductivity which Fulkerson calculates using the Backlund et al. relation [57] is in good agreement with the MD-SD lattice conductivity data. Fulkerson et al. experimentally shows that past the Curie temperature the electronic thermal conductivity data experiences a change in slope, as it begins to flatten out. This trend is also observed in our MD-SD simulations due to the loss of magnon thermal conductivity. As Buckland points out the inverse of the electronic thermal conductivity gives a measure of the electronic thermal resistivity, which can be separated into terms due to thermal scattering (W eT ), s-d scattering (W eS ), and impurity scattering (W eI ) [57]. Since at low temperatures the thermal scattering term dominates [57], we estimate that W eT % 2 cm-K/W, which is 33% higher than the value found by Backlund [57].
Note that the current SNAP potential is not trained to reproduce the bcc-fcc transition, so the discontinuity observed experimentally in the lattice thermal conductivity near 1180K cannot be reproduced with the current interatomic potential. However, this can be easily fixed in the future by including additional DFT training data [37].
The orange markers in Fig. 3b) represent the electronic contribution to the thermal conductivity obtained from ab initio calculations using the Kubo-Greenwood formalism. These calculations have been successfully used to evaluate the behavior of iron [58] and its alloys especially for the transport properties under earth-core conditions [59,60]. The Kubo-Greenwood evaluations are performed at the MD-SD predicted lattice constants. Fig. 3b) shows the DFT calculations are in great agreement with experimental data. Moreover, above 400K, the DFT calculations are also in reasonable agreement with the MD-SD data. Interestingly, below 400K, the MD-SD data overestimate the lattice conductivity (j p;MDÀSD þ j m;MDÀSD ) which in turn leads to a decrease in j e;MDÀSD .
To understand the temperature dependence and potential scattering mechanisms of either heat carrier type in our MD-SD simulations, we turn to a spectral analysis of the heat carrying modes in iron. The plot in Fig. 4a) shows phonon DOS data for three different temperatures (solid blue, green, and red) along with an experimental data set (black markers) gathered using a nuclear resonant inelastic X-ray scattering (NRIXS) technique at ambient conditions (values scaled for clarity) [61]. The experimentally predicted acoustic and optical peaks are near approx. 5.6 and 8.5 THz, respectively, and are matched reasonably well by the MD-SD data. The MD-SD simulations tend to overestimate the acoustic band edge peak by about 0.4 THz or 1.65 meV and do not show the small optical peak near 6.6 THz, observed by Minkiewicz [62] and Hu [63]. Increasing the frequency resolutions from 0.01 to 0.002857 THz did not lead to any discernible changes in the DOS curves. For the MD-SD results, as the temperature increases the DOS broadens as the acoustic peaks shift to lower energies with thermal expansion, consistent with previous findings [64][65][66]. The inset in Figure 4a) shows how the first peak (acoustic band edge) shifts at higher temperatures. Overall the MD-SD simulations do a good job of estimating the frequency of the acoustic peak within the range of 300-1000K indicating the PCC scheme captures thermal expansion at higher temperatures. Figure 4b) shows how the MD-SD magnon PSD in the frequency domain varies at different temperatures. As the temperature increases a uniform red shift and broadening of the magnon modes is observed, consistent with prior studies [67]. As the Curie temperature is approached, and the magnon thermal conductivity decreases, many higher frequency peaks diminish, owing to a loss of short ranged spin ordering. The inset of Figure 4b) shows the high frequency decay of the magnon PSD. It can be seen that the magnons remain only weakly identifiable at higher temperatures compared to the phonons [29] and only at low temperatures is a peak observed near 46-47 THz.
To better compare between the MD-SD, MD, and SD cases, we also examine how the phonon DOS and magnon spectra (U) compare at a temperature of T ¼ 570K (m ¼ 0:786), shown in Fig. 5. Fig. 5a) shows that while the phonon acoustic peaks remain at the same frequency, a small redshift of about 0.24 THz in the optical band edge peak is observed. The inclusion of spin dynamics also seems to dampen out excitations throughout the 0-9 THz range. Conversely, comparing the magnon spectra, U, between the MD-SD and SD data, no shifts at low frequencies (\15THz) are observed, as exemplified in Fig. 5b). However, at high frequencies, the addition of lattice motion into the dynamics seems to strongly dampen the high frequency peak near 41 THz.
After examining the spin orientation ( s ! ) spectra U, we turn to the spectra defined by the precession velocities of individual spins ( x ! ). Fig. 6 shows how the PSD of the precessional velocity in the z-direction (thermal gradient) (x z ) varies between MD-SD (solid green) and SD (solid orange) simulations for m ¼0.928. As can be seen for the isolated spins (SD), where the lattice is fixed, H reverts to a profile similar to U. Meanwhile, for the coupled MD-SD case, the phonon DOS is noticeably super-imposed in the x z signal. While this is noticeable along the longitudinal direction (z), it was found that the transerve directions (x and y components) of the precessional velocity, x ! , do not show this effect. As Figure 2 shows, the s ! vectors vary rapidly in time, showing a weak coherence with neighboring vectors, whereas the x ! vectors are more tightly bound to the z-axis and exhibit much stronger coherence with their neighbors.
What is particularly interesting in Fig. 6 are the pronounced magnon peaks below 5THz, which persist even when the lattice is thermalized (MD-SD). A non-magnetic simulation of iron would have only yielded a smooth increase in the phonon DOS up to the acoustic band edge, but this is now decorated with vibrations that show strong overlap in the spin procession and phonon spectra. In fact, many of the distinct peaks in the MD-SD vibrational spectra (except the optical band edge) contain a corresponding peak in the magnon spectra, highlighting strong magnon-phonon coupling in this ferromagnetic material.
The inset of Fig. 6 shows how the x z PSD varies with temperature. Interestingly, the strong magnonphonon coupled modes are not picked up by x z signal near the Curie temperature. The attenuation of the phonon DOS in H is indicative of the weakening spin-lattice coupling (relative to thermal fluctuations) at higher temperatures. This effect can be inferred also from Callen's law which describes the decay of the anisotropy energies (which largely arise from spin-orbit coupling) at higher temperatures [68]. At higher temperatures, we also observe the same red shift and peak broadening detected in U. No discernible peaks for H are found at frequencies above 20 THz. This contrasts the data for U where at low temperatures a high frequency peak near 41-43 THz can be observed. At T=950K (m=0.40), while still below the Curie temperature, the superposition of phonon vibrations in the magnon spectra is lost, indicating that at this elevated temperature the two sub-systems of heat carriers are decoupled. However, this separation of the lattice and magnetic degrees of freedom is not absolute, Fig. 3a) shows how the SD magnon conductivity is actually below the MD-SD conductivity near 950K indicating phonon driven thermal transport in the magnon sub-system. Furthermore, in all temperatures studied here, the phonon conductivity is lower in the coupled MD-SD case relative to the simulation conditions where spins are frozen. In effect, while the additional scattering between the two subsystems decreases the overall phonon conductivity, the strong magnon-phonon coupling gives rise to additional heat carrying spin wave modes that are mediated by persistent phonon modes at high temperatures. Furthermore, it is observed that coherency in the spin state (m [ 0:8) plays a drastic role in the ability for a material to conduct heat through spin waves, which is evidenced by the rapidly decaying features in the magnon spectra in Figs. 4, 5, 6.

Conclusion
The impact of phonon-magnon coupling on the thermal conductivity of iron was quantified up-to and beyond the Curie temperature in ferromagnetic iron. It was found that for the phonon contribution, the inclusion of spin dynamics leads to a small decrease in the thermal conductivity throughout the entire temperature range. Meanwhile, the inclusion of lattice vibrations significantly decreases the magnon thermal conductivity up to 700K, but remains non-negligible up to the Curie temperature. Both of Meanwhile, for the MD-SD data, the optical peak is shifted by 0.23 THz from the corresponding frozen spin, MD peak (at 8.26 THz). b) For T ¼ 570 K, the magnon spectra in the low frequency range (0-12 THz) does not change between the SD and MD-SD simulation conditions. However, as shown in the inset plot, near 40 THz the MD-SD magnon spectra does experience a strong attenuation, indicating that these high frequency vibrations are sensitive to lattice motion.  Figure 6 Comparison between MD-SD and SD angular velocity spectra for T ¼ 570 K. In the SD simulations, there are no lattice vibrations so the phonon DOS is not reflected in the angular velocity spectra H. Inset plot shows MD-SD data for H at different temperatures. At low temperatures, the phonon DOS peak is reflected in the angular velocity spectra, indicating that the precessional motions of the spins are reflective of phonon vibrations. Upon approaching the paramagnetic state the phonon DOS peaks (in the 5 to 9 THz range) vanish. these observed effects are due to the added scattering mechanisms that each of the two sub-systems provides to the other. However, it was observed that the spin contribution to thermal conductivity near the Curie temperature remains non-zero strictly due to spin wave coupling to lattice vibrations.
Utilizing the calculated lattice and spin conductivities along with experimental measurements of the total thermal conductivity, an inferred electronic conductivity of iron is obtained in good agreement with both experimental findings and ab initio calculations. While the classical MD simulations are unable to capture these electronic effects, it is encouraging that our ab initio trained parameterization reproduces a proper partitioning of the thermal conductivity among the phonon, magnon and inferred electronic degrees of freedom.
To supplement the observed trends in thermal conductivity, the phonon and magnon spectra were analyzed to obtain a deeper understanding of how temperature impacts each heat carrying mode. Our findings show that for the phonon DOS, our predictions are in good agreement with nuclear resonant inelastic X-ray scattering data, including peak red shifting with increasing temperature. As expected, approaching the magnetic order-disorder transition, peaks in the magnon spectra become largely featureless. From these spectral results, we were able to highlight changes in the spectra that are due to strong magnon-phonon coupling which helps interpret changes in the coupled MD-SD conductivity relative to the decoupled predictions. A shift in the phonon optical band edge due to spin contributions and the loss of the high-frequency magnon mode where lattice vibrations are present, exemplify these spectral changes.
The efforts described here outline a basis for characterizing the impact spin currents have on the heat flux in any magnetic material. These atomistic simulation methods and analyses can be easily applied to more complicated systems featuring grain boundaries, multiple magnetic domains, defects, and additional alloying elements. Thus, a more mechanistic and holistic view of thermal conductivity in this class of metals can be obtained within one set of simulation tools. Our future work is aimed at utilizing these methods to explore these complexities as well as phase transitions in magnetic materials. results and analysis. Any subjective views or opinions that might be expressed in the paper do not necessarily represent the views of the U.

Code Availability
The code which was used to train the SNAP potential used herein and reported in Ref. [68] is available from https://github.com/FitSNAP/FitSNAP.

Declarations
Conflict Of Interest Statement The authors declare that they have no known competing financial interests that could have appeared to influence the work reported in this paper.

Open Access
Here, the variables a and d are the interaction energy and interaction decay length and c is a dimensionless curvature parameter. Meanwhile, Ç R c À r ð Þ is Heaviside step function which switches at the radial cutoff R c . The corresponding values of a, d, and c for J r ij À Á and K r ij À Á along with the value of R c for the exchange interaction are shown in Table 1.
As shown in equation (10), the MD-SD Hamiltonian incorporates H SD from equation (8). The second term on the right side of the equation represents the kinetic energy and the third term (V SNAP ðr ij Þ) arises from the interatomic SNAP potential. The variables p and m i are the particle momentum and mass of atom i, respectively.
ðx i þ gðtÞÞ Â s i þ ::: Equations (11)-(12) describe how the derivatives of position and momentum change with time. Equation (13) is the Landau-Lifshitz-Gilbert equation which describes spin vector precession in a spin subsystem coupled to a thermal bath. The variable x i is the precessional angular velocity whose definition is based on the work of Evans et al. [74]. The definition for x i is given in equation (14) below, where l i is the atomic spin norm (set to 2.2 Bohr magneton) and b is the gyromagnetic ratio (b % 0.176 radÁTHzÁT À1 ).
The variable c L in equation (12) is the Langevin damping constant for the lattice which is coupled to the fluctuating force f (which follows Gaussian hf a ðtÞf b ðt 0 Þi ¼2k B T l c L d ab dðt À t 0 Þ ð 16Þ The variables k B and T l are the Boltzmann constant and target lattice temperature, respectively [33]. Similarly, in equation (13), we have the variable gðtÞ which is another fluctuating force (this time for the spin subsystem) that follows Gaussian statistics. Again here gðtÞ is coupled to the transverse dampening constant k via the fluctuation dissipation relation, shown in equations (17)-(18) [33].
hgðtÞi ¼0 ð17Þ Here, T s is the target spin temperature and h is Plank's constant.

B Computational performance
Simulations for the MD-SD, MD, and SD schemes were carried out using 20 dual-socket E5-2695 v4 @ 2.10GHz nodes (36 cores per nodes). The DFT-MD simulations, meanwhile, were carried out using 4 Tesla V100 GPUs. The performance metrics for our simulations are shown in Table 2. Due to the strong linear scaling of ML-IAP potentials like SNAP, we can greatly accelerate computational performance by using a large number of nodes/processors, thus keeping the atom/cpu-core count low [36]. This degree of parallelization is not typical for potentials like EAM or Lennard-Jones where a breakdown of 30,000 atoms /cpu-core is preferred. The atom/cpucore count in MD-SD and SD simulations is limited (to approx. 200-250 atoms/cpu-core) by the sectoring algorithm which is used within the Suzuki-Trotter decomposition of the spin equations of motion. Table 2 also shows that for the MD-SD and SD methods, the inclusion of spin dynamics leads to an approx. 40% drop in computational performance. This drop in performance is expected and can be attributed largely to the Suzuki-Trotter decomposition which requires quarter timestep updates. Meanwhile, the standard integration scheme in LAMMPS relies on the Velocity-Verlet algorithm which uses a half timestep updates. Even though the lattice is frozen in the SD scheme, the mechanical contribution from the SNAP potential is still computed in order to obtain the correct pressures/energies. For the DFT-MD method the computational performance scales as N 3 , whereas for the classical schemes a scaling of N 1 is observed, consistent with previous studies [36,75].