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

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.


I. INTRODUCTION
Ultrafast demagnetization was discovered by Beaurepaire and coauthors in 1996 [1].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 [2].In the same pioneering work of Beuarepaire [1], these experimental observations were interpreted using a three-temperature model (3TM) [1], which assumes three thermalised reservoirs, in particular, spin-, lattice, and electron reservoirs, that can exchange energy through coupling parameters (via electron-phonon, electron-spin, and spin-lattice coupling).The 3TM is often used to interpret ultrafast magnetization dynamics processes [1,3,4].Recently, many other models have been proposed [2,3,5,6] to describe possible mechanisms of ultrafast demagnetization, such as the importance of spin-dependent transport of laserexcited electrons [7], the optical inter-site spin transfer effect (OISTR) that was considered in Ref. [8], and the Elliott-Yafet electron-phonon spin-flip scattering, studied in Ref. [9].In [10] 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 the [10] is also establishing the relation between the dissipative parameters entering the Langevin equations for the lattice and spin degrees of freedom and the heat transfer coefficients of 3TM.
In addition, Zahn et al. [11,12] proposed an energyconserving model based on atomistic spin dynamics simulations.Using this model they studied ultrafast demagnetization of nickel [11], iron and cobalt [12].However, the model of Zahn et al. relies on the microscopic electron-lattice heat transfer coefficient which is hard to estimate.Literature values for the electron-phonon coupling parameter G ep can vary for nickel up to one order of magnitude [4,11].For iron, reported values of the electron-phonon coupling parameter differs from 7 • 10 17 to 5.48 • 10 18 W/m 3 K [12], and for cobalt from 6 • 10 17 to 4.05•10 18 W/m 3 K (please see [12] and references therein).
One of the reasons for this uncertainty of the electronlattice heat transfer coefficient is that the electronphonon coupling is often deduced from observables, without taking into account the energy cost of demagnetization.Also, these parameters are often estimated indirectly using, for example, a two-temperature model [13].It was demonstrated that G ep extracted from the data employing the two-temperature model is twice smaller than that calculated from ab initio theory [13].These differences in reported G ep values makes the interpretation of the experimental observations very challenging.
To address this, a previous work proposed a heat-arXiv:2308.08996v1[cond-mat.mtrl-sci]17 Aug 2023 conserving three-temperature model (HC3TM) [4] for calculation on spin, lattice and electron temperatures during atomistic spin-dynamics simulations.The benefit of this model is that it reduces the dependence on heattransfer parameters that are hard to estimate, such as electron-phonon, G ep , electron-spin, G es , and the spinlattice coupling, G sl .In addition, this model allows for a description of the nickel demagnetization rate on subpicosecond timescales.
In the present work, we apply the HC3TM model to study ultrafast demagnetization of bcc Fe and fcc Co and compare with previously published data on ultrafast dynamics of fcc Ni [4].The ultrafast demagnetization of these three ferromagnets was experimentally studied in Ref. [14].It was observed that the amplitude of the demagnetization value increases linearly with the fluence of the light for Fe, Ni, and Co.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 inter-atomic exchange.In the present work, we aim to study the impact of fluence on the demagnetization of 3d ferromagnets using the HC3TM.Moreover, we investigate in more detail the impact of various factors on ultrafast demagnetization, such as Gilbert damping, lattice damping, microscopic spin-lattice coupling, and laser fluence.

II. ATOMISTIC SPIN-LATTICE DYNAMICS SIMULATIONS
Spin and lattice dynamics is obtained from coupled atomistic spin-and lattice dynamics (SLD) simulations, using Langevin dynamics simulations, as proposed in Ref. [15]: where m i represents an atomic magnetic moment, m i and γ are the saturation magnetization and the gyromagnetic ratio correspondingly.Atomic displacements are denoted by u k , and velocities are v k .We obtain an effective exchange field B i = −∂H SLD /∂m i from the spin-lattice Hamiltonian, H SLD .The Hamiltonian used in this work includes magnetic, lattice and spin-lattice coupling parts following Ref.[15]: Concretely: where the force constant tensor Φ µν kl is a rank 2 tensor in real space, and M k is the mass of atom k.The magnetic Hamiltonian is described by where J αβ ij (0) is the exchange tensor at the equilibrium lattice positions.α, β ∈ {x, y, z} denote Cartesian components in spin space, while µ, ν ∈ {x, y, z} corresponds to Cartesian components in real space.
It was recently suggested that spin-lattice coupling can be represented in coupled spin-lattice simulations, by a dependence of the exchange interaction on atomic displacements u k , resulting in functions of the type J αβ ij (u k ) (see Ref. [15]).As the values of u k are usually small, the spin-lattice term in the Hamiltonian can be obtained by Taylor expanding the magnetic bi-linear Hamiltonian with respect to the lattice displacements.This results in a spin-lattice coupling term that is bi-linear in spin and linear in displacements [15].Concretely, where we introduce the coupling constant Γ αβµ ijk ≡ Γ αβµ ijk (u µ k ) = ∂J αβ ij (u µ k )/∂u µ k .All parameters for the spin-lattice dynamics simulations were obtained from ab-initio calculations (see Appendix B for details).This includes exchange interactions, magnetic moments, inter atomic forces, and spin-lattice couplings.
The force at site k is defined by F k = −∂H SLD /∂u k .Gilbert and lattice damping constants are denoted α and ν, respectively.In these types of Langevin simulations one employs stochastic fields, B f l i and F f l k , as white noise with properties In our simulations, we use D M = αk B T /γm, D L = νM k B T , where T and k B are temperature and Boltzmann constant respectively (for details see e.g.Ref. [16]).

III. HEAT-CONSERVING THREE-TEMPERATURE MODEL
The heat-conserving three-temperature model was proposed in Ref. [4] for the study of ultrafast demagnetisation of fcc Ni.Similar to Beurepaire's three-temperature model, it is assumed that lattice, and spin systems are connected to a finite electronic heat bath.Both spin and lattice subsystems have their own damping; Gilbert or lattice damping, correspondingly.The model does not explicitly rely on the heat transfer parameters of the conventional 3TM [1].The difference between HC3TM and the conventional 3TM is easy to see from the schematic picture in Fig. 1.Overall, instead of microscopic heat transfer coefficients, as used in 3TM, the HC3TM is formulated in terms of atomistic damping coefficients.Thus, the HC3TM calculates the heat transfer between the three reservoirs (spin, lattice and electron) [1] directly from the simulations.In fact, all parameters of this theory can be calculated using, for example, density functional theory, although we note that in the present work we varied the lattice damping parameter in a range of realistic values, as opposed to performing an explicit calculation of it.The relation between the lattice damping and G ep was established in [10] using spin-lattice-electron three-temperature model.It was shown that the lattice damping is directly proportional to G ep , which was previously also reported by Duffy et el. in [17,18].Another distinction of HC3TM is that (unlike 3TM) the temperature of the different subsystems is calculated at every time step of the spin-lattice simulation using the expression: where T l is calculated from the average kinetic energy of the lattice vibrations; ⟨E kin l ⟩/k B .The spin temperature, T s is calculated according to [19], i.e. using , where mi is the normalized local spin moment.The definition of an instantaneous measure of the temperatures in an out-of-equilibrium system is not obvious [20], but is still typically assumed to hold in these kinds of models.Introducing a time average of the temperatures, for narrow time windows, does however not change the results significantly.In addition, Eq. ( 8) contains the temperature dependent specific heats of the electron-, lattice-and spin systems, C e (T ), C l (T ) and C s (T ), respectively.The calculations details of the specific heat capacities are given in Appendix A. The possibility for an external stimulus, e.g. a laser, to momentarily increase the temperature of the electronic subsystem is captured by the source term W (t) of Eq. ( 8), which is modelled as a Gaussian function: , where t 0 is the center, σ = 0.02 ps is the width of the pulse, and W 0 is the magnitude (in Joules), correspondingly.More details concerning the HC3TM and its comparison with 3TM can be found in Ref. [4].

A. HC3TM for cobalt and iron
In this section we present results of spin-lattice dynamics simulations of fcc Co and bcc Fe.In both cases, we run simulations for cells with a 60 × 60 × 60 repetition of the fcc (bcc) unit cell, using periodic boundary conditions.In addition, we used N t = 1 × 10 6 time steps of dt = 10 −16 s, in combination with other parameters presented in Table I.
We start by describing results from ultrafast demagnetization of fcc Co.First, we present temperature profiles (Fig. 2) for spin, lattice, and electron subsystems, for a laser fluence of 45 J/m 2 , together with the corresponding magnetization dynamics.Other parameters relevant for these simulations are listed in Table I.It can be seen from the figure that similar to experimental observations [21] and theoretical studies [12] 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. [12], using the parametrized three-temperature model, features of the magnetization dynamics similar to the data of Fig. 2 were also observed, however for an adjusted Gilbert damping value of 0.1, which is unrealistically large.In contrast, in the calculations presented here we use a value of damping from ab-inito calculations (that is close to experiment), that as seen in Table I is much smaller than the value used in Ref. [12].Similar spin-lattice simulations using the HC3TM, were performed for bcc Fe.Just like for fcc Co we first present magnetization dynamics after the absorption of a laser pulse and the corresponding spin, lattice, and electronic temperatures, see Fig. 3. 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 [9], for the same pulse fluence (M/M 0 is equal around 0.9 in our calculations).

B. Results from coupled spin-lattice simulations
To study the impact of direct spin-lattice coupling on the magnetization dynamics, given by Eq. ( 7), we consider explicitly the spin-lattice coupling, where we obtain the elements of the tensor, Γ ijk in the spin-lattice coupling, from calculations based on density functional theory (e.g. as described in Ref. [15]).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 in Eq. ( 7) is included or neglected, and the results that follow were obtained from the HC3TM with Eq. ( 7) included explicitly in the simulations.
In Fig. 4 we show for bcc Fe the demagnetization profile from coupled spin-lattice dynamics simulations in combination with the HC3TM.The simulations are, as noted, done with and without the coupling term in Eq. (7).It can be seen from Fig. 4 that spin-lattice coupling impacts the magnetization dynamics very 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.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.(8).A possible explanation for the low influence of the spin-lattice couplings can be drawn as follows: even at the maximum spin temperature in the simulations (T s ≲ 600 K for Fe), the spin system is far away from T c , i.e. with a general strong coupling between spin moments.Thus, the small average displacements (in the order of 10 −2 a) and their distribution throughout the sample result in the summed effect of Eq. ( 7) being only a minor perturbation on the total exchange field.
Also, for bcc Fe the resulting profile is in rather good agreement with experimental observations.For materials with larger spin-lattice coupling [24] it can be expected that the importance of Eq. ( 7) increases, and our simulations show that for bcc Fe the spin-lattice coupling becomes more important for higher pulse fluences (data not shown).
While the impact of spin-lattice coupling, via Eq.( 7), on the magnetization dynamics is only marginal for the elemental ferromagnts, the impact of the dynamical properties of the lattice itself, and in particular, the lattice damping, is significant.We illustrate this for fcc Co, in Fig. 5.Here it is seen that a decrease of the lattice damping (ν in Eq. ( 3) leads to an increase of the demagnetization; both the magnetization minimum as well as the reduced magnetization ( M M0 ) are influenced by the lattice damping parameter.This can best be 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. 5, 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. 5) 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.One of the important conclusions 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 com-paring 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.

C. 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 [4].The results are summarized in Figs. 6 and 7. We begin with comparing ultrafast magnetization dynamics of fcc Co, fcc Ni, and bcc Fe, for the same laser pulse fluence, presented in Fig. 6.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. 7, due to very different T c values, we need a pulse fluence as high as 45 J/m 2 for cobalt to exhibit a demagnetization that is similar to nickel for a fluence of 4.2 J/m 2 .It can also be seen from Fig. 7 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. 7, we used the damping value from Ref. [22], see Table I.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. 5 obtained with the Gilbert damping set to 0.0014).
In Fig. 8 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 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 ex-FIG.8.
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 [4,11,12] lattice damping parameters obtained using the relation between lattice damping and Gep from [10].
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. 8.In particular, a higher Gilbert damping, which in the HC3TM corresponds 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 α values for iron are in the range of 0.0019 − 0.0072 [23].Nickel does have a consistently larger Gilbert damping than the other two elements, stemming from the 1/m N i scaling and the peak in the minority spin channel at the Fermi level [23].Varying the Gilbert damping in our HC3TM model results in curves that deviate slightly from those presented in Fig. 8 (data not shown) but the overall behaviour 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 for the same lattice damping value (ν = 1 × 10 −13 kg/s -solid lines in Fig. 8), we also used lattice damping values corresponding to the upper and lower limits of the electron-phonon coupling G ep found in the literature [4,11,12].These intervals for the reported electron-phonon couplings are shown as the shaded areas in Fig. 8 where it can be noticed that the trend of largest demagnetization for nickel and smallest demagnetization for cobalt holds over almost the full interval of the shaded areas.For all three elements, the lower bounds of the shaded areas in Fig. 8 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. [14] 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 inter-atomic exchange.In our model, as seen Fig. 8. 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, 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. 8 where the dashed curves indicate a lattice damping of ν = 5 × 10 −13 kg/s for iron and ν = 5 × 10 −14 kg/s for cobalt.However recent experimental [12] and first principles [27] 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. [14] 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. 8).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. [14].

V. CONCLUSIONS
Using a heat-conserving three-temperature model we have calculated spin, lattice, and electron temperatures in simulations of the ultrafast magnetization dynamics of cobalt and iron.We have studied in detail the impact of Gilbert and lattice damping and compared results of the magnetization dynamics for bcc Fe, fcc Co and fcc Ni.It was found that in studies of ultrafast magnetic phenomena, the lattice dynamics plays a surprisingly important role, even in cases when the direct spin-lattice coupling is not significant (i.e., when Eq. 7 is neglected).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 recent 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.
In our calculations we use temperature-dependent heat capacities for the lattice, electron, and spin subsystems.The corresponding curves are presented in Fig. 9.The lattice heat capacity is obtained from the Debye model.The spin heat capacity is calculated using a recently proposed approach from quantum statistics [28].The method more accurately describes heat capacities at low temperatures than what is obtained from Boltzmann statistics.However, its disadvantage is an overestimation of the heat capacity around T c , as can be seen from Fig. 9. Since in our calculations the spin temperature does not reach T c , this shortcoming of the method does not impact our results.Moreover, our calculations show that using this approach does not change any of our main conclusions or results.We can confirm that all results and conclusions such as linear dependence of demagnetization amplitude on fluence, etc remain the same regardless of heat capacities used in calculations (please see Appendix A in [4] for details of heat capacities calculations), including constant (i.e., temperature independent) heat capacities.The only thing that is affected by the change of capacities calculation is demagnetization amplitude.For the electronic subsystem, we assume that the electronic heat capacity is proportional to the electronic temperature C e = γ e T e , and for bcc Fe we use γ e = 4.9 Appendix B: Details of the ab-initio calculations Similarly to Ref. [4], the calculation of the interactions to parameterize Eqs. 1 and 3 were performed in the framework of Density Functional Theory (DFT), using two different methods: (i ) the real-space linear muffin-tin orbital method in the atomic sphere approximation (RS-LMTO-ASA) [29,30] and (ii ) the planewave pseudopotential-based Quantum ESPRESSO (QE) package [31] in combination with the PHONOPY [32] software.On one hand, this choice is motivated, essentially, on the well-known shortcoming of the LMTO-ASA method on obtaining accurate Hellmann-Feynman forces [33] -and, by consequence, phonon frequencies.On the other hand, to calculate the spin-lattice couplings (SLC), a supercell approach is necessary because of the breaking of inversion symmetry and long-range interactions [4].In this sense, RS-LMTO-ASA is more suitable for this calculation because it deals directly with real-space.
For the computation of the force constants using a combination of QE and PHONOPY, the fcc Ni case used the same definitions as described in Ref. [4].Analogously, in the case of bcc Fe and fcc Co, in QE the scalar relativistic schemes based on the projector augmentedwave method (PAW) [34] and the ultrasoft pseudopotentials (USPP) [35], respectively, are employed.As the exchange-correlation (XC) functional, we use the generalized gradient approximation (GGA) with the Perdew-Burke-Ernzerhof parametrization [36].Here, the choice of a spin-polarized GGA XC term is motivated by the fact that the calculated lattice constants as well as the phonon frequencies better reproduce the experimental data in comparison with the local spin density approximation [37].The iron and cobalt atoms were described by standard pseudopotentials from the QE library, with 3s, 3p, 4s, 4p and 3d valence electrons.A cutoff of 100 (1080) and 90 (500) Ry for the kinetic energy (charge density) was considered for fcc Co and bcc Fe, in this order.In turn, a Monkhorst-Pack (MP) [38] grid of 24 × 24 × 24 ⃗ k-points was set for the first Brillouin zone integration.The self-consistent calculations with a 10 −10 Ry convergence threshold were carried out using the Marzari-Vanderbilt [39] cold smearing with a spreading of 0.01 Ry.The force constants and phonon frequencies, computed by PHONOPY, were based on 6 × 6 × 6 supercells (216 atoms), where we considered a 4 × 4 × 4 ⃗ k-points mesh in the self-consistent cycle.
The exchange parameters, Gilbert damping (α), density of states and spin-lattice couplings (SLC) were obtained using the RS-LMTO-ASA method.In this realspace formalism, the eigenvalue problem is solved with the help of the recursion method [40], where the recursion chain is ended after LL steps by making use of the Beer-Pettifor terminator [41].For bcc Fe we considered LL = 31, while for fcc Co a much higher LL value (LL = 51) is needed to better describe the density of states and Green's functions at the Fermi level.In line with recent experience on finding a model to calculate reliable α's from first principles [23], the bulk systems consist of a big cluster in real-space containing ∼ 55000 (bcc) and ∼ 696000 (fcc) atoms, located in the perfect crystal positions.Differently from the QE calculations (in which the structures were optimized), the lattice parameters were fixed to the experimental values of a = 2.87 Å (bcc Fe) and a = 3.54 Å (fcc Co) [42].In this setup, the local density approximation (LSDA), with the von Barth and Hedin parametrization [43], was used.In particular, this choice of XC was based on the fact that some of quantities that parameterize Eqs.1-7 were already investigated using LSDA [23,[44][45][46].It is a known trait, however, that GGA produces quantities with similar quality for the 3d magnetic elements [23,46,47].Finally, the spinorbit coupling (SOC) is included as a l •s term, computed in each variational step [48,49].

FIG. 1 .
FIG. 1.The main idea behind the HC3TM (a) and 3TM (b) models.The text along the arrows indicates the parameters responsible for heat transfer in each model.For the parameters of the 3TM, see Ref. [1].

FIG. 2 .
FIG. 2. Spin, lattice and electron temperatures (a) and magnetization dynamics (b) of fcc Co, obtained with the HC3TM model.The pulse fluence is 45 J/m 2 .

FIG. 4 .
FIG. 4.Impact of spin-lattice coupling on magnetization dynamics of bcc Fe.The pulse fluence is 20 J/m 2 .

FIG. 6 .
FIG. 6. Ultrafast demagnetization in fcc Ni, fcc Co and bcc Fe for the same value of the laser pulse fluence 8 J/m 2 .

FIG. 9 .
FIG.9.Heat capacities for the electronic, spin, and lattice subsystems used in the calculations for bcc Fe.

TABLE I .
List of the parameters used in the simulations.