# Coupled atomistic spin-lattice simulations of ultrafast demagnetization in 3d ferromagnets

M. Pankratova, I P Miranda, D. Thonig, M. Pereiro, E. Sjöqvist, A. Delin, P.

Scheid, O. Eriksson, A. Bergman

► To cite this version:

M. Pankratova, I P Miranda, D. Thonig, M. Pereiro, E. Sjöqvist, et al. Coupled atomistic spin-lattice simulations of ultrafast demagnetization in 3d ferromagnets. Scientific Reports, 2024, 14 (1), pp.8138. 10.1038/s41598-024-58662-y. http://dx.doi.org/10.1038/s41598-024-58662-y

HAL Id: hal-05040351

https://hal.science/hal-05040351v1

Submitted on 19 Apr 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_2.png" /></figure>

OPEN

# Coupled atomistic spin-lattice simulations of ultrafast demagnetization in 3d ferromagnets

M. Pankratova¹,⁹,✉, I. P. Miranda¹,⁹,✉, D. Thonig¹,², M. Pereiro¹, E. Sjöqvist¹, A. Delin³,⁴,⁵, P. Scheid⁶,⁷, O. Eriksson¹,⁸ & A. Bergman¹

Despite decades of research, the role of the lattice and its coupling to the magnetisation during ultrafast demagnetisation processes is still not fully understood. Here we report on studies of both explicit and implicit lattice effects on laser induced ultrafast demagnetisation of bcc Fe and fcc Co. We do this using atomistic spin- and lattice dynamics simulations following a heat-conserving three-temperature model. We show that this type of Langevin-based simulation is able to reproduce observed trends of the ultrafast magnetization dynamics of fcc Co and bcc Fe. The parameters used in our models are all obtained from electronic structure theory, with the exception of the lattice dynamics damping term, where a range of parameters were investigated. It was found that while the explicit spin-lattice coupling in the studied systems does not impact the demagnetisation process notably, the lattice damping has a large influence on the details of the magnetization dynamics. The dynamics of Fe and Co following the absorption of a femtosecond laser pulse are compared with previous results for Ni and similarities and differences in the materials' behavior are analysed. For all elements investigated so far with this model, we obtain a linear relationship between the value of the maximally demagnetized state and the fluence of the laser pulse, which is in agreement with experiments. Moreover, we demonstrate that the demagnetization amplitude is largest for Ni and smallest for Co. This holds over a wide range of the reported electron-phonon couplings, and this demagnetization trend is in agreement with recent experiments.

Ultrafast demagnetization was discovered by Beaurepaire and coauthors in 1996¹. They observed demagnetization in a nickel film on picosecond timescales following the absorption of a femtosecond laser pulse. From the point of view of applications, ultrafast demagnetization is an important process in all-optical magnetization switching as well as for new applications in magnetic data storage and spintronics². In the same pioneering work of Beuarepaire¹, these experimental observations were interpreted using a three-temperature model (3TM)¹, which assumes three thermalized reservoirs, in particular, spin, lattice, and electron reservoirs, that can exchange energy through coupling parameters (electron-phonon, $G_{ep}$, electron-spin, $G_{es}$, and spin-lattice, $G_{sl}$). The 3TM is often used to interpret ultrafast magnetization dynamics processes¹,³,⁴. Recently, many other models have been proposed²,³,⁵,⁶ to describe possible mechanisms of ultrafast demagnetization, such as the importance of spin-dependent transport of laser-excited electrons⁷, the optical intersite spin transfer effect (OISTR) that was considered in Ref.⁸, and the Elliott-Yafet electron-phonon spin-flip scattering, studied in Ref.⁹. Recently, the competition between several ultrafast processes – same-site spin transfer, intersite spin transfer, and ultrafast spin

¹Department of Physics and Astronomy, Uppsala University, Box 516, 75120 Uppsala, Sweden. ²School of Science and Technology, Örebro University, 70182 Örebro, Sweden. ³Department of Applied Physics, School of Engineering Sciences, KTH Royal Institute of Technology, AlbaNova University Center, 10691 Stockholm, Sweden. ⁴Swedish e-Science Research Center (SeRC), KTH Royal Institute of Technology, 10044 Stockholm, Sweden. ⁵Wallenberg Initiative Materials Science for Sustainability (WISE), KTH Royal Institute of Technology, 10044 Stockholm, Sweden. ⁶LPCT, CNRS, UMR 7019, BP 70239, Université de Lorraine, 54506 Vandoeuvre-lés-Nancy Cedex, France. ⁷IJL, CNRS, UMR 7198, BP 70239, Université de Lorraine, 54000 Nancy Cedex, France. ⁸Wallenberg Initiative Materials Science for Sustainability, Uppsala University, 75121 Uppsala, Sweden. ⁹These authors contributed equally: M. Pankratova and I. P. Miranda. ✉email: maryna.pankratova@physics.uu.se; ivan.miranda@alumni.usp.br
flips—was studied in some detail¹⁰. In¹¹ a dynamic spin-lattice-electron model was proposed. Using this model, the authors calculated laser-induced demagnetization of iron thin film and obtained very accurate agreement with the experimental observations. One of the important outcomes of Ref.¹¹ is the establishment of the relations between the dissipative parameters entering the Langevin equations for the lattice and spin degrees of freedom and the subsystems coupling coefficients of 3TM.

In line with those developments, recently Zahn et al.¹²,¹³ proposed an energy-conserving model based on atomistic spin dynamics simulations. Using this model they studied ultrafast demagnetization of nickel¹², iron and cobalt¹³, showing a cohesive microscopic picture of the laser-induced dynamics in those 3d ferromagnets. Alternatively, our previous work⁴ proposed a heat-conserving three temperature model (HC3TM), where the heat distribution between spin and lattice subsystems is *measured* during the coupled atomistic spin dynamics simulations, in contrast with calculation of temperature in advance. Based on the HC3TM model, Ref.⁴ demonstrated that the fast interplay between spins and electrons allows for a shorter demagnetization time during the laser pulse, when compared to the regular 3TM with constant coefficients, i.e., HC3TM appears to be more consistent with experimental observations for fcc Ni compared to 3TM¹,¹⁴,¹⁵. In addition, as argued in Ref.⁴, the HC3TM has the advantage to rely to a much smaller degree on model parameters that can be difficult to estimate, such as, for example, *G*<sub>es</sub>, *G*<sub>ep</sub>, and *G*<sub>sl</sub> of 3TM. In fact several authors⁴,¹²,¹³ report on drastically different 3TM coupling parameters for the here studied systems, and when there is a scatter in such data the reliability of the model becomes questionable.

The impact of laser fluence of ultrafast magnetization dynamics was extensively studied both experimentally and theoretically¹⁶⁻²⁰. In particular, it was observed in¹⁹ for nickel and also for other materials in Refs.²¹,²² that the amplitude of the maximally demagnetized state depends linearly on laser fluence. New experimental results constantly bring new insights and adds to the discussion of existing demagnetization models. As a recent example, Scheid et al.²⁰ while experimentally studying maximally demagnetized state of iron, nickel, and cobalt, observed that nickel is the easiest element to demagnetize, however, demagnetization amplitudes for iron and cobalt were almost similar in spite of quite different electronic structures, saturation moments and Curie temperatures. Notably, a linear decrease of the magnetization with the laser fluence contradicts Bloch's law, stating that the magnetization is proportional to T<sub>m</sub><sup>3/2</sup>, where T<sub>m</sub> is the magnon temperature. A model to explain the experimental observations was proposed based on the assumption that the linear dependence of demagnetization on fluence is driven by the increase in temperature, the electron-phonon coupling, the electron-magnon scattering, and a reduction of the interatomic exchange.

In the present work, we further study ultrafast demagnetization of 3d ferromagnets with a focus on the role
of the lattice dynamics, spin-lattice interaction, and the performance of the recently suggested HC3TM with the
strong coupling to ab-initio electronic structure theory. We investigate in full detail the parameter space within
the HC3TM model, e.g. by analyzing the impact of the lattice damping, and microscopic spin-lattice coupling
during the ultrafast demagnetization. In addition, we apply HC3TM to study the experimentally observed (linear)
dependence of the demagnetization amplitude on laser fluence of 3d ferromagnets (bcc iron, fcc cobalt, and fcc
nickel). The choice of comparing these three elements is natural given the fact that they are the only ferromagnetic
transition metal elements, with distinctly different magnetic moments, ordering temperature and how they are
understood in terms of being strong or weak ferromagnets.

Results

Heat-conserving three-temperature model

The heat-conserving three-temperature model was proposed in Ref.⁴ for the investigation of ultrafast demagnetization of fcc Ni. Similar to Beurepaire's 3TM (see Fig. 1a,b), it is assumed that lattice and spin systems are connected to an electronic heat bath. If the electronic reservoir is not considered infinite, the electronic temperature

<figure><img src="image_3.png" /><figcaption>**Figure 1.** The main idea behind the HC3TM (a) and 3TM (b) models. The text along the arrows indicates the parameters that drive the heat transfer in each model. For the parameters of the 3TM, see Ref.¹.</figcaption></figure>
should be decreased as the spin and lattice temperatures increase. This can be modeled by considering the conservation of heat in the system, i.e., the heat that is transferred into the spin and lattice systems is deducted from the electronic system (see Supplementary materials). As mentioned in the Introduction, in the original 3TM formulation the temperature of each subsystem is governed by the coupling coefficients G<sub>ep</sub>, G<sub>es</sub>, and G<sub>sl</sub><sup>1</sup>.

In the HC3TM formulation, the heat absorbed in the spin and lattice subsystems is simulated on the fly using atomistic spin-lattice dynamics. Thus, the coupling between the subsystems is best described by the atomistic quantities v (lattice damping), α (Gilbert damping), and the explicit spin-lattice coupling Γ<sup>αβμ</sup><sub>ijk</sub> (see Eq. (8) in Section “Methods”) in the context given by Hellsvik et al.<sup>23</sup>, instead of the microscopic coupling parameters G<sub>ep</sub>, G<sub>es</sub>, G<sub>sl</sub> of 3TM (see Fig. 1a,b). In that sense, HC3TM allows for a consistent microscopic (and atomistic) description of the model parameters, in a similar fashion as in Ref.<sup>11</sup>.

In the current formulation, all parameters of this theory (α<sup>24</sup>, ν<sup>11,25</sup>, and Γ<sup>αβμ</sup><sub>ijk</sub><sup>23</sup>) can be calculated using, for example, density functional theory, although in the present work we varied the lattice damping for a range of realistic values, as opposed to performing an explicit calculation of it. The relation between the lattice damping and G<sub>ep</sub> has been established in<sup>11</sup> using a spin-lattice-electron three-temperature model. It was shown that the lattice damping is directly proportional to G<sub>ep</sub>, which was previously also reported by Duffy et al. in Refs.<sup>26,27</sup>.

In contrast to 3TM, where the temperatures can be calculated in advance¹, within the HC3TM, the temperature of the different subsystems is calculated on-the-fly—i.e., at every time step of the spin-lattice dynamics simulation—by considering the following relation (see Supplementary materials for motivation):

$$ \Delta T_e(t) = - \frac{C_l(T_l)}{C_e(T_e)} \Delta T_l(t) - \frac{C_s(T_s)}{C_e(T_e)} \Delta T_s(t) + \frac{W(t)}{C_e(T_e)}. \quad (1) $$

Here, Δ represents the difference with regard to the temperature, T<sub>l</sub> is the lattice temperature, calculated from the average kinetic energy of the lattice vibrations; ⟨E<sup>kin</sup><sub>l</sub>⟩/k<sub>B</sub>. The spin temperature, T<sub>s</sub>, is calculated according to<sup>28</sup>, i.e., using T<sub>s</sub> = (Σ<sub>i</sub> |m̂<sub>i</sub>×B<sub>i</sub>|<sup>2</sup>) / (2k<sub>B</sub>(Σ<sub>i</sub> |m̂<sub>i</sub>·B<sub>i</sub>|<sup>2</sup>)), where m̂<sub>i</sub> is the normalized local spin moment, B<sub>i</sub> an effective exchange field (see “Methods” for calculation details). The spin temperature, T<sub>s</sub>, used in our simulations is defined for an equilibrium situation where the energy flow out of the system is zero. During the ultrafast processes studied in this work, the energy flow is non-zero which raises the question of a reliable definition, and measurement of out-of-equilibrium temperatures<sup>29</sup>. However, in our simulations, we have access to the rate of change of the energy of the spin system, and comparing this rate with the terms that end up defining the spin temperature (Eqs. 15 and 16 in Ref.<sup>28</sup>) we find that it is actually not large enough to notably change the numerical value of the calculated spin temperature. Furthermore, we know from earlier tests that if we translate electron-phonon and electron-magnon couplings from previous 3TM studies to our HC3TM model, we can essentially reproduce the same temperature profiles and demagnetization curves. Together, these two observations support the validity of the use of the equilibrium spin temperature definition (please see Supplementary materials for more information). C<sub>e</sub>(T<sub>e</sub>), C<sub>s</sub>(T<sub>s</sub>), and C<sub>l</sub>(T<sub>l</sub>) in Eq. (1) denote the temperature dependent specific heats of the electronic, spin, and lattice subsystems correspondingly (see Supplementary materials for more information on the heat capacities of the spin, lattice, and electron subsystems and their units). W(t) represents the impact of a laser pulse, that increases the electronic temperature of Eq. (1), which is proportional to the power source term, modeled as a Gaussian function.

More details concerning the HC3TM and its comparison with 3TM can be found in Ref.<sup>4</sup> and Supplementary materials.

## HC3TM for cobalt and iron

In this section we present results for spin-lattice dynamics simulations of fcc Co and bcc Fe. In both cases, we run simulations for a 60 × 60 × 60 repetition of the fcc or bcc unit cell with periodic boundary conditions. The simulations are performed for N<sub>t</sub> = 1 × 10<sup>6</sup> time steps of dt = 10<sup>-16</sup>s. Other parameters of our simulations are presented in the Table I in Supplementary Materials.

We start by describing results from ultrafast demagnetization of fcc Co. First, we present temperature profiles (Fig. 2a) for spin, lattice, and electron subsystems, for a laser fluence of 45 J/m<sup>2</sup>, together with the corresponding magnetization dynamics. Other parameters relevant for these simulations are listed in Table I in Supplementary materials. It can be seen from the figure that similar to experimental observations<sup>30</sup> and theoretical studies<sup>13</sup> demagnetization happens on subpicosecond timescales, followed by a remagnetization that has a faster recovery after ~1 ps and a slower remagnetization after that. We obtain values of the position of the magnetization minima that are close to experimental values, and we also obtain demagnetization/remagnetization timescales that are consistent with experimental observations. It is important to note that in Ref.<sup>13</sup>, using the parametrized three-temperature model, features of the magnetization dynamics similar to the data of Fig. 2b were also observed, however for an adjusted Gilbert damping value of 0.1, which is somewhat unrealistically large<sup>24</sup>. In contrast, in the calculations presented here we use a damping 0.0024, close to experimental value (see Table I in Supplementary materials).

Similar coupled spin-lattice simulations using the HC3TM, were performed for bcc Fe, and shown in Fig. 2c,d. Just like for fcc Co we first present in Fig. 2d magnetization dynamics after the absorption of a laser pulse and the corresponding spin, lattice, and electronic temperatures, see Fig. 2c. Using the HC3TM model we obtain realistic demagnetization/remagnetization times in comparison with experimental studies and in this case even the maximally reduced magnetization is quite close to the experimentally observed value<sup>9</sup>, for the same pulse fluence (M/M<sub>0</sub> ~ 0.9 in our calculations).
<figure><img src="image_4.png" /><figcaption>**Figure 2.** Spin, lattice and electron temperatures of fcc Co (a) and bcc Fe (c) and magnetization dynamics of fcc Co (b), and bcc Fe (d) obtained in the HC3TM with fluence of 45 J/m² and 15 J/m² correspondingly.</figcaption></figure>

## Results from coupled spin-lattice simulations

To study the impact of direct spin-lattice coupling on the magnetization dynamics (see Eq. (8) in "Methods") we consider explicitly the spin-lattice coupling, where we obtain the elements of the tensor, $\Gamma_{ijk}$ in the spin-lattice coupling, from calculations based on density functional theory (e.g., as described in Ref. 23). Calculated values of these parameters were used to study the impact of direct spin-lattice coupling for bcc Fe, and fcc Co. This explicit, calculated spin-lattice coupling allows for a direct exchange of heat between the spin and lattice subsystems, outside of the channel provided by the HC3TM. Note that the HC3TM can be applied whether or not the direct spin-lattice term (see Eq. (8) in "Methods") is included or neglected. The results that follow below were obtained from the HC3TM with spin-lattice coupling included explicitly in the simulations and compared with the case without spin-lattice coupling.

In Fig. 3 we show the demagnetization profile for bcc Fe from coupled spin-lattice dynamics simulations in combination with the HC3TM. The simulations are, as noted, done with and without the spin-lattice coupling term (given by expression Eq. (8) in "Methods"). It can be seen from Fig. 3 that spin-lattice coupling impacts the magnetization dynamics weakly for bcc Fe. If any effect can be observed, it would be a slightly slower remagnetization rate when the spin-lattice coupling is included. The difference compared to the scenario without the spin-lattice coupling is however almost of the same order of magnitude as the numerical noise of the simulations. When artificially increased by three times, the spin-lattice coupling has a more prominent effect, however, it remains small, as shown in Supplementary materials. We note that a simulation which completely ignores this term still has ability to transfer heat from the spin system to the lattice, via Eq. (1). A possible explanation for the minor influence of the spin-lattice couplings can be sketched as follows: even at the maximum spin temperature in the simulations ($T_s \lesssim 600$ K for Fe), the spin system is far away from the Curie temperature ($T_C$), i.e., with a general strong coupling between spin moments. Thus, the small average displacements (in the order of $\sim 10^{-2}a$)

<figure><img src="image_5.png" /><figcaption>**Figure 3.** Impact of spin-lattice coupling on magnetization dynamics of bcc Fe. The pulse fluence is 20 J/m².</figcaption></figure>
and their distribution throughout the sample result in the summed effect of spin-lattice coupling being only a minor perturbation on the total exchange field. In addition, the spin-lattice couplings considered in this study does not allow the transfer of angular momentum between the spins and the lattice²³ but since the magnitudes of the calculated spin-lattice couplings are found to be small we do not expect that allowing for a full angular momentum transfer would result a significant deviation from our results.

Also, for bcc Fe the resulting profile is in rather good agreement with experimental observations. For materials with larger spin-lattice coupling (such as CrI₃³¹) it can be expected that the importance of spin-lattice coupling increases, and our simulations show that for bcc Fe the spin-lattice coupling becomes more important for higher pulse fluences (data not shown).

## Impact of lattice damping

While the impact of spin-lattice coupling, via Eq. (8), on the magnetization dynamics is only marginal for the elemental 3d ferromagnets, the impact of the dynamical properties of the lattice itself, and in particular, the lattice damping, is significant. This is in line with Refs.¹²,¹³,³²,³³ (among others), which highlight the importance of the lattice degree of freedom in the ultrafast demagnetization process. We illustrate this for fcc Co, in Fig. 4. Here it is seen that a decrease of the lattice damping (ν in Eq. (4)) leads to an increase of the demagnetization; both the magnetization minimum as well as the reduced magnetization (M/M₀) are influenced by the lattice damping parameter. This can be better understood from the fact that the coupled spin and lattice system allows for heat to be dissipated from both reservoirs, when the lattice (and spin) damping parameter is larger. This causes the temperature profile to reach lower values and to equilibrate quicker both for the spin- and lattice system, as shown in Fig. 4, with a resulting stronger impact on the magnetization profile for lower values of the lattice damping. This is detailed clearly when comparing the spin temperatures (Fig. 4) for various values of lattice damping. Lower lattice damping values lead to higher spin temperature and therefore, to larger magnetization drop and a longer remagnetization time (see the dynamics of lattice and electronic temperatures for this case in the Supplementary materials). One of the important outcomes from this observation is that it is essential to take the lattice dynamics properly into account while studying ultrafast demagnetization. If the lattice is not considered at all, then one needs to note, while comparing simulations with experimental data, that this will lead to an overestimation of the magnetization drop and to longer remagnetization times, for a given laser pulse.

Both lattice and Gilbert damping parameters can change as a function of the electronic temperature, and this variation can in principle be captured by *ab-initio* calculations (or, at least, theoretically by using ground-state *ab-initio* quantities), as fair approximations. Let us roughly analyze the magnitudes of these changes. Refs.¹²,¹³ show that $G_{ep}$ increases ~ 20% (Ni), ~ 30% (Co), and ~ 40% (Fe) between 300 K $\le T_e \le T_e^{max}$ (and so do $\nu$). In turn, $\alpha(T_e)$ changes ~ 400% in the same temperature interval for Co³⁴ (we are assuming that analogous change also happens to Co fcc). The increase of lattice damping will make the heat flow from the electronic system to the lattice system more intensively, so that less heat is redirected to the spin system. This means that we can expect an effect much smaller than the one shown in Fig. 4 coming from $\nu$, and a more intense effect coming from $\alpha(T_e)$ in Co, in the direction of increasing the demagnetization rate (but counterbalanced by $\nu$, as discussed). Therefore, in the investigated systems and fluences, our simple analysis suggests that assuming $\nu$ and $\alpha$ constants is a decent approximation for iron and nickel, however with some changes expected for cobalt.

## Comparison of magnetization dynamics: bcc Fe, fcc Co, fcc Ni

In this section we discuss magnetization dynamics of bcc Fe and fcc Co. We compare the results from the two elements to previous results reported for fcc Ni, which also has been studied using the HC3TM⁴. The results are

<figure><img src="image_6.png" /><figcaption>**Figure 4.** Spin temperature (a) and magnetization dynamics (b) of fcc Co for various lattice damping values used in HC3TM ($\nu$ in Eq. (4)), with temperature-dependent heat capacities (see Ref.⁴ for details). The pulse fluence is 45 J/m².</figcaption></figure>
summarized in Fig. 5. We begin with comparing ultrafast magnetization dynamics of fcc Co, fcc Ni, and bcc Fe, for the same laser pulse fluence, presented in Fig. 5a. It can be seen that in this case, the demagnetization for fcc Ni is the most prominent among the three systems, while the smallest effect is found for fcc Co. This difference in demagnetization properties is consistent with the difference in observed values of $T_C$, which is highest for fcc Co and smallest for fcc Ni.

The next step in our comparison is to choose values of the pulse fluence for all three materials resulting in the same (or very similar) magnetization dynamics. As one can see in Fig. 5b, due to very different $T_C$ values, we need a pulse fluence as high as 45 J/m² for cobalt to exhibit a demagnetization that is similar to nickel for a fluence of 4.2 J/m². It can also be seen from Fig. 5b that in this case the magnetization dynamics of all three systems are very similar. The remaining differences are minute, for example, the minima of cobalt's magnetization curve appears somewhat later than the one of iron and nickel.

It is important to take into account that this very similar dynamics is observed for specific parameter values, such as the Gilbert and lattice damping. For example, to obtain the curve for cobalt presented in Fig. 5b, we used the damping value from Ref.³⁵, see Table I in Supplementary materials. However, with the damping value obtained from *ab-initio* calculations, the magnetization dynamics will be somewhat different. In particular, we then observe much slower remagnetization of cobalt (as, for example, shown in Fig. 4 obtained with the Gilbert damping set to 0.0014). Overall, our findings are in line with studies based on the Landau-Lishits-Bloch approach¹⁹, where the importance of the coupling between spin and electron systems was concluded.

In Fig. 6 we plot the maximally demagnetized state (minimum value of the magnetization curve, $M_{min}/M_0$) as function of absorbed laser pulse fluence for bcc Fe, fcc Co, and fcc Ni for the same value of the lattice damping (solid lines). The impact of fluence on the demagnetization was reported before for several materials.¹⁶⁻²¹ Here we systematically compare the impact of fluence on 3d ferromagnets and study the performance of the HC3TM and the impact of system parameters on the relation between demagnetization and laser fluence. The first observation is that one needs the lowest laser fluence to demagnetize fcc Ni and the strongest for fcc Co. This can, as mentioned earlier in this paper, be connected to the corresponding $T_C$. Secondly, $M_{min}/M_0$ in the simulations depends linearly on laser fluence, something which was observed experimentally in Ref.²⁰ and also for other materials in Refs.²¹,²².

In our work we use Gilbert damping that is the highest for nickel and the lowest for cobalt, which might add to the observed trend in Fig. 6. In particular, a higher Gilbert damping, which in the HC3TM can be correlated to the electron-magnon coupling, leads to a stronger heat transfer from the electronic to the spin subsystem, and therefore faster demagnetization, as in the case of nickel.

In the literature, the reported Gilbert damping parameters of iron, nickel, and cobalt deviate between different experiments. For example, reported $\alpha$ values for iron are in the range of 0.0019–0.0072²⁴. Nickel does have a consistently larger Gilbert damping than the other two elements, stemming from the $1/m_{Ni}$ scaling and the peak in the minority spin channel at the Fermi level²⁴. Varying the Gilbert damping in our HC3TM model results in curves that deviate slightly from those presented in Fig. 6 (data not shown) but the overall behavior still remains the same.

To demonstrate the impact of lattice damping on the trend of the demagnetization more in detail, we have performed spin-lattice simulations for all three elements using several different values of the lattice damping. In addition to the simulations using almost the same lattice damping value ($v \approx 1 \text{ ps}^{-1}$—solid lines in Fig. 6), we also used lattice damping values corresponding to the upper and lower limits of the electron-phonon coupling $G_{ep}$ found in the literature⁴,¹²,¹³.

These intervals for the reported electron-phonon couplings are shown as the shaded areas in Fig. 6 where it can be noticed that the trend of largest demagnetization for nickel and smallest demagnetization for cobalt holds

<figure><img src="image_7.png" /><figcaption>**Figure 5.** (a) Ultrafast demagnetization in fcc Ni, fcc Co and bcc Fe for the same value of the laser pulse fluence 8 J/m². (b) Ultrafast demagnetization in fcc Ni, fcc Co and bcc Fe for a different values of laser pulse fluence. The fluence is chosen to obtain the same demagnetization amplitude.</figcaption></figure>
<figure><img src="image_8.png" /><figcaption>**Figure 6.** The dependence of the maximally demagnetized value versus laser fluence for bcc Fe, fcc Co and fcc Ni. The shaded area demonstrate a range of reported⁴,¹²,¹³ lattice damping parameters obtained using the relation between lattice damping and $G_{ep}$ from¹¹.</figcaption></figure>

over almost the full interval of the shaded areas. For all three elements, the lower bounds of the shaded areas in Fig. 6 correspond to the lowest lattice damping while the upper bounds indicate the result of the largest lattice damping of the considered values. This further exemplifies the general picture that a stronger electron-lattice coupling transfers more heat from the electron to the lattice sub-system thus effectively reducing the amount of heat that gets pumped into the electron system.

In Ref.²⁰ Scheid et al. experimentally showed that the maximum demagnetization amplitude in Fe and Co are surprisingly similar. Various reasons behind the distribution of demagnetization rates in iron, cobalt, nickel were discussed, including the impact of electron-phonon coupling, the electron-magnon scattering rate, and the ultrafast light-induced quenching of the interatomic exchange. In our model, as seen in Fig. 6, this similarity could be explained by a very large lattice damping in Fe relatively to Co. Indeed, in HC3TM, and by using first principle values of Gilbert damping, to obtain similar slopes for iron and cobalt, we need to use a lattice damping for iron 10 times larger than for cobalt. This is illustrated in Fig. 6 where the dashed curves indicate a lattice damping of $v \approx 5.0 \text{ ps}^{-1}$ for iron and $v \approx 0.5 \text{ ps}^{-1}$ for cobalt. However recent experimental¹³ and first principles²⁵ results suggest that the lattice damping in cobalt is in fact larger than in iron. This would amplify the demagnetization in iron, while minimizing the one of cobalt. That conclusion of Scheid et al.²⁰ is therefore in line with our findings, because even for the same lattice damping for all systems (which in our model corresponds to electron-phonon coupling, see Fig. 1) we clearly obtain the highest demagnetization rate for nickel followed by iron and then cobalt (see solid lines in Fig. 6). One may further improve the agreement with the experimental observations by considering additionally ultrafast light-induced quenching of the interatomic exchange as suggested in Ref.²⁰.

## Discussion and conclusion

Using a heat-conserving three-temperature model we have calculated spin, lattice, and electron temperatures in simulations of the ultrafast magnetization dynamics of fcc cobalt and bcc iron. We have studied in detail the impact of lattice and Gilbert damping and compared results of the magnetization dynamics for bcc Fe and fcc Co with new data for fcc Ni, a system that was previously reported in Ref.⁴. It was found that in studies of ultrafast magnetic phenomena, the lattice dynamics plays a surprisingly important role in the HC3TM, even in cases when the direct spin-lattice coupling is not significant (i.e., when Eq. (8) is neglected). As analyzed here this is caused by how heat is transferred between the different reservoirs in the simulations, and how this heat transfer is influenced by dissipation mechanisms. To be concrete, we find here that the spin temperature used in the simulations depend on several mechanisms, one being how energy is dissipated from the lattice system. An experimental detection of this could possibly be to compare the magnetization dynamics of different isotopes of the systems analyzed here.

Simulations for various laser fluences were performed to study and compare ultrafast demagnetization of iron, cobalt, and nickel in various regimes. We have demonstrated that the simulations are consistent with observations, which show a linear trend between the maximally demagnetized state and laser fluence. We furthermore show that the experimentally found trend that for the same fluence, fcc Ni demagnetizes the most and fcc Co the least, holds for a wide range of realistic choices of the electron-phonon coupling strength.

## Methods. Atomistic spin-lattice dynamics simulations

As reported in Ref.²³, coupled spin and lattice dynamics can be obtained by solving the following Langevin equations of motion
$$ \frac{d\mathbf{m}_i}{dt} = - \frac{\gamma}{(1 + \alpha^2)} \mathbf{m}_i \times (\mathbf{B}_i + \mathbf{B}_i^\text{fl}) - \frac{\gamma}{(1 + \alpha^2)} \frac{\alpha}{m_i} \mathbf{m}_i \times (\mathbf{m}_i \times [\mathbf{B}_i + \mathbf{B}_i^\text{fl}]) \quad (2) $$

$$ \frac{d\mathbf{u}_k}{dt} = \mathbf{v}_k \quad (3) $$

$$ \frac{dv_k}{dt} = \frac{\mathbf{F}_k}{M_k} + \frac{\mathbf{F}_k^\text{fl}}{M_k} - \nu \mathbf{v}_k, \quad (4) $$

where $\mathbf{m}_i$ represents an atomic magnetic moment, $\mathbf{m}_i$ and $\gamma$ are the saturation magnetization and the gyromagnetic ratio correspondingly. Atomic displacements are denoted by $\mathbf{u}_k$, and velocities are $\mathbf{v}_k$. We obtain an effective exchange field $\mathbf{B}_i = -\partial H_{\text{SLD}}/\partial \mathbf{m}_i$ from the spin-lattice Hamiltonian, $H_{\text{SLD}}$. The Hamiltonian used in this work includes magnetic, lattice and spin-lattice coupling parts following Ref. 23:

$$ H_{\text{SLD}} = H_{\text{LL}} + H_{\text{SS}} + H_{\text{SSL}} \quad (5) $$

Here the lattice Hamiltonian reads

$$ H_{\text{LL}} = \frac{1}{2} \sum_{kl} \Phi_{kl}^{\mu\nu} u_k^{\mu} u_l^{\nu} + \frac{1}{2} \sum_k M_k v_k^{\mu} v_k^{\nu}, \quad (6) $$

where $\Phi_{kl}^{\mu\nu}$ is the force constant tensor, and $M_k$ is the mass of atom k. The magnetic Hamiltonian is described by

$$ H_{\text{SS}} = -\frac{1}{2} \sum_{ij} J_{ij}^{\alpha\beta}(0) m_i^{\alpha} m_j^{\beta}, \quad (7) $$

where $J_{ij}^{\alpha\beta}(0)$ is the exchange tensor at the equilibrium lattice positions. $\alpha, \beta$ denote Cartesian components in spin space, while $\mu, \nu$ corresponds to Cartesian components in real space. The magnetic anisotropy is not taken into account in our simulations since for the considered systems it is known to be on $\mu eV/\text{atom range}$, and verified to have only a negligible impact on the magnetization dynamics. We note, however, that for materials with larger spin-orbit coupling, the inclusion of an anisotropy term in Eq. (7) may be relevant.

It was shown in 23 that one of the ways to include spin-lattice coupling is to account for the dependence of exchange interaction values on atomic displacements $\mathbf{u}_k$, resulting in functions of the type $J_{ij}^{\alpha\beta}(\mathbf{u}_k)$ (see Ref. 23). As the displacements $\mathbf{u}_k$ are usually small (see Supplementary materials), one can write the Taylor expansion of the bi-linear magnetic Hamiltonian with respect to the lattice displacements up to the first order, resulting in spin-lattice coupling term bi-linear in spin and linear in displacements 23. The spin-lattice part of the Hamiltonian can then be written as

$$ H_{\text{SSL}} = -\frac{1}{2} \sum_{ijk} \Gamma_{ijk}^{\alpha\beta\mu} m_i^{\alpha} m_j^{\beta} u_k^{\mu}, \quad (8) $$

where we introduce the coupling constant $\Gamma_{ijk}^{\alpha\beta\mu} = \partial J_{ij}^{\alpha\beta}/\partial u_k^{\mu}$. We note that here only the isotropic part ($\alpha = \beta$) of the resulting tensor $\Gamma_{ijk}$ is considered for sites $k = i, j$, as the terms for $k \neq i, j$ are at least one order of magnitude smaller$^4$. Also note that to allow transfer of angular momentum off-diagonal terms should be included, diagonal terms of tensor $\Gamma_{ijk}^{\alpha\beta\mu}$ lead only to exchange striction$^{23}$. More complete and rotationally invariant Hamiltonian expressions for the atomistic spin-lattice coupling can be found in, e.g.$^{36,37}$.

All parameters for the spin-lattice dynamics simulations were obtained from ab-initio calculations (see Supplementary materials) for details). This includes exchange interactions, magnetic moments, inter atomic forces, and spin-lattice couplings. The force at site k is defined by $\mathbf{F}_k = -\partial H_{\text{SLD}}/\partial \mathbf{u}_k$. Gilbert and lattice damping constants are denoted $\alpha$ and $\nu$, respectively. In these types of Langevin simulations one employs stochastic fields, $\mathbf{B}_i^\text{fl}$ and $\mathbf{F}_k^\text{fl}$, as white noise with properties $(\mathbf{B}_{i,\mu}^\text{fl}(t)\mathbf{B}_{i,\nu}^\text{fl}(t')) = 2D_M\delta_{ij}\delta_{\mu\nu}\delta(t-t')$ and $(\mathbf{F}_{i,\mu}^\text{fl}(t)\mathbf{B}_{i,\nu}^\text{fl}(t')) = 2D_L\delta_{kl}\delta_{\mu\nu}\delta(t-t')$. We utilise $D_M = \alpha k_B T_e/\gamma m$, $D_L = \nu M k_B T_e$, where $k_B$ is Boltzmann constant, and $T_e$ is electronic temperature respectively (for details see e.g., Ref. 38). The formalism above is implemented in the UppASD$^{39}$ code which was used for all simulations in this work.

## Data availability

The data supporting the findings of this study are available in the paper and upon reasonable request from the corresponding authors (M.Pa. and I.M.).

Received: 19 October 2023; Accepted: 2 April 2024

Published online: 07 April 2024

## References

1. Beaurepaire, E., Merle, J.-C., Daunois, A. & Bigot, J.-Y. Ultrafast spin dynamics in ferromagnetic nickel. Phys. Rev. Lett. 76, 4250 (1996).

2. Kirilyuk, A., Kimel, A. V. & Rasing, T. Ultrafast optical manipulation of magnetic order. Rev. Mod. Phys. 82, 2731-2784. https://doi.org/10.1103/RevModPhys.82.2731 (2010).

3. Scheid, P., Remy, Q., Lebègue, S., Malinowski, G. & Mangin, S. Light induced ultrafast magnetization dynamics in metallic compounds. J. Magn. Magn. Mater. 560, 169596. https://doi.org/10.1016/j.jmmm.2022.169596 (2022).
4. Pankratova, M. et al. Heat-conserving three-temperature model for ultrafast demagnetization in nickel. *Phys. Rev. B* **106**, 174407. https://doi.org/10.1103/PhysRevB.106.174407 (2022).

5. Locht, I. L. M., Di Marco, I., Garnerone, S., Delin, A. & Battiato, M. Ultrafast magnetization dynamics: Microscopic electronic configurations and ultrafast spectroscopy. *Phys. Rev. B* **92**, 064403. https://doi.org/10.1103/PhysRevB.92.064403 (2015).

6. Shim, J.-H. et al. Role of non-thermal electrons in ultrafast spin dynamics of ferromagnetic multilayer. *Sci. Rep.* **10**, 6355 (2020).

7. Battiato, M., Carva, K. & Oppeneer, P. M. Superdiffusive spin transport as a mechanism of ultrafast demagnetization. *Phys. Rev. Lett.* **105**, 027203. https://doi.org/10.1103/PhysRevLett.105.027203 (2010).

8. Willems, F. et al. Optical inter-site spin transfer probed by energy and spin-resolved transient absorption spectroscopy. *Nat. Commun.* **11**, 1-7. https://doi.org/10.1038/s41467-020-14691-5 (2020).

9. Carpine, E. et al. Dynamics of electron-magnon interaction and ultrafast demagnetization in thin iron films. *Phys. Rev. B* **78**, 174422. https://doi.org/10.1103/PhysRevB.78.174422 (2008).

10. Ryan, S. A. et al. Optically controlling the competition between spin flips and intersite spin transfer in a Heusler half-metal on sub-100 fs timescales. *arXiv: 2305.16455* (2023).

11. Ma, P.-W., Dudarev, S. L. & Woo, C. H. Spin-lattice-electron dynamics simulations of magnetic materials. *Phys. Rev. B* **85**, 184301. https://doi.org/10.1103/PhysRevB.85.184301 (2012).

12. Zahn, D. et al. Lattice dynamics and ultrafast energy flow between electrons, spins, and phonons in a 3d ferromagnet. *Phys. Rev. Res.* **3**, 023032. https://doi.org/10.1103/PhysRevResearch.3.023032 (2021).

13. Zahn, D. et al. Intrinsic energy flow in laser-excited 3d ferromagnets. *Phys. Rev. Res.* **4**, 013104. https://doi.org/10.1103/PhysRevResearch.4.013104 (2022).

14. Koopmans, B. et al. Explaining the paradoxical diversity of ultrafast laser-induced demagnetization. *Nat. Mater.* **9**, 259-265 (2010).

15. You, W. et al. Revealing the nature of the ultrafast magnetic phase transition in ni by correlating extreme ultraviolet magneto-optic and photoemission spectroscopies. *Phys. Rev. Lett.* **121**, 077204. https://doi.org/10.1103/PhysRevLett.121.077204 (2018).

16. Hohlfeld, J., Matthias, E., Knorren, R. & Bennemann, K. Nonequilibrium magnetization dynamics of nickel. *Phys. Rev. Lett.* **78**, 4861 (1997).

17. Mueller, B. et al. Feedback effect during ultrafast demagnetization dynamics in ferromagnets. *Phys. Rev. Lett.* **111**, 167204 (2013).

18. Roth, T. et al. Temperature dependence of laser-induced demagnetization in Ni: A key for identifying the underlying mechanism. *Phys. Rev. X* **2**, 021006 (2012).

19. Atxitia, U., Chubykalo-Fesenko, O., Walowski, J., Mann, A. & Münzenberg, M. Evidence for thermal mechanisms in laser-induced femtosecond spin dynamics. *Phys. Rev. B* **81**, 174401. https://doi.org/10.1103/PhysRevB.81.174401 (2010).

20. Scheid, P. et al. Uncovering the role of the light-induced reduction of the inter-atomic exchange in the ultrafast demagnetization of Fe, Co, and Ni. *arXiv preprint arXiv:2301.01055* (2023).

21. Windsor, Y. W. et al. Exchange scaling of ultrafast angular momentum transfer in 4f antiferromagnets. *Nat. Mater.* **21**, 514-517 (2022).

22. Mishra, K. et al. Ultrafast demagnetization in a ferrimagnet under electromagnetic field funneling. *Nanoscale* **13**, 19367-19375. https://doi.org/10.1039/D1NR04308K (2021).

23. Hellsvik, J. et al. General method for atomistic spin-lattice dynamics with first-principles accuracy. *Phys. Rev. B* **99**, 104302. https://doi.org/10.1103/PhysRevB.99.104302 (2019).

24. Lu, Z. et al. Influence of nonlocal damping on magnon properties of ferromagnets. *Phys. Rev. B* **108**, 014433. https://doi.org/10.1103/PhysRevB.108.014433 (2023).

25. Ritzmann, U., Oppeneer, P. M. & Maldonado, P. Theory of out-of-equilibrium electron and phonon dynamics in metals after femtosecond laser excitation. *Phys. Rev. B* **102**, 214305. https://doi.org/10.1103/PhysRevB.102.214305 (2020).

26. Duffy, D. & Rutherford, A. Including the effects of electronic stopping and electron-ion interactions in radiation damage simulations. *J. Phys.: Condens. Matter* **19**, 016207 (2006).

27. Duffy, D. & Rutherford, A. Including electronic effects in damage cascade simulations. *J. Nuclear Mater.* **386**, 19-21 (2009).

28. Ma, P.-W., Dudarev, S. L., Semenov, A. A. & Woo, C. H. Temperature for a dynamic spin ensemble. *Phys. Rev. E* **82**, 031111. https://doi.org/10.1103/PhysRevE.82.031111 (2010).

29. Casas-Vázquez, J. & Jou, D. Nonequilibrium temperature versus local-equilibrium temperature. *Phys. Rev. E* **49**, 1040-1048. https://doi.org/10.1103/PhysRevE.49.1040 (1994).

30. Unikandanunni, V. et al. Anisotropic ultrafast spin dynamics in epitaxial cobalt. *Appl. Phys. Lett.* **118**, 232404. https://doi.org/10.1063/5.0049692 (2021).

31. Sadhukhan, B., Bergman, A., Kvashnin, Y. O., Hellsvik, J. & Delin, A. Spin-lattice couplings in two-dimensional $cri_3$ from first-principles computations. *Phys. Rev. B* **105**, 104418. https://doi.org/10.1103/PhysRevB.105.104418 (2022).

32. Waldecker, L., Bertoni, R., Ernstorfer, R. & Vorberger, J. Electron-phonon coupling and energy flow in a simple metal beyond the two-temperature approximation. *Phys. Rev. X* **6**, 021003. https://doi.org/10.1103/PhysRevX.6.021003 (2016).

33. Maldonado, P. et al. Tracking the ultrafast nonequilibrium energy flow between electronic and lattice degrees of freedom in crystalline nickel. *Phys. Rev. B* **101**, 100302. https://doi.org/10.1103/PhysRevB.101.100302 (2020).

34. Hiramatsu, R., Miura, D. & Sakuma, A. Temperature variation of gilbert damping constant in 3d transition metals from first principles. *J. Phys. Soc. Jpn.* **92**, 044704. https://doi.org/10.7566/jpsj.92.044704 (2023).

35. Schoen, M. A. W. et al. Magnetic properties in ultrathin 3d transition-metal binary alloys. II. experimental verification of quantitative theories of damping and spin pumping. *Phys. Rev. B* **95**, 134411. https://doi.org/10.1103/PhysRevB.95.134411 (2017).

36. Weißenhofer, M. et al. Rotationally invariant formulation of spin-lattice coupling in multiscale modeling. *Phys. Rev. B* **108**, L060404. https://doi.org/10.1103/PhysRevB.108.L060404 (2023).

37. Sadhukhan, B., Bergman, A., Hellsvik, J., Thunström, P. & Delin, A. Spin-lattice couplings in a skyrmion multilayers of Pd-Fe/Ir(111). *arXiv:2309.03074* (2023).

38. Eriksson, O., Bergman, A., Bergqvist, L. & Hellsvik, J. *Atomistic Spin Dynamics: Foundations and Applications* (Oxford University Press, 2017).

39. Uppsala atomistic spin dynamics. https://github.com/UppASD/UppASD.

## Acknowledgements

This work was financially supported by the Knut and Alice Wallenberg Foundation (Grant Numbers 2018.0060, 2021.0246, and 2022.0108), Vetenskapsrådet (Grant Numbers 2019-03666, 2016-05980, and 2019-05304), the European Research Council (Grant Number 854843-FASTCORR), the foundation for Strategic Research SSF, and Ölle Engkvist foundation. Support from STandUP and eSSENCE is also acknowledged. O.E. and A.D. also acknowledge support from the Wallenberg Initiative Materials Science for Sustainability (WISE) funded by the Knut and Alice Wallenberg Foundation (KAW). Computations were enabled by resources provided by the National Academic Infrastructure for Supercomputing in Sweden (NAISS, projects 2023/1-10, 2023/5-454, and 2023/1-44) at both the National Supercomputing Centre (NSC, Tetralith) and the PDC Centre for High Performance Computing (PDC-KTH, Dardel), partially funded by the Swedish Research Council through Grant Agreement No. 2018-05973. P.S. acknowledges support from the ANR-20-CE09-0013 UFO.
Author contributions

M.Pa., I.M., and A.B. initially designed the project. M.Pa. and I.M. performed the simulations and calculations, M.Pa., I.M., P.S., O.E., and A.B. analyzed the results and their interpretations. M.Pa., and I.M. wrote the initial version of the manuscript, M.Pa., I.M., D.T. M.Pe., A.D., E.S. P.S., O.E., A.B. critically discussed the results and contributed to preparing the paper.

Funding

Open access funding provided by Uppsala University.

Competing interests

The authors declare no competing interests.

Additional information

**Supplementary Information** The online version contains supplementary material available at https://doi.org/10.1038/s41598-024-58662-y.

**Correspondence** and requests for materials should be addressed to M.P. or I.P.M.

Reprints and permissions information is available at www.nature.com/reprints.

**Publisher's note** Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

<figure><img src="image_9.png" /></figure>

**Open Access** This article is licensed under a Creative Commons Attribution 4.0 International

License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/.

© The Author(s) 2024