Intrinsic energy flow in laser-excited 3$d$ ferromagnets

Ultrafast magnetization dynamics are governed by energy flow between electronic, magnetic, and lattice degrees of freedom. A quantitative understanding of these dynamics must be based on a model that agrees with experimental results for all three subsystems. However, ultrafast dynamics of the lattice remain largely unexplored experimentally. Here, we combine femtosecond electron diffraction experiments of the lattice dynamics with energy-conserving atomistic spin dynamics (ASD) simulations and ab-initio calculations to study the intrinsic energy flow in the 3d ferromagnets cobalt (Co) and iron (Fe). The simulations yield a good description of experimental data, in particular an excellent description of our experimental results for the lattice dynamics. We find that the lattice dynamics are influenced significantly by the magnetization dynamics due to the energy cost of demagnetization. Our results highlight the role of the spin system as the dominant heat sink in the first hundreds of femtoseconds. Together with previous findings for nickel [Zahn et al., Phys. Rev. Research 3, 023032 (2021)], our work demonstrates that energy-conserving ASD simulations provide a general and consistent description of the laser-induced dynamics in all three elemental 3d ferromagnets.


I. INTRODUCTION
Ultrafast manipulation of magnetic order with light promises pathways to new applications in magnetic data storage and spintronics [1].Femtosecond laser excitation can change magnetic order in various ways; for example, it can induce ultrafast demagnetization [2], switch the magnetization direction [3,4], and induce spin reorientation [5].The microscopic mechanisms governing the response of magnetic materials to laser excitation continue to be a topic of current research [6][7][8][9][10][11][12][13][14].An important factor governing the response of a material to laser excitation is the intrinsic energy flow between electronic, magnetic, and lattice degrees of freedom.
When Beaurepaire et al. discovered ultrafast demagnetization in Ni, they introduced a phenomenological three-temperature model (3TM) to describe the observed magnetization dynamics [2].While the 3TM offers an intuitive explanation for the observed dynamics, recent studies suggest that it falls short of a full description of ultrafast demagnetization.In particular, there is experimental and theoretical evidence that the spin system is not in a thermal state on ultrafast timescales [12,15,16], suggesting that a more detailed description of the magnetic degrees of freedom is necessary.
To obtain a full quantitative description of a material's response to laser excitation, any proposed model must be verified by comparison to experimental data of the responses of electronic, magnetic, and lattice degrees of freedom.The lattice plays a major role in the dynamics of 3d ferromagnets, since it drains energy from the electrons via electron-phonon coupling on similar timescales compared with the demagnetization, thus reducing the temperature of the electron system.On the other hand, lattice dynamics are also influenced by magnetization dynamics, even if the coupling is only indirect via the electron system.Our previous work on Ni demonstrated that energy flow into and out of the spin system leads to a significant slow-down of the lattice dynamics [17].This suggests that accounting for this energy flow is integral to any model quantitatively describing the responses of all three subsystems in 3d ferromagnets.
Despite their significant role in the energy flow dynamics, the lattice dynamics of 3d ferromagnets are less studied compared with electron and spin dynamics [6,10,11,16,18,19].Time-resolved diffraction offers the most direct way to study lattice dynamics since it is only sensitive to the lattice.Hitherto, only two studies of the sub-picosecond lattice dynamics of Co or Fe with time-resolved diffraction exist [20,21] and neither of them focuses on the lattice heating in the ferromagnet.Furthermore, literature values for the electron-phonon coupling parameter G ep vary significantly, from 6 × 10 17 to 4.05 × 10 18 W m 3 K for Co [6,[22][23][24] and from 7 × 10 17 to 5.48 × 10 18 W m 3 K for Fe [23,[25][26][27][28][29].In addition, there are several literature values for the electronphonon coupling parameter λ, which is related to G ep (see, for example, Ref. 30) and also varies significantly [31][32][33][34].
In ferromagnets, extracting the electronphonon coupling solely from experiments is particularly challenging because three different subsystems contribute to the observed dynamics.
Here we measure the lattice dynamics of Co and Fe directly using femtosecond electron diffraction.Instead of extracting G ep from experiments, we perform spinresolved density functional theory (DFT) calculations, which yield G ep as well as the heat capacities of the electrons and the lattice [17,35].
Based on the experimentally measured lattice dynamics and the DFT results, we study the intrinsic energy flow between electronic, magnetic, and lattice degrees of freedom.We employ energy-conserving atomistic spin dynamics (ASD) simulations [17,36], a hybrid model which combines conventional ASD simulations with a description of the energy flow between all subsystems.By directly simulating the evolution of the spin system, ASD simulations have the advantage that they are not constrained to thermal descriptions of the spin system.Previously, we applied this approach to Ni with excellent agreement between theory and experiment [17].Here we demonstrate that the same considerations hold also for Co and Fe, thus generalizing our approach to all three elemental 3d ferromagnets.To demonstrate the strong influence of the magnetization dynamics on the lattice dynamics, we compare results of the conventional two-temperature model (TTM), which does not consider the spin system, to results of the energy-conserving ASD simulations.With the latter, we obtain excellent agreement with the lattice dynamics of Co and Fe as well as a good description of the magnetization dynamics.This demonstrates that ASD simulations offer a consistent description of the laser-induced dynamics in all three elemental 3d ferromagnets.
In Section II, we describe the experiment and the data analysis.Section III presents both experimental results for the lattice dynamics and model results.Based on the ASD simulation results, in Section IV we discuss the intrinsic energy flow between electrons, spins, and the lattice in detail.Section V summarizes the main findings.

II. TIME-RESOLVED DIFFRACTION EXPERIMENT
The samples of our experiments are freestanding thin films of Co or Fe with a thickness of 20 nm, sandwiched between 5 nm-thick layers of silicon nitride.They were grown on a single crystal of NaCl by magnetron sputtering, resulting in polycrystalline films.Next, they were transferred onto a transmission electron microscopy (TEM) grid by floating the films on water.The samples were not subjected to a magnetic field before the experiment, i.e., different magnetic domains are likely present in the sample.
To directly access the lattice dynamics after laser excitation, we employ femtosecond electron diffraction using the setup described in Ref. 38.A schematic illustration of the experiment is presented in Fig. 1(a).In the electron diffraction experiment, the thin films are excited with an ultrashort laser pulse.The lattice response to laser excitation is probed using an ultrashort high-energy electron pulse.The electrons diffract off the sample and are recorded in transmission.The electron energy was 70 keV for the experiments on Co and 60 keV for the experiments on Fe.All experiments were performed at room temperature (295 K).Since the samples are polycrystalline, the diffraction patterns consist of Debye-Scherrer rings, as shown exemplarily in Fig. 1(a) for our Co sample.Our main observables are changes in the intensities of the diffraction rings following laser excitation.These are directly related to the change in atomic mean-squared displacement (MSD) [39]: Here, q is the scattering vector of the diffraction ring is the intensity as a function of pumpprobe delay, and I 0 is the intensity before laser excitation.
To extract the MSD dynamics from the diffraction patterns, we employ a global two-step fitting routine [37].In brief, the first step is a fit to the diffraction pattern before laser excitation.
The fit function consists of a background function plus Lorentzians for the diffraction rings, all convolved with a Gaussian to account for the finite coherence of the electron beam.In the second step, the time-dependent changes are extracted.For this, we fix most parameters of the fit function and allow only changes of the lattice constant (i.e., expansion/contraction of the lattice), changes of the MSD, and changes of the background parameters.The lattice dynamics are extracted from the full diffraction pattern instead of individual diffraction rings, which increases the reliability of the results.Further information on the global fitting routine is available in Ref. 37.

A. Experimental results for the lattice dynamics
Experiments were performed on Co and Fe for several excitation densities each.For every excitation density, several delay scans were recorded and the results were Here the background-subtracted pattern is shown for illustrational purposes; however, note that in the analysis of the diffraction patterns, fits are performed to background and rings simultaneously.The lower part shows the differences of the RPs compared with the RP before laser excitation for several pump-probe delays (solid curves).The dashed black lines show the fit results of the global fitting routine described in Sec.II and in detail elsewhere [37].(c) Same as (b) but for Fe.
averaged before applying the two-step fitting routine.Examples for the resulting MSD dynamics of Co and Fe are presented in Figs.2(a) and 2(b), respectively.For the conversion of MSD to lattice temperature, we calculated the temperature-dependent Debye-Waller factors for Fe and Co based on the phonon density of states (DOS) from DFT (see Appendix A).We performed fits to the experimental data using a single exponential function, convolved with a Gaussian of 250 fs (FWHM) to account for the time resolution.The time constant of the single exponential function, the amplitude, and the onset (time zero) were fit parameters and the fit range was from -0.5 to 4 ps.The results for the time constant are shown in Fig. 2(c) for different excitation densities.For Co, we find that the time constant increases with increasing excitation density.For Fe, no clear trend is observed.

TTM
In the next step, our goal is to analyze the intrinsic energy flow between electronic, magnetic, and lattice degrees of freedom.
For this, we compare our experimental data to models for the energy flow.In order to minimize the number of free parameters in the models, we use spin-resolved DFT to obtain the (electron-temperature-dependent) electron-phonon coupling parameter as well as the electron and lattice heat capacities.The results for the heat capacities and the electron-phonon coupling parameters are presented in Figs.3(a) and 3(b).All electronic heat capacity and G ep curves presented in Fig. 3 are the sum of majority and minority carrier contributions.Details about the DFT calculations are described in Appendix A.
Having obtained the input parameters for the models from DFT, we start with the conventional TTM [40, 41], which considers only the electronic and lattice degrees of freedom but disregards the spin system.The system is modeled as two heat baths, electrons and the lattice, which are coupled by the electron-phonon coupling parameter G ep .The evolution of the electron temperature (T e ) and the lattice temperature (T l ) is then described by two coupled differential equations: Here c e and c l are the electronic and lattice heat capacities, and P (t) is the source term, i.e., the energy input to the electronic system due to the laser excitation.The laser excitation is modeled as a Gaussian with a FWHM of 80 fs.Its maximum (time zero) is determined from the single exponential fits described earlier.The energy deposited by the laser is determined from the lattice temperature after electron-lattice equilibration (in the range from 1.5 to 4 ps after laser excitation) and the heat capacity (sum of electron and lattice contribution).Hence, there are no fit parameters in this TTM.The comparison between the TTM and the experimental results for the lattice dynamics is shown in Fig. 4 for both materials and several fluences each (dashed curves).We find that for both Fe and Co, the lattice temperature rise predicted by the TTM is faster compared with our experimental results.This finding agrees with previous results on Ni [17].A major source of this disagreement is the fact that the TTM does not consider magnetic degrees of freedom.Therefore, also the energy associated with magnetization dynamics is neglected.However, as we showed previously for the case of Ni, energy flow into and out of magnetic degrees of freedom has a profound influence on lattice dynamics [17].Hence a model that takes the spin system into account is needed.

ASD simulations
In order to include the spin system in our model of the energy flow dynamics, we use energy-conserving ASD simulations, which simulate the dynamics of the spin system based on a Heisenberg model and the stochastic Landau-Lifshitz-Gilbert (s-LLG) equation.The coupling of electron and phonon system is described with a TTM based on the DFT results [see Figs.3(a) and 3(b)], as in the previous subsection.Energy conservation is achieved by monitoring the energy content of the spin system and subtracting/adding the change in spin energy from/to the electron system at each time step of the simulation.The TTM equation for the electron temperature (Eq. 3) is thus modified as follows: Here ∆E s corresponds to the change of spin energy in the time step ∆t.It is calculated as follows: Here the S i are the individual spins of the ASD simulation and the factor s 2 / [s(s + 1)] accounts for the quantized nature of the spins (s ≈ 3 2 for Co and s ≈ 2 for Fe).Note that direct spin-lattice coupling is not included in the model.
The fast demagnetization timescales in 3d ferromagnets suggest that the magnetization dynamics are dominated by electron-spin coupling.In nickel, spin-lattice coupling was estimated to be an order of magnitude smaller than other coupling constants [2].More details about the energy-conserving ASD simulations are described in Ref. 17 and the material-specific simulation parameters for Co and Fe are stated in Appendix B. With this model, both the nonequilibrium spin dynamics and the energy flow between electrons, spins, and the lattice can be described.
The coupling between electrons and spins in the ASD simulations is governed by the damping parameter α.It determines how fast the spins react to the stochastic field of the s-LLG equation, whose amplitude in turn depends on the electronic temperature.Here, we use α = 0.01 for Co and α = 0.005 for Fe, which yield a good description of the experimentally measured lattice dynamics at low excitation densities as well as realistic magnetization dynamics.These values are in good agreement with recent experimental results for α [42].
Figure 4 presents the ASD simulation results for both Co and Fe.First, we focus on the results for Co, shown in Figs.4(a)-4(e).We find excellent agreement with the experimentally measured lattice dynamics for all excitation densities.Clearly, the agreement is much better than that obtained with the TTM.This finding highlights the importance of considering energy flow into and out of magnetic degrees of freedom, in agreement with our previous results for Ni [17].
The ASD simulation results for the magnetization dynamics of Co are presented in Fig. 4(e).The general shape of the magnetization dynamics, in particular the pronounced drop and relatively fast recovery of the magnetization, agrees well with recent experimental results [24,43].
Regarding the demagnetization dynamics in the first hundreds of femtoseconds, the ASD simulation results reach the minimal magnetization roughly 100-200 fs faster than in reported experiments [6,24,43].
This could be due to deviations of the electronic distribution from a Fermi-Dirac distribution at early times after laser excitation, and due to the phenomenological electron-spin coupling in the ASD simulations.In addition, the ASD simulations describe an idealized system without defects or surface effects and assume homogeneous excitation, which can also contribute to the observed discrepancies.Regarding the magnetization recovery, we observe good agreement with results from Ref. 24 while the recovery measured by Refs.6 and 43 is slower than the ASD simulation results.
It should be noted that there is some spread in the experimental results for the magnetization dynamics, even when only thin films on non-metallic substrates are considered [6,24,[42][43][44][45].On short timescales, the measured results can contain artifacts from statefilling effects when probing optically [46,47].On longer timescales, magnetization dynamics can be influenced by transport effects (of electrons and phonons out of the probed region), which depend on the sample geometry.Also other macroscopic sample properties may play a role in the magnetization response.A recent study found differences in the ultrafast response depending on the orientation of the magnetization relative to the crystal lattice [24].In principle, both the demagnetization and the magnetization recovery contain valuable information on the coupling strength between electrons and spins.For example, reducing α in the ASD simulations leads to a slower demagnetization but also to a less pronounced magnetization recovery because the spin system heats less (and thus absorbs less energy) during the time when the lattice is still cold.A more precise comparison of model results to the responses of all subsystems could be obtained by measuring the lattice, magnetization, and electron dynamics on identical samples.
Next, we focus on the ASD simulation results for Fe, shown in Figs.4(f)-4(j).For low fluences, we obtain excellent agreement with the experimentally measured lattice dynamics, again corroborating the strong influence of the spin dynamics on the lattice dynamics.However, the quality of agreement is not as high as for Co.Specifically, for high fluences, the simulations predict lattice dynamics that are slower than the experimental observations.
In the following, we discuss possible reasons for these deviations.In our ASD simulations, the strength of the electron-spin coupling, described by the damping parameter α, is constant.At higher excitation densities, however, the electron-spin coupling could react to the laser-induced changes of the electronic structure.Since Fe has the largest spin heat capacity of all three elemental 3d ferromagnets at room temperature in combination with a rather low electronic heat capacity, its lattice dynamics are most sensitive to energy flow into and out of the spin system.Therefore, it is plausible that deviations between ASD simulations and experiments performed at high fluences are larger for Fe compared with Ni or Co.Furthermore, transient nonthermal electron and phonon distributions could contribute to the observed lattice dynamics for both Fe and Co [23,48].Experimentally, we observed a small apparent shift in time zero by tens of fs for high excitation densities.This could be caused by electron thermalization, which is more efficient at high excitation densities and typically enhances energy transfer to the lattice [49].Nonthermal distributions of electrons and phonons are not accounted for by our models and including them might change the optimal α toward lower values.In addition, direct spin-lattice coupling is not included in our model, as mentioned above.Even though we expect this coupling to be weak, it constitutes another channel for energy flow to the lattice and could enhance in particular the energy flow out of the spin system.Finally, DFT calculations are ground-state calculations.After laser excitation, band structure changes (for example a transient reduction of the exchange splitting) can occur [12], which lead to changes of the electronic heat capacity and the electronphonon coupling, especially for higher fluences.Hence, ASD simulations are expected to be most accurate for low excitation densities in general, which we observed also for Ni [17].Nevertheless, for low and moderate fluences, our ASD simulations offer an excellent description of the laser-induced lattice dynamics for all three 3d ferromagnets.
Regarding the magnetization dynamics of Fe, the ASD simulation results are presented in Fig. 4(j).The initial demagnetization rate agrees well with experimental results [43].For the magnetization recovery, different results are reported in the literature [10,16,19,42,43,47,50], from very little or no recovery [43,47] to almost complete recovery [50] on few-picosecond timescales.Only thin films on non-metallic substrates are considered here, which are expected to have the least transport effects.Due to the large spread of literature results, as in the case of Co, a more precise comparison of the model to the results of all three subsystems would require measuring their responses on identical samples.Based on the available experimental data, we conclude that our simulations provide a realistic description of the magnetization dynamics.Energyconserving ASD simulations thus offer a description that is consistent with the responses of the lattice and the magnetization, which is an important step toward a complete, consistent description of the laser-induced dynamics of 3d ferromagnets.

IV. DISCUSSION
The good agreement of the ASD simulations with our experiments and the disagreement of the TTM show that energy flow into and out of the spin system has a significant impact on the lattice dynamics of Co and Fe.Based on the ASD simulation results, we are now able to analyze the intrinsic energy flow between electronic, magnetic, and lattice degrees of freedom in detail.The distribution of the absorbed energy between the three subsystems is presented in Fig. 5.
After laser excitation, the total energy in the system stays constant, which is visualized by the black curve.From then on, energy is only redistributed between the different degrees of freedom.The laser pulse excites the electrons (blue curve), which initiates the energy flow from the electrons to the spin system (green curve).Already shortly after excitation, the spin system contains more of the additional energy than the electron system.Once spins and electrons have equilibrated and the electrons cool down further due to electron-phonon coupling, energy starts flowing back from the spin system to the electrons.In addition, energy also flows from the electrons to the lattice, such that in total, energy flows out of the electron system, although at a lower rate than during the demagnetization.Finally, thermal equilibrium is established after several picoseconds.Similar to our previous results for Ni, we find that also for Fe and Co, the ASD simulations predict a nonthermal spin system on short timescales after laser excitation.This is presented in Fig. 6.The additional spin energy in the system is shown as solid curves (the absorbed energy densities are the same as in Fig. 4).In addition, the dashed curves show how a thermalized spin system would behave.The thermalized case is based on the equilibrium properties of the spin system and the magnetization dynamics from the nonequilibrium The additional energy density ∆E in the system is displayed.After laser excitation, the total energy (black) stays constant and energy is redistributed between electronic (blue), magnetic (green), and lattice (red) degrees of freedom.
simulation.We use the equilibrium relationships between magnetization and spin energy, shown in the insets of Fig. 6, to translate the magnetization dynamics from the simulations into spin energy dynamics.Comparing this result to the spin energy dynamics obtained directly from the simulations allows to identify deviations from thermal behavior: Whenever the two quantities do not coincide, the spin system is in a nonthermal state.On short timescales below ca. 1 ps, dashed and solid curves differ, which indicates that the spin system is in a nonthermal state during this period.This nonthermal state is characterized by a relatively high spin energy content compared with the demagnetization amplitude, as the comparison between dashed and solid curves directly shows.This is analogous to our ASD simulation results for Ni and indicates that relatively many spin excitations with significant misalignment of neighboring spins are present compared with thermal equilibrium, which cost a lot of energy per magnetization reduction.The finding is corroborated by inspecting the simulated spin configuration at short time delays (see Appendix B), which exhibits disorder on small length scales, i.e. also between spins that are close to each other.For the fluences reached in our experiments, the spin system thermalizes within the first picosecond after laser excitation.In contrast, for Ni, we observed a prolonged nonthermal behavior for higher fluences.In contrast, the dashed curves show the additional spin energy content of a hypothetical, thermalized spin system with the magnetization dynamics from the ASD simulations, which was calculated using the equilibrium relationship between spin energy and magnetization.Differences between the solid and the dashed curves indicate a nonthermal spin system.Note that the very small differences that persist on timescales larger than 3 ps are numerical artifacts.They could stem from the finite time steps in the nonequilibrium simulations or the larger α employed in the simulations of quasistatic heating.The insets show the equilibrium relationships between spin energy and magnetization.
These differences between Fe, Co, and Ni are caused by their different Curie temperatures.Ni has a Curie temperature of only 631 K, while the Curie temperatures of Fe and Co are 1044 K and 1390 K, respectively [51].As a consequence, for the same absorbed energy density, Ni demagnetizes much more than Fe or Co [43].For stronger demagnetization of Fe or Co, a prolonged nonthermal behavior is observed as well, which is shown exemplarily for Fe in Appendix B. The prolonged nonthermal behavior is found to be caused by domain formation during the remagnetization process, in agreement with previous results by Kazantseva et al. [15].The ASD simulation results thus suggest that in particular for strong demagnetization, a thermal description of the spin system is not adequate to describe nonequilibrium dynamics of ferromagnets.
In summary, according to the ASD simulations, two key aspects of ultrafast energy flow dynamics in ferromagnets are apparent: the slowdown of the lattice response caused by energy flow into and out of the spin system, and the nonthermal behavior of the spin system on short timescales.The first effect can be qualitatively reproduced with a simple 3TM as well (depending on the coupling constants, which are typically fit parameters).However, the second key aspect, the transient nonthermal excitation of the spin system, cannot be modeled with a 3TM.Therefore, we expect worse quantitative agreement of a 3TM with the ultrafast dynamics of all subsystems.
Besides the ASD simulations presented here and the 3TM, another model describing the coupled energy flow dynamics of electrons, spins, and the lattice is the spin-lattice-electron dynamics (SLED) model by Ma et al. [52].
In contrast with our simulations, the SLED model employs the effective spin temperature [53] to calculate the energy flow dynamics.In order to investigate differences between the two descriptions, a comparison analogous to Fig. 6 can be made for the effective spin temperature.In our simulations, we find that on short timescales the energy content calculated using the effective spin temperature can differ from the energy of the spin system calculated using the spin Hamiltonian (see Eq. 5), in particular for iron and nickel and for high fluences.Such deviations from equilibrium relationships indicate a nonthermal behavior of the spin system, as discussed earlier.An additional difference between our model and the SLED model is that the latter includes direct spin-lattice coupling.For these two reasons, we expect qualitatively similar results for the energy flow dynamics but quantitative differences between these two models.

V. SUMMARY AND CONCLUSIONS
In this work, we investigated the ultrafast lattice dynamics of ferromagnetic Co and Fe using femtosecond electron diffraction.To model the intrinsic energy flow between electronic, magnetic, and lattice degrees of freedom, we combined spin-resolved DFT calculations of the electron-phonon coupling with energy-conserving ASD simulations.We found that for both Co and Fe, the ultrafast spin dynamics have a profound impact on the lattice dynamics, slowing down the lattice heating due to energy transfer into and out of magnetic degrees of freedom.These findings generalize our previous results for Ni [17], highlighting the prominent role of the spin system in the energy flow dynamics of all three elemental 3d ferromagnets.
For a full description of the laser-induced dynamics, it is thus essential to take energy flow into and out of the spin system into account.This is achieved with energy-conserving ASD simulations, which simulate the spin dynamics while also accounting for the intrinsic energy flow between electrons, spins, and the lattice.For low and moderate fluences, the ASD simulations yielded excellent agreement with the measured lattice dynamics, as well as a good description of the magnetization dynamics for both Co and Fe.They are thus an important step toward a model for ultrafast demagnetization that is consistent with the responses of electronic, magnetic, and lattice degrees of freedom.
Furthermore, we found that the ASD simulations predict a nonthermal spin system for both Co and Fe on short timescales after laser excitation.For high fluences, the nonthermal state of the spin system can last for several picoseconds, suggesting that particularly for strong excitations, a thermal description of the laserinduced spin dynamics is not sufficient.
Our findings are also of relevance for other demagnetization models since they enable the comparison with the experimental lattice dynamics for all three elemental 3d ferromagnets and highlight the importance of a consistent description of the energy flow dynamics.The combination of experiment and theory presented in this work can also be applied to gain insight into the ultrafast energy flow dynamics in other technologically relevant magnetic materials, e.g.magnetic oxides and layered van der Waals materials.In addition, the incorporation of the energy exchange of the spin system in the ASD simulations may prove to be invaluable for the explanation of the behavior of more complex materials and heterostructures in the future.
The femtosecond electron diffraction data and the ASD simulation results presented here are available on a data repository [54,55].The code to calculate the electronphonon coupling parameter and the heat capacities from DFT results is also available, along with results for iron, cobalt, and nickel [56].
The optimized norm-conserving Vanderbilt pseudopotentials were generated using the method of Ref. 62 and are of generalized-gradient approximation Perdew-Burke-Ernzerhof type [63].Sixteen electrons were treated explicitly for Fe, and 17 electrons were explicitly taken into account for Co.The plane-wave expansion of the electronic wave function had a cutoff of 40 Ha for Fe and 50 Ha for Co; 22 electronic bands were calculated for Co and 15 for Fe.These bands were calculated with Fermi occupation featuring a smearing of 0.001 Ha.An unshifted k-point grid of 32 × 32 × 32 points was used for both elements.The lattice constant for body-centered cubic (bcc) Fe was set to 2.756 Å, which was obtained by relaxing the structure.For hcp Co, we used the experimental lattice constants a = 2.5071 Å and c = 4.0695 Å [64].To obtain the electron-phonon coupling G ep , the spin-resolved electron-phonon matrix elements were computed as described in Ref. 34 for a 8 × 8 × 8 grid of q points.From the results, we extracted the Eliashberg functions for majority and minority electrons.The electron-phonon couplings and the electronic heat capacities were then calculated as described in Ref. 35.Following Ref. 17, we take chemical potential shifts into account and assume particle conservation within each spin type.Band shifts according to the Stoner model are not considered, since our description of magnetization dynamics with ASD simulations is based on the Heisenberg model.
The results for the electron DOS, the Eliashberg functions, and the electron-phonon couplings are presented in Fig. 7.The magnetic moments calculated from the spin-resolved electron DOS, 1.95 µ B per atom for Co and 2.40 µ B per atom for Fe, are larger than the experimental results of 1.72 µ B and 2.22 µ B per atom [65].
Based on the phonon DOS, we also calculated the MSDs as functions of temperature, as described in Ref. 39.To increase the accuracy of the calculation, we replaced the phonon DOS in the region below 5 meV by a fit with the function f (x) = ax 2 + bx 3 .This ensures that the dominating term for very small phonon wavevectors is quadratic, which corresponds to the correct long-wavelength limit.The results were used to convert transient MSD changes to lattice temperatures [see Figs.
with S i representing a classical, normalized spin vector at site i.Each spin couples to its neighboring spin vectors S j via the coupling constant J ij .The second term of the Hamiltonian (Eq. 6) describes the on-site anisotropy with an easy axis along the z axis and constant anisotropy energy.All parameters are material-dependent and listed in Table I.Except α, they are based on Ref. 68.The  simulations are performed on a simple cubic (sc) lattice.This has no significant effect on the energy flow or magnetization dynamics, which was verified directly by comparing simulations of bcc and sc Fe.The reason is that for a given Curie temperature, the different number of neighbors in each case is compensated by a different value of J.Note that in contrast with the samples employed in the diffraction experiments, the ground state in the ASD simulations is a single-domain state.Due to the typical time and energy scales of domain wall dynamics, we don't expect a significant influence of the domain structure on the intrinsic energy flow dynamics studied here.By solving the s-LLG equation numerically, the dynamics of the system are calculated [69].Here γ = 1.76 • 10 11 1 Ts refers to the gyromagnetic ratio and H i describes the effective field derived via H i = − ∂H ∂Si .The material-dependent and phenomenological damping parameter α determines the coupling strength of the spin system to the electron system and thus the energy transfer rate between the two heat baths.In order to simulate the effects of finite temperatures, a Langevin thermostat is included by adding a field-like stochastic term ζ i to the effective field The added noise term has white noise properties [70]: In order to better reproduce experimentally measured equilibrium properties such as magnetization and heat capacity, we make use of a rescaled temperature model, which utilizes a slightly modified electron temperature T sim for the noise generation.
Figure 8 presents simulation results for equilibrium properties and a comparison to literature values.In addition to the equilibrium properties, the rescaled temperature model also yields a better description of nonequilibrium dynamics [71].Further details are available in Refs.17 and 71.
A major advantage of ASD simulations is that they   9. Details about the spin configuration at selected time delays from an ASD simulation of Fe at a high excitation density.The upper panel shows the additional spin energy as a function of pump-probe delay (solid curve) for an absorbed energy density of 1930 J/cm 3 .The dashed curve corresponds to the energy content of a thermalized spin system with the magnetization dynamics from the ASD simulations (shown directly below), analogous to Fig. 6.Note that for illustration purposes, we used a higher damping parameter of α = 0.1 for the simulations shown here, and a spin system consisting of 100 × 100 × 100 spins.Due to the relatively high fluence, the spin system exhibits a nonthermal behavior also on few-picosecond timescales.The lower part of the figure shows the spin configurations at different pump-probe delays.The spin components Sx and Sz are displayed, normalized to 1. Sy behaves analogously to Sz for symmetry reasons.Here only the surface of the cube is visible.The inside of the cube displays an analogous behavior to the surface.
are not limited to a thermal description of the spin system, since the spins are simulated directly.Figure 9 shows the evolution of the spin energy content and a direct visualization of the simulated spin dynamics for a relatively high fluence of 1930 J/cm 3 .In addition, for illustration purposes, a simulation with a higher damping of α = 0.1 is shown.A higher damping leads to a larger disorder of the spin system directly after excitation, however, the qualitative behavior displayed in Fig. 9 is also observed for lower values of the damping parameter at high fluences.
The comparison of the simulated spin energy content (solid curve) to the energy content of a thermalized spin system with the simulated magnetization dynamics (dashed curve) reveals that the spin system remains in a nonthermal state for several picoseconds.At short time delays, the nonthermal state is characterized by a relatively large spin energy content compared with the demagnetization amplitude.The behavior reverses on longer timescales.
Further insights on these nonthermal states can be gained from the visualization of the spin dynamics.The instantaneous spin configuration is illustrated exemplarily for several delays after excitation.
During and directly after excitation, e.g., at t 1 = 0 ps, there is significant short-range disorder in the spin system, i.e., significant misalignment between neighboring spins.According to the Heisenberg Hamiltonian, this comes with a significant energy cost and thus leads to the relatively large energy content of the spin system.On longer timescales, the magnetization recovers.However, for high fluences/strong demagnetization, domains form.This is already visible at t 2 = 0.75 ps.There are areas with significant magnetization in which the spins point predominantly in the x or y direction (note that in Fig. 9, each spin is normalized such that S 2 x +S 2 y +S 2 z = 1).This behavior is similar to spin simulation results reported in Ref. 15. Within a domain, spins are parallel.Therefore, the energy cost of this spin configuration is relatively low.Nevertheless, due to the different directions of the magnetization, the global magnetization is reduced.In the beginning, the domains are relatively small.As time progresses, the domains become larger (see t 3 = 1.5 ps and t 4 = 4.5 ps) and eventually disappear (see t 5 = 9 ps) as the spin system approaches thermal equilibrium.Note that for low fluences, domain formation as illustrated in Fig. 9 does not occur, since it requires significant initial disordering of the spin system.The initial nonthermal disorder of the spin system, visualized in Fig. 9 for t 1 and characterized by a large number of high-energy spin excitations, occurs for all fluences (see also Fig. 6).

FIG. 1 .
FIG. 1.The femtosecond electron diffraction experiment.(a) Schematic illustration of the measurement.The samples are thin, freestanding films on TEM grids, which are excited by ultrashort laser pulses.The lattice response is probed using ultrashort electron pulses, which diffract off the sample.Diffraction patterns are recorded in transmission.In the case of polycrystalline samples such as the samples studied in this work, the diffraction patterns consist of rings.A diffraction pattern of our Co sample is shown.(b) Diffraction pattern of Co and time-resolved changes.The upper part shows the azimuthally averaged diffraction pattern (radial profile, RP) of Co.Here the background-subtracted pattern is shown for illustrational purposes; however, note that in the analysis of the diffraction patterns, fits are performed to background and rings simultaneously.The lower part shows the differences of the RPs compared with the RP before laser excitation for several pump-probe delays (solid curves).The dashed black lines show the fit results of the global fitting routine described in Sec.II and in detail elsewhere[37].(c) Same as (b) but for Fe.

FIG. 2 .
FIG. 2. Experimental results for the lattice dynamics and single exponential fits.(a) Evolution of the atomic MSD and corresponding lattice temperature for Co.(b) MSD evolution and corresponding lattice temperature for Fe.The solid lines in (a) and (b) are the results of fits to the experimental data with a single exponential function, convolved with a Gaussian (250 fs FWHM) to account for the time resolution of the experiment.The excitation wavelength was 2300 nm.(c) Fit results for the time constant of the single exponential function for different excitation densities, yielding different final lattice temperatures.

FIG. 3 .
FIG. 3.Heat capacities, electron-phonon coupling parameters, and schematic illustrations of the employed energy flow models.(a) Electronic (dashed curves) and lattice (solid curves) heat capacities, and (b) electron-phonon coupling parameters Gep as a function of electron temperature calculated from spin-resolved DFT results.

3 FIG. 4 .
FIG. 4. Experimentally measured lattice dynamics and model predictions.(a)-(d) Lattice temperature as a function of pumpprobe delay in Co for different absorbed energy densities.The experimental data are shown as black circles.The results of the TTM are shown as dashed curves and the results of the ASD simulations are shown as solid curves.The TTM and ASD results were convolved with a Gaussian with a FWHM of √ 250 2 − 80 2 fs ≈ 237 fs, which accounts for the temporal broadening induced by the probe pulse (the effect of the pump pulse width is already included in the model itself).The gray-shaded areas represent the standard errors of the experimental data, obtained from the fitting routine described in Section II.The displayed energy densities correspond to the absorbed energy densities in the ASD simulations.(e) Magnetization dynamics of Co predicted by the ASD simulations.(f)-(i) Experimental results for the lattice temperature in Fe alongside results of the TTM and the ASD simulations.(j) Magnetization dynamics of Fe predicted by the ASD simulations.

FIG. 5 .
FIG. 5. Intrinsic energy flow between electrons, spins, and the lattice for (a) Co and (b) Fe.Here results with the same excitation density as in Fig. 4(b) (Co) and Fig. 4(g) (Fe) are presented.The additional energy density ∆E in the system is displayed.After laser excitation, the total energy (black) stays constant and energy is redistributed between electronic (blue), magnetic (green), and lattice (red) degrees of freedom.

3 ]FIG. 6 .
FIG.6.Nonthermal spin dynamics for (a) Co and (b) Fe.The solid curves show the additional energy content of the spin system ∆Es as a function of pump-probe delay for the same fluences as in Fig.4.In contrast, the dashed curves show the additional spin energy content of a hypothetical, thermalized spin system with the magnetization dynamics from the ASD simulations, which was calculated using the equilibrium relationship between spin energy and magnetization.Differences between the solid and the dashed curves indicate a nonthermal spin system.Note that the very small differences that persist on timescales larger than 3 ps are numerical artifacts.They could stem from the finite time steps in the nonequilibrium simulations or the larger α employed in the simulations of quasistatic heating.The insets show the equilibrium relationships between spin energy and magnetization.
APPENDIX B: ASD SIMULATIONSAtomistic spin dynamics simulations use a classical Heisenberg spin model:

FIG. 7 .
FIG. 7. Results of the spin-resolved DFT calculations.(a) Spin-split electron DOS of Co.The position of the Fermi level is shown as a gray line.The majority (majo) DOS is shown in dashed blue and the minority (mino) DOS is shown in red.Note that hcp Co has two atoms per primitive unit cell.(b) Spin-split electron DOS of Fe.(c) Majority and minority Eliashberg functions α 2 F for Co (blue) and Fe (green).The dashed curves correspond to the majority Eliashberg functions and the solid curves represent the minority Eliashberg functions.(d) Majority and minority electron-phonon coupling parameter Gep for Co and Fe.

FIG. 8 .
FIG. 8. ASD simulation results for equilibrium relationships and comparison to literature results.(a) Spin heat capacity of Co.The ASD simulation result (ASD sim.) is shown as a solid blue curve and the experimental result (exp.) is shown as a dashed black curve.The experimental result was obtained based on measurements of the total heat capacity from Ref. 66, which were dilation-corrected using the expansion coefficients from Ref. 67.To obtain the spin heat capacity, the DFT results for the electronic and lattice contributions were subtracted.(b) Same as (a) but for Fe.(c) Magnetization as a function of temperature for Co.The solid curve shows the ASD simulation result.The dashed black curve is a literature result from Ref. 51.The magnetization is normalized to its value at 0 K.(d) Same as (c) but for Fe.

TABLE I .
ASD simulation parameters for Co and Fe.