Accounting for Quantum Effects in Atomistic Spin Dynamics

Atomistic spin dynamics (ASD) is a standard tool to model the magnetization dynamics of a variety of materials. The fundamental dynamical model underlying ASD is entirely classical. In this paper, we present two approaches to effectively incorporate quantum effects into ASD simulations, thus enhancing their low temperature predictions. The first allows to simulate the magnetic behavior of a quantum spin system by solving the equations of motion of a classical spin system at an effective temperature relative to the critical temperature. This effective temperature is determined a priori from the microscopic properties of the system. The second approach is based on a \semi model where classical spins interact with an environment with a quantum-like power spectrum. The parameters that characterize this model can be calculated ab initio or extracted from experiments. This semi-classical model quantitatively reproduces the absolute temperature behavior of a magnetic system, thus accounting for the quantum mechanical aspects of its dynamics, even at low temperature. The methods presented here can be readily implemented in current ASD simulations with no additional complexity cost.


I. INTRODUCTION
Magnetic materials modeling is crucial for the understanding of groundbreaking experiments [1], as well as for a variety of new technologies [2,3].The development of ultra-fast laser pulses led to the discovery of ultra-fast demagnetization induced by thermal fluctuations [1] and magnetization reversal [4,5].These discoveries paved the way for the development of new technologies such as heat assisted magnetic recording [6].Furthermore, thermal effects play a crucial role in helicity-dependent magnetization switching [7,8].
Atomistic spin dynamics (ASD) simulations are among the most effective tools for studying magnetization dynamics in materials with thermal fluctuations [9][10][11][12].They model the atoms in the materials as localized classical magnetic moments interacting through an effective Hamiltonian.This Hamiltonian can be fully parameterized with ab initio calculations and numerically solved, including Markovian thermal fluctuations, to obtain the spin dynamics.ASD simulations have been successful in several situations e.g., in predicting the experimentally observed heat-assisted magnetization reversal [4,5], modeling [13,14] the domain wall dynamics induced by thermal gradients [15,16] and the current-induced skyrmion dynamics [17] observed in [18].Despite their versatility and predictive capability, ASD simulations, being predominantly classical in nature, exhibit inaccuracies in modeling low-temperature demagnetization driven by quantum excitations like spin waves [19].Accurately accounting for quantum fluctuations would require simulating the full quantum dynamics.The dynamics of coupled quantum spins can be studied with different numerical methods, e.g.quantum Monte Carlo simulations [20] and tensor-network methods [21,22].However, the quantum nature of the problem limits the system size while restricting to high temperatures and short time scales [23].In fact, given the current hardware limitations, the dynamics of more than a few dozens of quantum spins is already inaccessible.In contrast, one of the appealing characteristics of ASD simulations is their capability to model the dynamics of millions of coupled classical spins [11,12,24].This makes them suitable for studying systems where the translational symmetries are broken, e.g. by defects or local perturbations.
To improve the accuracy of ASD simulations in the low temperature regime, in Ref. [25] the authors employ a heuristic approach to reproduce the experimental temperature dependence of a material's magnetization.In this approach, they postulate that the experimental temperature T exp and the atomistic spin dynamics simulation temperature T sim are related by a power-law, as τ sim = τ α exp , where τ sim/exp = T sim/exp /T sim/exp c and T sim/exp c is the Curie temperature of the material.The required rescaling exponent α was obtained by comparing the results of the ASD simulations with experimental measurements.Importantly the temperature rescaling not only reproduces the systems equilibrium, but it also captures the behavior of the ultra-fast magnetization dynamics in pump-probe experiments [1].While this approach seems effective across various situations, it is purely heuristic, lacking a microscopic understanding of its underlying mechanisms.
An alternative approach to account for quantum effects in ASD simulations was proposed in Ref. [26], where the authors solve (Markovian) Landau-Lifshitz-Gilbert dynamics, but include stochastic noise in the spin dynamics which has a quantum power spectrum.Thanks to this approach, it is possible to reliably model the low-temperature behavior of the magnetization of materials.However, the high temperature behavior was not fully captured, i.e.Curie temperatures are somewhat overestimated.
Building on these previous works, in this paper we first provide an analytical expression and a physical explanation for the power-law rescaling of Ref. [25].This is done in Sec.III, where we compare the classical spin model to its quantum mechanical counterpart in the mean field framework.By means of this comparison, we derive rescaling power-laws, which turn out to be in good agreement with the ones found in Ref. [25].Strikingly, we do not need preliminary ASD simulations to calibrate the rescaling exponent.Our predictive method can be readily applied to ASD simulations to account for dynamical quantum effects at no additional computational cost.
Secondly, in Sec.IV of this paper we follow the quantum noise spectrum approach of [26] and improve it to better capture the high temperature regime and predict critical temperatures.To achieve this we directly solve the steady-state of a classical spin system interacting with a quantized non-Markovian bath.This approach does not require any rescaling and quantitatively reproduces the low temperature behavior of the magnetization as well as its critical temperature T c .This model can be fully parameterized by means of ab initio calculations or experiments and constitutes a highly promising candidate for improving the accuracy of ASD simulations by accounting for quantum effects.
Before reporting on these two methods to include quantum effects in ASD, we here introduce the general setting.

II. GENERAL SETTING
We consider the well-known Heisenberg model [27], which is described by the Hamiltonian where J is the exchange parameter which can be calculated ab initio [28], γ is the electron gyromagnetic ratio, H ext is an external magnetic field pointing in the z direction, the sum over i runs over all the lattice sites, and ⟨i, j⟩ over nearest-neighbor pairs.For a periodic system made of spins of length S 0 = nℏ/2 with n a natural number, the mean field approximation (MFA) self-consistent equation for the expectation value of the normalized spins [29] [30].Here, β = 1/k B T with k B the Boltzmann constant and T the temperature, F cl/qu (x) is the Langevin function F cl (x) = coth(x)−1/x for the classical Heisenberg model, and the Brillouin function F qu (x) = {(n + 1) coth [(n + 1)x/n] − coth (x/n)}/n for the quantum Heisenberg model with n = 2S 0 /ℏ [30].We show the normalized spin magnetization sz for a spin system with S0 = ℏ/2 (see Eq. ( 2)) evaluated with the MFA as function of the relative temperature τ cl/qu = T /T cl/qu c for classical spin vectors (blue dotted line) and quantum spin operators (red dashed line).Note that the only effect of J and z is to determine Tc.Since here we calculate the magnetization as a function of the relative temperature τ cl/qu = T /T cl/qu c , J and z do not need to be specified.We also plot the classical (blue solid line) and quantum (red solid line) single-spin sz curves in a field H 0 eff = zJS0/γ .
Furthermore, H eff (s z ) = H ext + zJS 0 s z /γ is the effective magnetic field felt by each spin in the lattice, with z being the coordination number.
It is insightful to express the magnetization curves s = zJ⟨S 2 ⟩ cl/qu /(3k B ), where ⟨S 2 ⟩ cl = S 2 0 for the classical case and ⟨S 2 ⟩ qu = S 0 (S 0 +ℏ) for the quantum one.These critical temperatures are typically overestimated since the MFA does not account for the system's fluctuations [11].Note that, in terms of the relative temperatures τ cl/qu , the expectation value of the normalized spin takes the self-consistent form Therefore, in the MFA for the Heisenberg model the dependence of s cl/qu z is reduced to two free parameters, namely the relative temperature τ cl/qu and the spin length S 0 .
The deviation of classical ASD simulations from experiments is most prominent at low temperatures [25], particularly so for magnetic materials with atoms carrying a small magnetic moment.Indeed, it is in this regime that quantum contributions play an important role in spin systems causing discrepancies between the classical and quantum results [31][32][33].To see how these deviations arise, we can consider the spins at low temperature to be approximately in the ground state.The spins then fully align along the z direction (s z ≃ 1) and the effective field becomes independent of s z , H T ≃0 eff ≃ zJS 0 /γ.Thus, in this limit, the dependence of the right hand side of Eq. ( 2) on s z can be neglected.This makes the expectation value of each spin independent of the expectation values of the other spins in the system, reducing the multi-spin problem to a single-spin problem in the low-T limit (see Fig. ( 1)).While a single classical spin has a continuous energy spectrum, the quantum counterpart does not.This implies that, already at the smallest finite temperature, the classical spin can experience infinitesimally excited states such that the expectation value s z deviates from the ground state s z = 1.On the other hand, the expectation value of the quantum spin remains in the ground state s z = 1 as long as the thermal energy k B T is smaller than the energy gap between the ground state and the first excited state, ∆E = γH T =0 eff ℏ.In Fig. 1, we show the comparisons between the expectation value of a single-spin in an effective field H 0 eff and the expectation value of a multi-spin system obtained with the MFA as functions of the relative temperature τ cl/qu .This comparison is shown for both quantum and classical cases.In the limit of zero temperature, the single-spin calculations (solid lines) closely resemble those of the multi-spin case obtained by means of the MFA (dotted and dashed lines).Note also that, in both the single-spin as the multi-spin case, while the classical spin (cyan solid and dotted lines) has a finite gradient at T = 0, the quantum spin (red solid and dashed lines) is characterized by a plateau which only drops once τ qu ≳ 0.4.That is, the classical spin is susceptible to arbitrarily small temperatures while, for the quantum case, there is a minimum amount of thermal energy needed to excite the system.
We conclude that the difference between the quantum and classical MFA can be tracked down to the singlespin qualitative behaviour in an external field.This observation is key in designing the approach described in Sec.IV.There the assumption is used that it is possible to account for quantum effects by just coupling each individual classical spin of the system to an environment with a quantum-like power spectrum.

III. POWER LAW TEMPERATURE RESCALING IN THE MFA
The experimental behaviour of the magnetization at low temperature is known to follow the Bloch's law M (T ) ≃ M 0 (1 − (τ exp ) 3/2 ).On the other hand, classical ASD simulations, in the same temperature range, are known to produce a linear magnetization dependence on the temperature M (T ) ≃ M 0 (1 − γτ sim ).To overcome this discrepancy, in Ref. [25], Evans et.al. introduced a heuristic temperature rescaling to ASD simulations to successfully reproduce the experimental temperature dependence of the magnetization.This method proves effective in several contexts and it is now a staple in the ASD community's toolkit.Nevertheless, it lacks a detailed physical interpretation which we will provide below.This temperature mapping from the classical to the quantum behavior was done by means of the rescaling for a power α ∈ R that was fixed by fitting the simulations to experimental data.The estimation of α required the comparison of classical ASD simulations with the experimental M (T ) curve or alternatively the Kuz'min law.The latter is a phenomenological equation able to accurately reproduce the experimental M (T ) while also complying with the Bloch's law [34].This same rescaling was found successfully matching the experimentally observed ultra-fast demagnetization in nickel [1,25].Despite the successes of this method two main questions persist: what is the microscopic origin of such rescaling?Is there a strategy that allows to calculate the rescaling exponents that does not require prior knowledge of experimental results?Here, we provide a physical interpretation for the mechanism behind the success of such rescaling.We propose a method for calculating the rescaling exponent α without recurring to any prior experimental knowledge about the material's behaviour.Additionally, we give an analytical expression to evaluate α suggesting its sole dependence on the spin length S 0 .Our method can be readily used in ASD simulations to improve the low-temperature description without additional computational cost.
To understand the origin of this rescaling, we first note that the quantum-classical difference is already apparent in the MFA expression (2).The key idea in this section, is to fit the curves of the classical MFA solution s cl z to the curves of the quantum MFA solution s qu z , assuming a temperature rescaling of the type shown in Eq. ( 3), that is τ cl = (τ qu ) α for various exponents α.When expressing the temperature in terms of τ = T /T c , the only free parameter left in Eq. ( 2) is the spin length S 0 .Therefore, to find the best fitting exponent α for a given material will require an appropriate choice of S 0 .To do so, we write S 0 = nℏ/2 with n = 1, 2, . . ., and note that a spin of length ℏ/2 carries one Bohr magneton.We then choose the smallest integer n that obeys n ≥ µ[µ B ], where µ[µ B ] is the atomic magnetic moment of the material expressed in units of Bohr magnetons.These magnetic moments are obtained e.g. by means of the augmented spherical wave (ASW) implementation of density functional theory [35].Finally, we determine α in Eq. ( 3) by minimizing the mean squared error where s cl/qu z are evaluated using Eq. ( 2).The temperature rescaling exponent determined by this MFA based method is applicable within a specific temperature range where the relevant microscopic characteristics of the materials remain constant.For instance, in the case of nickel and gadolinium, this method is expected to be effective within their temperature ranges [0, T c ].However, for materials like cobalt, which exhibit a coexistence of hcp and fcc phases, it becomes challenging to uniquely determine the coordination number z. Indeed we expect our approach for calculating the rescaling exponent to work as long as the material in question has a single coordination number z across the entire temperature range.Additionally, when the Heisenberg model proves inadequate due to factors such as the presence of magnetic anisotropy, as in the case of iron, this method may require extensions.
To illustrate the results obtained with the procedure discussed above, we look at the case of Ni and Gd whose atomic magnetic moments are µ Ni = 0.66µ B and µ Gd = 7.63µ B respectively, translating into a spin length S Ni 0 = ℏ/2 and S Gd 0 = 8ℏ/2.The coordination numbers for both Ni and Gd are z = 12, which are constant in [0, T c ]. (Note that the numerical value isn't practically needed for calculations in terms of the relative temperature, see Eq. (2).) In Fig. 2, we compare the magnetization curves for Ni and Gd: s qu z (τ qu ) predicted by the quantum Heisenberg model, the classical prediction with rescaled temperature s cl z ((τ qu ) α ), and the Kuz'min phenomenological model [34] which closely matches experiments.We further plot the bare (without rescaling) classical prediction s cl z (τ cl ).As expected from the previous discussion, the standard classical and quantum MFA calculations disagree, particularly at low temperatures.This discrepancy is notably more pronounced for materials with low S 0 values, such as Ni, compared to those with higher spin values like Gd.The higher spin in Gd leads to behavior that more closely aligns with classical expectations [33].
Crucially, comparing in Fig. 2 the quantum MFA and the classical temperature-rescaled magnetizaion we observe a strong agreement between the two, particularly at low temperatures.We further quantify the difference between these two predictions by the global relative mean square error, , and find that it is at most 2.5%, and decreases as S 0 increases as expected [33].We further note that in the low temperature regime the agreement between the rescaled classical magnetization and the experimental data is even more pronounced.This is likely due to the fact that the temperature rescaling adheres to a power law, which accurately reproduces the Bloch's law, a feature not captured by the bare quantum MFA due to the fact that the bare MFA does not accurately account for fluctuations.Repeating the procedure outlined by Eq. ( 4), we obtain the different scaling exponents α as a function of S 0 .The values for different S 0 are shown in Fig. 3 (black diamonds).Remarkably, we find that the values obtained by this numerical minization procedure can be described to very high accuracy by the simple function (green solid line in Fig. 3) Indeed, this extremely simple expression matches the rescaling exponent α(S 0 ) from Eq. ( 4) with less than 0.05% error.Remarkably, our derived rescaling exponent α, predicts a value very similar to the one used in Ref. [25] for nickel (we obtain α Ni = 2.31 vs α Ni = 2.32 found by the authors in Ref. [25]) and gadolinium (we obtain α Gd = 1.38 vs α Gd = 1.28 found by the authors in Ref. [25]).In Ref. [25], finding the exponent α required running ASD simulation for a wide range of temperatures and fitting  3), obtained by minimizing Eq. ( 4), as a function of spin length S0 (black diamonds).An analytical function, Eq. ( 5), that very closely matches the obtained exponents is also plotted (green line).
the results to experimental data.Here, instead, we determined α (Eq.( 5)) by matching the quantum and classical Heisenberg mean-field predictions.Using the analytical expression in Eq. ( 5), we can find α for a specific material in a self-contained fashion, that is, eliminating the need for fitting with experimental data.Including the rescaling described by Eq. ( 5) in ASD simulations, one can theoretically predict equilibrium magnetization and magnetization dynamics without previous reliance on experimental data.
In conclusion, to account for quantum contributions at a given relative experimental temperature τ exp = T exp /T exp c , one can use standard ASD simulations at a temperature τ sim = (τ exp ) α(S0) , where α(S 0 ) is fully determined by Eq. ( 5).

IV. ACCOUNTING FOR ENVIRONMENTAL FLUCTUATIONS
The temperature rescaling strategy just discussed is based on the relative temperature τ .Therefore, to predict the magnetization as a function of temperature with real units, it requires prior knowledge of the experimental critical temperature T exp c of the material.Here, we propose a second approach, that predicts the magnetization against the absolute temperature with no need for prior knowledge of T exp c .This new method is fundamentally different from the rescaling discussed in Sec.III, and is based on a physically motivated microscopic model where the spins are additionally strongly coupled to a quantum environment.This model can further be fully parameterized through ab initio calculations.As we will see, this approach accurately reproduces the temperature dependence of nickel's magnetization, including lowtemperature behavior and a close approximation of the experimental critical temperature.
In studying magnetic materials, it is essential to consider the interaction between atomic magnetic moments (spins) and non-magnetic degrees of freedom like phonons.To incorporate these degrees of freedom, we linearly couple the spins in the Heisenberg model H H (described by Eq. ( 1)) to independent bosonic baths as in Ref. [36], where describes a continuous frequency reservoir at each lattice site j with X (j) ω being the positions and momentum operators of the 3D reservoir oscillators with frequency ω.The interaction term between the magnetic system and the reservoir is where C ω determines the coupling between the reservoir oscillators and the spins S (j) at lattice position j.This approach serves as a microscopic model for the Landau-Lifshitz-Gilbert (LLG) equation, including inertial terms [36].We assume isotropic coupling of the spins to these external modes i.e.C ω is a scalar, not a tensor [37], and that the spin-phonon coupling C ω is the same for all spins (i.e. it does not depend on the lattice index j).From the theory of open systems, the effects of the reservoir are summarized by the spectral density J (ω) = C 2 ω /(2ω).The spectral density has a well-defined microscopic meaning and, in the case of phonons, is directly linked to the phonon density of state [38].The phonon density of state can further be obtained from first principle calculations [39].For our simulations on nickel, we parameterize the spectral density as a Lorentzian ), which we fit to the experimental phonon density of states reported in Ref. [40].This sets a peak frequency at ω 0 = 6 THz and its broadening at Γ = 8 THz.The coupling strength is set by the dimension-free Gilbert damping η [36], taken here as η = 0.01, in line with values typically reported in the literature [11,41].
The full quantum treatment of the Heisenberg model is unfeasible [42], even more so in the case of an Heisenberg model strongly coupled to a bosonic bath.Our main aim here is to develop a method that is capable to capture relevant quantum effects of the magnetization dynamics and equilibrium properties, without dramatically increasing the computational complexity of standard spin vectors ASD methods.Therefore, in line with standard ASD models, we will assume that the spins are classical and furthermore they interact with a bath of harmonic oscillators.Key, however, is that we will impose the environment oscillators obey a quantum-like power spectrum [19,26].This is reminiscent of the techniques employed in Ref. [43] for the exact modeling of quantum Brownian motion via a semi-classical model.In Ref. [26], a quantum environmental power spectrum with J (ω) = η G ℏω, where η G is the Gilbert damping, was used, which eliminates zero point fluctuations of the environment.Here, we employ the same strategy but where J (ω) is now the Lorentzian parametrized for Ni as described above.
To study the effects of this quantum environment on the equilibrium properties of the magnetization, we extend the MFA to account for the interaction with the environment (EMFA).We perform the usual decomposition of the spins as S (j) = ⟨S (j) ⟩ + (S (j) − ⟨S (j) ⟩), where the first term represent the ferromagnetic order and the second term captures the fluctuations.Inserting this expression in the total Hamiltonian of Eq. ( 6) and neglecting the quadratic terms in the fluctuations we obtain where the effective magnetic field generated by the spin S (j) 's neighbor and the harmonic oscillator bath is ω where the sum is over the i neighbours of j.The Gibbs state of the total system is ρ = e −βHH+B /Tr[e −βHH+B ].The first term in Eq. ( 9) is constant and cancels in the Gibbs state.Hence the expectation value of the spin S (j) can now be calculated using the Gibbs distribution as where the trace Tr[•] is performed over the bath and the spin degrees of freedom.Assuming that all spins behave in the same way, i.e. S = s S 0 for all j, we obtain the quantum phononic Environment+spin neighbours Mean Field Approximation (EMFA) equations for the case of an open Heisenberg model as where β = 1/k B T and Eq. ( 11) defines a self-consistent equation for the expectation value of the spins in the system.In summary, to solve the EMFA equation (Eq.( 11)) we proceed as follows (see Fig. 4): (a) We parametrize the spectral density J(ω) in the power spectrum in Eq. ( 8) for the external reservoir (in our case phonons) with experimental data Ab initio or experimental characterization of the material's phonon spectral density Solve numerically the self-consistent EMFA equation for each temperature including the field from neighbor spins Run SpiDy.jl[36,44] with quantum environmental power spectrum to find the dynamical steady-state for a single spin at all temperatures Apply the SCGA according to Ref. [38] or with ab initio calculations [39].
(b) We simulate the dynamics of a single classical spin in a magnetic field H 0 and interacting with a bath with the quantum power spectrum (8) for each temperature using the method described in Ref. [36,44].Solving for a long enough time, we find the steady state of the single spin system which aligns along the direction of H 0 .Thereby, assuming H 0 = H 0 êz , we find the functional form of F EMFA (H 0 /k B T ) = ⟨s z (T )⟩ H0 in Eq. ( 11) describing the steady-state of a single spin in an external magnetic field H 0 and interacting with a phononic environment at temperature T .
(c) We numerically solve the self-consistent equation ⟨s z (T )⟩ = F EMFA (H eff [⟨s z (T )⟩]/k B T ) using the fixed point method for different T , obtaining s EMFA z (T ) ≡ ⟨s z (T )⟩.Note here that the stochastic part of the effective magnetic field (second term of the right-hand-side of Eq. ( 12)) has already been accounted for in the previous step, therefore H eff [⟨s z ⟩] = J γ zS 0 ⟨s z ⟩.The solution of this equation gives the temperature dependence of the magnetization of the material described by the spin in the Heisenberg model which is additionally coupled to a phonon bath with spectral density J (ω).The EMFA method produces a critical temperature, T EMFA c , beyond which the magnetization vanishes s z (T > T EMFA c ) = 0.
Remarkably, for Ni we find that the solution of Eq. ( 11) as a function of the relative temperature τ (relative to critical temperatures) compares exceptionally well with the Kuz'min curve which approximates the experimental curve (gray).Moreover, our goal with the EMFA method is to provide a good estimate for the absolute critical temperature T c .Note, that such prediction is not possible with the heuristic rescaling method discussed in Sec.III.In contrast, the EMFA method we described in this section is suitable to predict critical temperatures.For Ni we find T EMFA c = 745 K, which compares reasonably well with the experimental critical temperature of T exp c = 632 K [34].At this point it is worth remembering that the EMFA method relies on the MFA, which only takes into account first order fluctuations around the mean.But it is possible to go to higher orders.To do so we make use of the so called self-consistent Gaussian approximation (SCGA).The SCGA is a field theoretical technique that allows to account for the fluctuations S (j) − ⟨S (j) ⟩ at the second order self-consistently similarly to the random phase approximation (RPA) [45].Based on Garanin's work [45], Hinzke et al. showed [46], that the effect of including second order fluctuations can be reduced to a linear rescaling of the temperature in the s z (T ) function as obtained by means of the MFA as T → γ SCGA T , where γ SCGA = 0.79.Therefore, to account for higher order fluctuations and reproduce the absolute temperature dependence of the magnetization we perform a further step.
(d) To account also for second order fluctuations of the spin, we rescale the EMFA temperature as reported in Refs.[45,46] as T → γ SCGA T with γ SCGA = 0.79.
In Fig. 5, the circles line shows the temperature dependence of the magnetization for nickel as predicted by our approach (EMFA+SCGA).The plot here uses physical temperature in Kelvin (T [K]) instead of relative temperature (τ cl/qu ).We compare the (EMFA+SCGA) approach with bare EMFA (triangles line) and the phenomenological Kuz'min curve [34] (gray solid line), which is known to match experimental behavior.We further plot the bare classical prediction in absolute temperature s cl z (T ) (squares line).As shown in the figure, both EMFA and EMFA+SCGA give a good prediction of the experimental s z (T ) trend as a function of the absolute temperature highlighting its qualitative discrepancy with the experiment.Differently to Fig. 2, these magnetization curves are evaluated against the physical temperature in units of Kelvin.If the EMFA and EMFA+SCGA curves were shown over their respective relative temperature τ , they would be the same.We note that both the EMFA and EMFA+SCGA accurately reproduce the low-temperature behavior of the experimental magnetization.Furthermore, EMFA+SCGA gives a relatively accurate prediction of the critical temperature.
in K.In particular, the EMFA gives a prediction for the critical temperature of T EMFA c = 745 K while the more accurate EMFA+SCGA provides T EMFA+SCGA c = 590 K which is in relatively good agreement with the experimental critical temperature T exp c = 632 K [34].We further note that the EMFA accurately reproduce the low-temperature behavior of the experimental magnetization which is not captured by the classical MFA.Importantly, this is a direct consequence of the introduction of the environment.Indeed, in the absence of the environment i.e. setting C ω = 0, the EMFA should reduce to the standard MFA, and both classical MFA and MFA+SCGA (without environment nor temperature rescaling) fail to capture the qualitative behavior in this regime.Importantly, in contrast to the method described in Sec.III, no phenomenological temperature rescaling is employed here.
Finally, we compare our approach to other methods.Our equations of motion, obtained using Eq. ( 6) with the power spectrum in Eq. ( 8), resemble the quantum thermostat approach described in Ref. [26].In both approaches, the environmental fluctuations follow a quantum power spectrum where the zero-point fluctuations are removed.The main difference is that, in our approach, we account for an environment with a structured spectral density (that is with non-Markovian dynamics) instead of white noise.The structured spectral density approach is microscopically motivated and is suited to represent the coupling between spins and environmental degrees of freedom, see Refs.[38,40].In addition, we solve the equations of motion for a single spin and use such solution to predict the M (T ) curve in the EMFA approximation.In Ref. [26], the authors solve numerically all the coupled equations of motion for all the spins in the system.However, it is worth noting that the approach to including the environment presented here, without the mean-field approximation, could be adapted for ASD simulations.In such case, it would naturally lead to non-Markovian effects in the dynamics, which have been experimentally observed in ultra-fast magnetization dynamics [47] and have otherwise been introduced ad-hoc into ASD calculations.

V. SUMMARY AND CONCLUSIONS
In this article, we presented two different methods to account for quantum corrections in ASD.
Firstly, we showed that, introducing a power-law temperature rescaling, the classical MFA can accurately reproduce the temperature dependence of the quantum MFA magnetization.This power-law rescaling takes the form τ cl = (τ qu ) α(S0) where α(S 0 ) is an intrinsic property of the system and depends only on the spin length S 0 .The rescaling exponents α(S 0 ) obtained for Ni and Gd appear to be in very good agreement with the values of Ref. [25], used to reproduce experimental results via dynamical ASD simulations.Here, we give a prediction of the rescaling exponent α as a function of S 0 , which does not require a priori knowledge of the microscopic description of the system and is not a fit of experimental results.Our approach can be easily combined with ASD to account for quantum effects.This result gives a clear prescription on how to account for quantum effects in ASD simulations.To do this, it is sufficient to calculate the rescaling exponent α(S 0 ) by using Eq. ( 5) and run ASD simulations at a temperature τ sim = τ α(S0) exp .To avoid the need to "artificially" rescale while still accounting for quantum effects, we presented a second method (see Sec.(IV)).In this section, we extended the MFA to the study of a semi-classical open Heisenberg model in contact with environmental degrees of freedom, whose quantum statistics bring in quantum effects, e.g., that of phonons.The parameters of this quantum en-vironment+MFA approach can be directly calculated by means of ab initio calculations or extracted from experiments.We further introduced corrections to the (E)MFA by means of the SCGA.We found very good agreement between the magnetization curves of nickel predicted by this method and the experimental results, particularly at low temperatures.The predicted magnetization is furthermore capable of reproducing the absolute experimental T c to good accuracy.This suggests that the quantum corrected classical spin model considered here can be used to faithfully predict the dynamics of magnetic materials at all temperatures.

Figure 1 .
Figure 1.Quantum and classical MFA magnetization.We show the normalized spin magnetization sz for a spin system with S0 = ℏ/2 (see Eq. (2)) evaluated with the MFA as function of the relative temperature τ cl/qu = T /T cl/qu z as a function of the relative temperatures τ cl/qu = T /T cl/qu c .Here, T cl/qu c are the critical temperatures marking the transition from the trivial solution s cl/qu z = 0 to the non-trivial solutions s cl/qu z ̸ = 0.In the absence of an external magnetic field, H ext = 0, the critical temperatures are T cl/qu c

Figure 2 .
Figure 2. Temperature-rescaled equilibrium magnetization for nickel and gadolinium.We show the comparison between the results for Ni (left panel) and Gd (right panel) of the classical MFA s cl z with a temperature rescaling according to Eq. (3) (blue dotted lines) and the quantum MFA s qu z (red dashed lines).Also shown are the phenomenological Kuz'min curves (gray solid lines).Additionally, the classical MFA predictions are shown (black dotted lines) to highlight their discrepancy to the experiments.The insets zoom into the low temperature behavior of the same curves.

Figure 3 .
Figure 3. Rescaling power α as a function of spin length S0.Here we plot the dependence of the rescaling parameter α, see Eq. (3), obtained by minimizing Eq. (4), as a function of spin length S0 (black diamonds).An analytical function, Eq. (5), that very closely matches the obtained exponents is also plotted (green line).

Figure 4 .
Figure 4. Operational flowchart of the EMFA approach combined with the self-consistent gaussian approximation (SCGA).

Figure 5 .
Figure5.Magnetization predicted with the EMFA approach.EMFA (orange triangles) is compared to the EMFA plus corrections from the self-consistent Gaussian approximation (EMFA+SCGA) (red circles), and the Kuz'min curve (gray solid line).The classical MFA prediction in terms of the physical temperature T is also shown (blue squares) highlighting its qualitative discrepancy with the experiment.Differently to Fig.2, these magnetization curves are evaluated against the physical temperature in units of Kelvin.If the EMFA and EMFA+SCGA curves were shown over their respective relative temperature τ , they would be the same.We note that both the EMFA and EMFA+SCGA accurately reproduce the low-temperature behavior of the experimental magnetization.Furthermore, EMFA+SCGA gives a relatively accurate prediction of the critical temperature.