Ultrafast Demagnetization through Femtosecond Generation of Non-thermal Magnons

,


I. INTRODUCTION
The discovery that magnetic order can be manipulated on sub-picosecond timescales by femtosecond laser pulses [1][2][3] has fueled the emergence of intensive experimental and theoretical research efforts in the field of ultrafast magnetization dynamics.What makes this field particularly interesting, apart from its technological potential in future memory and spintronic devices [4,5], is that many well-established physical paradigms cannot be simply transferred from the equilibrium to the ultrafast regime, due to its highly non-equilibrium nature.Relatedly, albeit more than 25 years of intense research, the underlying mechanisms of ultrafast demagnetization are still heavily debated [6][7][8]: while some works [9][10][11][12][13][14] lean towards longitudinal excitations -i.e., the reduction of the magnetic moment carried by each atom due to the decrease of exchange splitting -others [15][16][17][18][19] hint at transverse spin excitations -a reduction of the average magnetization due to the mutual tilting of the moments carried by different atoms -as the main contribution.Non-local contributions due to superdiffusive spin currents [20,21] are relevant in certain situations [22][23][24][25].However, it has become evident that they are most likely not the only mechanism of ultrafast demagnetization [26,27].
Theoretical models describing ultrafast magnetization dynamics typically rely on a separation of electronic, phononic and -if magnetization dynamics are to be considered separately -spin degrees of freedom.Beaurepaire et al. [1] introduced the three-temperature model (3TM) to explain the flow of the energy transferred by the laser * markus.weissenhofer@fu-berlin.de by assuming that each subsystem is internally in thermal equilibrium and the system can hence be described by three temperatures (for electrons, phonons and spins), together with the respective distributions (Fermi-Dirac and Bose-Einstein).However, it was pointed out in numerous investigations that the distributions are non-thermal on ultrafast timescales [28][29][30][31][32][33][34][35][36][37].Also, the 3TM discards completely the transfer of angular momentum due to demagnetization, which, according to recent experiments [38,39], appears to be primarily to the lattice.
Transverse demagnetization is often studied using atomistic spin dynamics simulations based on the stochastic Landau-Lifshitz-Gilbert (LLG) equation together with an extended Heisenberg model [40][41][42], which can successfully reproduce experimentally measured demagnetization curves [43,44].The stochastic LLG is a Langevin-type equation with a coupling to a heat bath with given temperature via a single parameter, the Gilbert damping parameter.This parameter includes all possible contributions -Fermi surface breathing, crystal defects, coupling to phonons, s−d coupling, etc. [45][46][47][48][49][50][51][52] to damping and while it can in principle be obtained from ab initio calculations, in practice it is typically taken from experimental measurements of ferromagnetic resonance (FMR) [53].On the one hand, this ensures the versatility of atomistic spin dynamics simulations, but on the other hand, it obscures the details of the underlying microscopic energy and angular momentum transfer processes -which are crucial for understanding the fundamentals of ultrafast demagnetization.For this reason, steps have been taken in recent years to explicitly consider the coupling of spins to phonons [54][55][56][57][58][59][60][61][62] and electrons [63][64][65].Also, due to the classical nature of the commonly used stochastic LLG, the equilibrium magnon occupations cal-culated by it follow Rayleigh-Jeans rather than Bose-Einstein statistics, henceforth leading to the wrong temperature scaling of the magnetization [66,67].Implementation of quantum statistics in the spin-dynamics simulations can however provide the correct low-temperature scaling of the magnetization [68,69].
In this work, we investigate the laser-induced generation of magnons, the low energy transverse excitations of the spin system, due to electron-magnon scattering.We develop a quantum kinetic approach, which will be termed (N+2)-temperature model [(N+2)TM], to perform quantitative simulations of the time evolution of the non-thermal magnon dynamics in bcc iron.Being based on ab initio parameters and considering also nonthermal magnon distributions, our work goes well beyond what has been done in Refs.[63,64,70] and the conventional 3TM.In addition, we show that the 3TM and its relevant parameters can be obtained from our (N+2)TM and, with that, from ab initio calculations.Importantly, using ab initio calculated input parameters, our quantum kinetic theory predicts a sizable and ultrafast demagnetization of iron within 200 fs, in excellent agreement with experiments [15].

II. OUT-OF-EQUILIBRIUM MAGNON DYNAMICS MODEL
To describe the time evolution of the ultrafast nonthermal magnon occupation dynamics, we assume that their creation and annihilation is dominated by electronmagnon scattering processes.In this work, we use the sp − d model [71,72] to describe such processes.The basic idea of both s − d model and sp − d model is the separation of electrons in localized (d band) electrons and itinerant (s band or s and p bands) electrons.The magnetic moments of the d electrons make up the Heisenbergtype [73] magnetic moments of constant length, the small energy excitations of which are the magnons.The itinerant electrons are described within a Stoner-type model [74].While an unambiguous identification of sp and d electrons as localized and itinerant is strictly speaking not possible, it has nonetheless been established in literature that these models provide a suitable framework for the description of electron-spin interaction in many phenomena relevant for spintronics, e.g.magnetic relaxation [75][76][77], ultrafast demagnetization [63-65, 70, 78-80] and spin torques [81].
We assume local exchange between the itinerant and localized spins, as given by the Hamiltonian Ĥem ∼ N i=1 ŝitin • Ŝloc i , with N being the number of atoms, and ŝitin and Ŝloc i the spin operators for itinerant (sp) electrons and localized (d) electrons at atom i.In second quantization and second order in magnon variables (details in Method Section V A), the Hamiltonian reads Here, ∆ is the sp − d exchange parameter, S is the absolute value of the localized spins, k and q are vectors in reciprocal space, ĉ( †) kνσ is the fermionic electron annihilation (creation) operator for the itinerant electrons -with ν being the band index and σ ∈ {↑, ↓} -and b( †) q is the bosonic magnon annihilation (creation) operator.The first term in Equation ( 1) describes the spin-splitting of the itinerant electrons due to the exchange with the localized magnetic moments, the second one the excitation (annihilation) of a magnon due to a spin flip process and the third one the spin-conserving scattering of a magnon and an electron from one state to another.It is worth noting that the second term leads to a transfer of both energy and angular momentum (i.e., spin) -since it can change the total number of magnons -while the third term can only transfer energy.For this reason, this term was discarded earlier works [63][64][65], however, our quantitative analysis reveals that the energy transferred by this term can exceed the energy transferred by the term first order in magnon operators.
We complete our Hamiltonian H = Ĥe + Ĥm + Ĥem by considering Ĥe = kνσ ε kνσ ĉ † kνσ ĉkνσ and Ĥm = q ℏω q b † q bq , with ε kνσ = ε kν − 2∆δ σ↑ being the mode and spin dependent electron energies that are calculated from first-principles calculations and ℏω q being the magnon energies.Note that we have absorbed the term zero-th order in magnon variables in Equation (1) in the otherwise spin-independent Ĥe .
Next, we use the Hamiltonian introduced above to construct a quantum kinetic approach for the description of the out-of-equilibrium dynamics of electrons and magnons.We define the rates of energy exchange between both subsystems as where the dot represents temporal derivative and with the electron (f kνσ ) and magnon (n q ) occupation numbers.The equivalence in Equation (3) results from the conservation of total energy.The time derivatives of the occupation numbers can be calculated by applying Fermi's golden rule to the scattering Hamiltonian (1).
To simplify the calculations, we further assume a thermal electron distribution and can hence introduce a single electronic temperature T e that relates to the occupation of electronic states via the Fermi-Dirac distribution.This allows us to apply and also extend (by including terms second order in the bosonic operators) the ideas laid out in Allen's seminal work on electron-phonon interaction [82] to electron-magnon scattering, yielding ṅq = n BE (ω q , T e ) − n q γ q + q ′ (n q + 1)n q ′ n BE (ω q − ω q ′ , T e )+(q ↔ q ′ ) Γ qq ′ , with n BE (ω q , T e ) = [e ℏωq k B Te −1] −1 being the Bose-Einstein distribution evaluated at the electron temperature.The scattering rates are given by with ε F being the Fermi energy.The functions I σσ ′ (T e ) have the property lim Te→0 I σσ ′ (T e ) = 1 and account for the smearing of the Fermi-Dirac distribution at high electron temperatures, similar to what has been derived for electron-phonon scattering [35].The expression for I σσ ′ (T e ) and details of the derivation of Equations ( 4)-( 5) are in the Method Section V A. Note that a comparison with linear spin-wave theory in the framework of the Landau-Lifshitz-Gilbert equation [83] reveals that γ q /ω q = α q can be viewed as a mode-dependent Gilbert damping parameter.Due to the assumption that the electron occupation numbers follow the Fermi-Dirac distribution at all times, the change in electron energy is determined by the change in T e , i.e., Ėe = kνσ ε kνσ (∂f kνσ /∂T e ) Ṫe = C e Ṫe , with the electronic heat capacity C e = kνσ ε kνσ (∂f kνσ /∂T e ).By additionally considering the absorption of a laser pulse with power P (t) by the electrons and a coupling of the electrons to a phonon heat bath as in the 2TM, we finally obtain our out-ofequilibrium magnon dynamics model: Here, T p , C p and G ep are the phonon temperature and heat capacity and electron-phonon coupling con-stant, respectively.Note that we do not consider direct magnon-phonon coupling, which has been shown to be a reasonable approximation for 3d ferromagnets [43,44].We would like to point out that the nonthermal magnon occupations n q can be translated to mode-specific temperatures via the Bose-Einstein distribution, T q := ℏω q /(k B ln(n −1 q + 1)).Based on this -and in distinction from the 3TM -we term the framework provided by Equations ( 6)-( 8) the (N+2)-temperature model ((N+2)TM).Below, we reveal by solving these coupled equations numerically that they provide a viable framework to describe laser-induced ultrafast magnetization dynamics and the generation of non-thermal magnons, going beyond the well-established 3TM.
Before doing so, we want to shortly discuss the relation between the (N+2)TM introduced here and the 3TM.Albeit their phenomenological nature, the 2TM (T e and T p ) and the 3TM (T e , T p and T m ) have been successfully applied to explain a plethora of phenomena [84], perhaps most prominently by Beaurepaire et al. to describe the ultrafast demagnetization of Ni [1].Allen [82] and Manchon et al. [78] demonstrated that the 2TM and the 3TM can be derived from a microscopic out-of-equilibrium approach similar to the one used here.By assuming instantaneous relaxation of the magnon occupation numbers to the Bose-Einstein distribution with a single magnon temperature T m , our (N+2)TM reduces to the 3TM (in absence of magnon-phonon coupling), with the magnon heat capacity C m = q C q = q ℏω q (∂n q /∂T m ) and the electron-magnon coupling constant Details of the derivation are found in Method Section V B. The above expression goes beyond what was derived in Ref. [78] by including terms second order in magnon variables and allows us to determine the electron-magnon coupling fully based on ab initio parameters.We would like to point out that it can be extended further by going to higher order in the magnon variables.

A. Magnon lifetimes and Gilbert damping
We apply the (N+2)TM model defined by Equations ( 6)-( 8) to bcc iron.To obtain a full solution of the outof-equilibrium dynamics, it is required to calculate material specific quantities.First, we estimate ∆ ≈ 0.75 eV from the band structure and with that we compute the given as color code, shown along high-symmetry lines of the BZ.The lifetimes are due to the first-order contribution to the electron-magnon scattering.
quantities γ q , Γ qq ′ and I σσ ′ (T e ), both using the fullpotential linear augmented plane wave code ELK [85] (details can be found in the Method Section V C).For bcc iron it turns out that I σσ ′ (T e ) only scales weakly with temperature and hence we use the low temperature limit I σσ ′ (T e ) = 1 hereinafter.The parameters governing the magnon energies ℏω q = S(2d ) were taken from earlier works: the exchange constants J ij are from first-principles calculations [86] and the magneto-crystalline anisotropy energy d = 6.97 µeV per atom is from experiments [87].Further, we used the saturation moment µ s = 2.2 µ B and spin S = 2.2/2.Based on these parameters and the formulas derived above, we get C m = 5.720 × 10 4 Jm −3 K −1 and G em = 6.796 × 10 17 Wm −3 K −1 at T m = 300 K. Notably, the term first order in magnon variables leads to a contribution to G em that is one order of magnitude smaller than the second-order term.We further use the roomtemperature values C e = 1.013 × 10 5 Jm −3 K −1 , C p = 3.177 × 10 6 Jm −3 K −1 and G ep = 1.051 × 10 18 Wm −3 K −1 that were obtained in Refs.[35,37] from first-principles calculations.The influence of a temperature dependent electronic heat capacity C e on the demagnetization is discussed in the Supporting Information.
Both ω q and the inverse of γ q , i.e., the lifetime of magnons due to the contribution to electron-magnon scattering linear in the magnon variables, are shown in Figure 1 along high-symmetry lines of the Brillouin zone (BZ).It can readily be observed that the lifetimes of high-frequency magnons are drastically reduced as compared to the low energy ones.The q-dependent lifetimes give rise to mode-specific Gilbert damping α q (= ω q /γ q ).Our finding of mode-dependent Gilbert damping is consistent with experiments [88] and also with a recent fieldtheory derivation [89].The computed α q values, shown in Method Section V C, range between 1.5 × 10 −3 and 1.08 × 10 −2 .These values are close to the experimentally obtained ones (via FMR measurements) for Fe ranging from 1.9 × 10 −3 to 7.2 × 10 −3 [90][91][92][93][94][95], however with a − nq)/(N S − q n init q ), with n init q = n BE (ωq, 300 K) being the occupation number before the laser pulse.(c) Demagnetization max(|∆M/M0|) versus laser fluence computed for a ferromagnetic layer with a thickness of 20 nm.The dotted line serves as a guide to the eye.somewhat larger variation with q as compared to what was reported in Ref. [83].
We note that the q-dependent Gilbert damping goes beyond the conventional LLG description which assumes one single damping parameter for all spin dynamics.Moreover, a further distinction between the current theory and the LLG framework is that, in the latter, there is a single damping term that governs both the energy and angular momentum transfer [96], whereas the current theory has two terms [see Equation (1)], one that transfers energy and angular momentum and one that transfers only energy.As shown in the following, this 2 nd term is found to be important for non-thermal magnon generation.

B. Ultrafast dynamics
Based on the above-given parameters, we calculate the coupled out-of-equilibrium magnon, electron, and phonon dynamics induced by a Gaussian laser pulse 2] with A = 9.619 × 10 7 Jm −3 and ζ = 60 fs for N = 20 3 magnon modes.Note that this value of A translates to an absorbed fluence of 0.19 mJ/cm 2 for a ferromagnetic layer with thickness of 20 nm, which is a typical thickness in ultrafast demagnetization experiments [1].
Figure 2(a) depicts the time evolution of electron, phonon and average magnon temperature -together with the temperature range of all magnon temperatures -calculated using the (N+2)TM.The electron temperature reaches a maximum of 685 K at around 52 fs after the maximum of the laser pulse (located at t = 0) and converges to the phonon temperature in less than 1.5 ps.The maximum of the average magnon temperature of 520 K is reached only slightly after the electronic one at around 136 fs, followed by a convergence to the electronic and phononic temperature to a final temperature of around 329 K at 3 ps, in agreement with what can be estimated from the energy supplied by the laser pulse and the individual heat capacities via ∆T = A/(C m + C e + C p ) = 28.8K. Notably, the magnon temperatures still cover a range of around 50 K at this point in time.Our results clearly demonstrate the shortcomings of the conventional 3TM (shown as dotted lines): While the initial increase of temperatures is comparable to the (N+2)TM, magnon thermalization happens much faster in the 3TM.
In Figure 2(b), we show the laser-induced change in magnetization (associated with the localized magnetic moments) due to the creation of additional magnons.We observe ultrafast transversal demagnetization of around one percent in less than 300 fs, demonstrating that the timescales obtained by our ab initio based calculations are in reasonable agreement with experimental measurements (see, e.g., [15,[97][98][99]).Notably, the minimum of the magnetization and the maximum in the average magnon temperature computed by the (N+2)TM are at different points in time.Also, the drop in the (localized) magnetization is much less pronounced than expected from the increase in average temperature: in thermal equilibrium, a temperature increase from 300 K to above 500 K approximately leads a demagnetization of 20% for iron [100].These observations clearly demonstrate the shortcomings of the 3TM -where a thermal magnon distribution at all times is assumed -and underline the importance of treating the full, non-thermal magnon distribution in the ultrafast regime.
Figure 2(c) depicts the maximum of the demagnetization versus laser fluence for an iron layer of 20 nm.We find a nonlinear dependence, which is a result of the nonlinearity of our (N+2)TM, and a substantial demagneti- zation of around ten percent at 0.95 mJcm −2 .We note that for high fluences, higher-order magnon-magnon scattering terms that are not included in the current model could start to play a role.
The obtained amount of demagnetization and the magnetization decay time (below 200 fs) for this fluence are comparable with experiments, which suggests that ultrafast magnon excitation [15][16][17] provides a viable mechanism for ultrafast laser-induced demagnetization.It is also consistent with time-resolved extreme ultraviolet magneto-optical and photoemission investigations that detected magnon excitations during ultrafast demagnetization of elemental ferromagnetic samples [18,101].
For a more precise examination of the predictions of the (N+2)TM, we compare the calculated timedependent demagnetization with experimental data for Fe in Figure 3.The experimental data were measured by Carpene et al. [15] on a 7-nm thin film, using the time-resolved magneto-optical Kerr effect for two different pump laser fluences of 1.5 mJcm −2 and 3 mJcm −2 .In the calculations we used an absorbed laser fluence that is about five times lower, as the exact value of the absorbed fluence in experiments is difficult to estimate (due to influence of optical losses, sample reflection, etc).Specifically, in the simulations we used absorbed laser energies of 433 Jcm −3 and 693 Jcm −3 in a 7-nm Fe film.Figure 3 exemplifies that not only the amount of demagnetization but also the full time dependence of the demagnetization predicted by the (N+2)TM is in remarkable agreement with experiments.

C. Non-thermal magnon dynamics
Next, we analyze the non-thermal magnon dynamics in more detail in Figure 4. There, we show the magnon temperatures versus frequency (a) and along high-symmetry lines of the BZ (b) at different points in time.The laser pulse primarily heats up high energy magnons, while the temperature of low energy magnons barely changes and even decreases slightly in the vicinity of the Γ point (the temperatures drop by up to around 2.5 K).This surprising observation is caused by a redistribution of magnons from this region to other parts of the BZ due to the term second order in the magnon operators in Equation (1); the effective second order scattering rate γ q := q ′ Γ qq ′ is negative for low magnon frequencies (more details can be found in the Method Section V C).It is also observed that although the magnon temperatures reached after the laser pulse are generally higher at higher frequencies, however, there is not necessarily a monotonous increase of temperature with frequency at all times: e.g., at 100 fs after the laser pulse [Figure 4(b)], the temperatures at the points H, N, and P is higher than in between these points.Notably, the position of the maximum magnon temperature in the BZ also varies with time.

D. Discussion
Different physical mechanisms have been proposed for ultrafast demagnetization in elemental 3d ferromagnets [9,11,15,18].The preeminent mechanisms are Elliott-Yafet (EY) electron-phonon spin-flip scattering [11,13] and ultrafast magnon generation [15].In the former, a Stoner-type picture is used to model the longitudinal reduction of the atomic moment due to electron-phonon spin-flip scattering, whereas the latter is based on lengthconserving transverse spin-wave excitations.Experimental indications of electron-phonon scattering [38,39] as well as of electron-magnon scattering have been reported [18,101].
The strength of the different demagnetization channels is an important issue in the on-going discussion on the dominant origin of ultrafast demagnetization [7].Ab initio calculated quantities such as the EY spin-flip probability are essential to achieve reliable estimates [102][103][104].Griepe and Atxitia [14] recently employed the microscopic 3TM [11] and obtained quantitative agreement with measured demagnetizations for the elemental 3d ferromagnets.They compared the fitted EY spinflip probability α sf with ab initio calculated values [104] and found these to be in good agreement, in support of an electron-phonon mechanism of ultrafast demagnetization.A drawback of their employed approach is however that only magnetization reducing spin flips are included.EY spin flips that increase the magnetization are also possible and, including these would lead to a significantly smaller demagnetization amplitude [104].This in turn would question again what amount of demagnetization is precisely due to EY electron-phonon spin flip scattering.Conversely, in our non-thermal magnon approach we employ ab initio calculated quantities without fit parameter.We find that the ab initio predicted ultrafast demagnetization agrees accurately with experiments, which provides a strong support for the prominence of the non-thermal magnon channel to the ultrafast demagnetization process.

IV. CONCLUSIONS
We have developed an ab initio parameterized quantum kinetic approach to study the laser-induced generation of magnons due to electron-magnon scattering, which we applied to iron.Our results clearly demonstrate that on ultrafast timescales the magnon distribution is non-thermal and that henceforth the simple relation between magnetization and temperature via the M (T ) curves computed at equilibrium does not hold: since predominantly high-energy magnons are excited the energy transferred from the laser-excited electrons creates relatively few magnons and hence the demagnetization (proportional to the total number of magnons) is much less pronounced than expected from the increase of the average magnon temperature.Notably, the number of magnons actually decreases near the center of the Brillouin zone, which is due to the scattering from low to high energy magnons by a previously neglected scattering term that can transfer energy but not angular momentum.This term, which is not included in LLG simulations, is a crucial quantity for out-of-equilibrium magnon dynamics.
Our ab initio-based calculations of the induced demagnetization in iron furthermore provide strong evidence that non-thermal magnons are excited fast and lead to a sizable demagnetization within 200 fs.The resulting time-dependent demagnetization agrees remarkably well with experiments, which establishes the relevance of magnon excitations for the process of ultrafast optically induced demagnetization.

V. METHOD A. Derivation of electron-magnon scattering rates
In this Method Section we derive the (N+2)TM for the description of non-thermal magnons from a microscopic Hamiltonian for electron-magnon scattering.We start with a local sp − d model Hamiltonian, with J sp−d being the sp − d volume interaction energy, ŝitin = σ being the spin operators of itinerant (s and p) electrons and S loc i being the localized (d) spins located at r i .For now, we treat the latter as classical vectors.The expectation value for a given spin wave function Ψ(r) is given by Here, we have introduced S ± i = S x i ± iS y i .Next, we perform a plane wave expansion of the wave functions (for a single band of itinerant electrons), and a Holstein-Primakoff transformation of the localized spins, together with introducing the Fourier transform of the magnon amplitudes Insertion of ( 15)-( 17) into ( 14) and keeping terms up to second order in magnon variables, we get For multiple itinerant bands and in second quantization we obtain where we have introduced ∆ = J sp−d SN V .Note that due to the plane wave ansatz we have implicitly assumed that the itinerant electrons are completely delocalized and interband scattering (from ν to ν ′ ̸ = ν) fully contributes to the electron-magnon scattering.
Next, we use Fermi's golden rule to get the change of the magnon occupation number n q = ⟨ b † q bq ⟩.Fermi's golden rule computes the probability W (i → f ) for a small perturbation term in the Hamiltonian, Ĥ′ (in our specific case, Ĥem ) via where |i⟩ and |f ⟩ denote the initial and final state, respectively.We start with the term first order in the magnon variables, ṅ with f kνσ = ⟨ĉ † kνσ ĉkνσ ⟩ and ε kνσ and ℏω q being the eigenenergies of electrons and magnons, respectively.Hereinafter, we make the assumption that due to the fast equilibration processes for electrons, they always follow the Fermi-Dirac distribution, f FD (ε kνσ , T e ) = [e (ε kνσ −εF)/kBTe + 1] −1 , with a single electron temperature T e .Before we continue we need the following relation, with n BE (ω q , T e ) = [e ℏωq k B Te − 1] −1 being the Bose-Einstein distribution evaluated at the electron temperature.Now we can simplify Equation (22), yielding With γ q being the linewidth -i.e., the inverse lifetime -of the magnon due to the first order contribution to electronmagnon scattering.Following the ideas laid out by Allen [82] and Maldonado et al. [35], it can be computed as with ε F being the Fermi energy, the spin-dependent density of states is g σ (ε) = kν δ(ε − ε kνσ ) and the thermal correction factor given by It is obvious that lim Te→0 I σσ ′ (T e ) = 1.Note that we have used that the energy scale of magnons is much smaller than the one of electrons, i.e., that ℏω q ≪ ε, ε ′ .The contribution of the term second order in magnon variables to the occupation number can be calculated analogous and reads (n q + 1)n q ′ n BE (ω q − ω q ′ , T e ) + q ↔ q ′ Γ qq ′ (T e ) (37) with B. Derivation of the three temperature model In what follows, it is demonstrated that the three temperature model (3TM) can be obtained from the (N+2)temperature model derived in the main text, ṅq = n BE (ω q , T e ) − n q γ q + q ′ (n q + 1)n q ′ n BE (ω q − ω q ′ , T e ) + (q ↔ q ′ ) Γ qq ′ , ( by assuming instantaneous relaxation of the magnon occupation numbers to the Bose-Einstein distribution with a single magnon temperature T m , i.e., n q = n BE (ω q , T m ).For the sake of readability we rewrite n BE (ω q , T m ) = n q (T m ).
We start with the first order scattering term: Here we have introduced the mode-dependent magnon heat capacity C q = ℏω q ∂nq(Tm) ∂T .In order to calculate the scattering term second order in the magnon variables, we first introduce the following relation Now we calculate ṅ(2) q = q ′ (n q (T m ) + 1)n q ′ (T m )n q−q ′ (T e ) + (q ↔ q ′ ) Γ qq ′ Using the expressions for ṅ(1) q and ṅ(2) q , the change in total energy of the magnons can then be calculated as With that, the (N+2)TM transforms into the 3TM (in the absence of magnon-phonon coupling), which is given by

C. Ab initio calculations
To obtain a full solution of the (N+2)TM, it is necessary to compute the material specific quantities ∆, γ q , Γ qq ′ and I σσ (T e ).For this purpose, we use the full-potential linear augmented plane wave code ELK [85].
As a first step, we determine the coupling parameter ∆ of the sp − d model, which sets the general scale of the electron-magnon scattering.As shown in the main text, the first term (zeroth order in magnon variables) in the electron-magnon scattering Hamiltonian reads Ĥ(0) em = −∆ kν (ĉ † kν↑ ĉkν↑ − ĉ † kν↓ ĉkν↓ ), with ν ∈ {s, p}.Based on this, ∆ can be estimated from the projected density of states (DOS), since it is one half of the spin-dependent energy splitting of the s-and p-bands.In general, this splitting may vary for different electronic states.This is not accounted for in the model used here, where instead a single parameter is used to model the spin splitting.We find, however, that for bcc iron this is justified, since the shift in both s-and p-bands around the Fermi energy -the relevant region for electron-magnon scattering -between spin up and down states is approximately constant with a value of ∆ ≈ 0.75 eV, see left panel of Figure 5. Now we calculate the first and second order scattering rates using the formulas derived above, The calculation of both quantities requires a spin-dependent summation over the Fermi surface, analogous to what was done in Ref. [103] for the evaluation of the spin-dependent Eliashberg function for electron-phonon scattering.As in Ref. [103] we use a Gaussian broadening of the Dirac delta distributions by 0.03 eV.Also, since we only include the contribution of s-and p-states (indicated by ν, ν ′ ) to the scattering, we have to project the Kohn-Sham states (indicated by n, n ′ ) onto the spherical harmonics with P nν kσ being projector functions.The functions I σσ ′ (T e ) describe corrections to the scattering rate at high electron temperatures and are given by with g σ (ε) = kν δ(ε − ε kνσ ) = kν n P nν kσ δ(ε − ε knσ ) being the cumulative DOS of both s-and p-states.We find that they increase monotonously with the electron temperature (see right panel of Figure 5).However, even for temperature up to 2000 K, the I σσ ′ (T e ) functions are below two.Hence, we concluded that the approximation I σσ ′ = 1 is reasonable for the laser fluences -heating the electrons up to around 700 K -considered in the main text.
Figure 6 depicts the numerically calculated scattering rates using I σσ ′ = 1 and ∆ = 0.75 eV as obtained above.In the left panel, we show the scattering rate γ q that is first order in the magnon variables through color code on the magnon dispersion.It is strictly positive and tends to increase with magnon frequency.The right panel shows the effective scattering rate γ (2) q = q ′ Γ qq ′ due to the scattering term second order in magnon variables.Notably, this quantity is negative for low frequencies and positive for high frequencies, indicating that it leads to a depopulation of magnons at low energies due a scattering from low to high energies (the total magnon number is kept constant).In general, the values of the effective second order scattering rate are comparable to the one first order in magnon variables.They are, however, distributed differently: e.g., for magnons close to the Γ point the second order scattering rate is by far the dominating one.This is the reason why, as demonstrated in the main text, a laser pulse can in fact lead to a cooling of low energy magnons, i.e., to a decrease of their occupation numbers.
Lastly, we show in Figure 7 the ab initio computed mode-dependent Gilbert damping, α q = ω q /γ q .Interestingly, the Gilbert damping α q is large (∼ 0.01) at the BZ center and at the high-symmetry points H, N and P at the BZ edge.There is also a noticeable directional anisotropy in the Gilbert damping for modes along Γ−H and Γ−P.We emphasize that the Gilbert damping is here due to the electron-magnon scattering term that is first order in the magnon variables.Other scattering mechanisms as phonon-magnon scattering could contribute further to the mode-specific Gilbert damping.

Figure 2 .
Figure 2. Laser-induced ultrafast non-equilibrium dynamics of iron calculated from an ab initio parameterized model.(a) Temporal evolution of electron temperature Te, phonon temperature Tp and average magnon temperature ⟨Tq⟩ = 1/N q Tq obtained by the (N+2)TM (solid lines).The blue shaded region indicates the temperature range within which all magnon temperatures are contained.Dashed lines show the results of the 3TM solved with ab initio calculated input parameters.(b) Relative change of total magnetization of the localized magnetic moments ∆M/M0 = q (n init q

Figure 3 .
Figure 3.Comparison between experiment and the (N+2)TM theory for ultrafast demagnetization in iron.The experimental data (symbols) are those of Carpene et al. [15] and the solid lines are calculated from the ab initio parameterized (N+2)TM.

Figure 4 .
Figure 4. Magnon temperatures of iron during ultrafast laser excitation at different points in time (w.r.t. the maximum of the laser pulse) calculated from the ab initio parameterized (N+2)TM.(a) Magnon temperatures (dots) versus frequency.The solid line indicates the electron temperature.(b) Magnon dispersion and their temperatures, depicted by the color code, shown along high-symmetry lines of the BZ.

Figure 5 .
Figure 5. Left: Projected spin-polarized DOS for bcc iron.Spin-minority density is shown by positive values, spin-majority density by negative values.The exchange splitting is 2∆ ≈ 1.5 eV in a large interval around the Fermi energy and for both sand p-states.Right: Thermal correction factors I σσ ′ versus electron temperature Te calculated from the projected DOS.

Figure 6 .Figure 7 .
Figure 6.Magnon dispersion of bcc iron along high-symmetry lines of the Brillouin zone.The color coding describes (left) the scattering rates γq due to the electron-magnon scattering term first order in magnon variables γq and (right) the effective scattering rate γ (2) q= q ′ Γ qq ′ due to the term second order in magnon variables, calculated with I σσ ′ = 1 and ∆ = 0.75 eV.