# Realistic finite temperature simulations of magnetic systems using quantum statistics

## Abstract

We have performed realistic atomistic simulations at finite temperatures using Monte Carlo and atomistic spin dynamics simulations incorporating quantum (Bose-Einstein) statistics. The description is much improved at low temperatures compared to classical (Boltzmann) statistics normally used in these kind of simulations, while at higher temperatures the classical statistics are recovered. This corrected low-temperature description is reflected in both magnetization and the magnetic specific heat, the latter allowing for improved modeling of the magnetic contribution to free energies. A central property in the method is the magnon density of states at finite temperatures and we have compared several different implementations for obtaining it. The method has no restrictions regarding chemical and magnetic order of the considered materials. This is demonstrated by applying the method to elemental ferromagnetic systems, including Fe and Ni, as well as Fe-Co random alloys and the ferrimagnetic system GdFe .

## I Introduction

Computer simulations of materials have developed to a very powerful alternative and complement to experimental methods. Electronic structure methods are now commonplace to predict novel ground state properties of existing and yet to be synthesized materials, sometimes in combination with screening of materials databases and applying artificial intelligence tools Ortiz et al. (2009); Jain et al. (2013); Pizzi et al. (2016). These tools are only useful if the underlying methods are accurate enough to properly describe the materials. A common method to simulate finite temperature magnetic properties is to employ a two-step process where in the first step all material specific parameters are calculated within zero-temperature electronic structure theory to construct a simplified statistical physics model which is then solved in the second step by employing atomistic simulations. This methodology has been very successful of qualitatively predicting Curie temperatures and magnon spectra of transition metals, alloys and diluted magnetic semiconductors Rosengaard and Johansson (1997); Pajda et al. (2001); Bergqvist et al. (2004); Şaşıoğlu et al. (2005); Ležaić et al. (2007); Kudrnovský et al. (2008); Sato et al. (2010); Bergman et al. (2010); Bergqvist et al. (2013); Etz et al. (2015); Eriksson et al. (2017). An overwhelming majority of the reported simulations of this kind have hitherto been based on classical Boltzmann statistics at finite temperature instead of the proper Bose-Einstein quantum statistics. In principle, quantum statistics are expected to enhance the critical temperature with the factor , where is the spin quantum number. For large enough and at high enough temperatures, classical statistics has been shown to work reasonable well but at the same time, the limitations are well known. For instance, when using classical statistics the magnetization at low temperatures decreases too rapidly and the simulated low temperature magnetic specific heat is non-zero. Up until recently, there has been a lack of practical theories taking into account quantum statistics in atomistic simulations. An early effort was done by Körmann et. al. Körmann et al. (2011) where reducing a full atomistic description to a minimal nearest-neighbour model allows for a quantum Monte Carlo treatment of the reduced model. More recently, Woo. et. al. Woo et al. (2015) demonstrated how it is possible to incorporate quantum Bose-Einstein statistics in existing atomistic simulation framework and obtained qualitatively good results for Fe in both low and high temperature region. In this work, we have extended and modified this latter methodology to more complex situations, such as random alloys and anti- and ferri- magnetic systems. We have also examined how the magnetic density of states, that is a central property of the model as discussed below, can be modelled and how the resulting quantum statistics is affected. Attempts to create a completely parameter-free self-consistent theory are discussed.

The paper is organized as follows: In Section II we introduce the methodology and give the details of the calculations, in Section III we present our results and compare the different methodologies used for obtaining the magnon density of states. A summary and outlook of the study is finally provided in Section IV.

## Ii Theory

### ii.1 Control of temperature in atomistic simulations

A statistical physics treatment using either Atomistic spin dynamics (ASD)Skubic et al. (2008); Eriksson et al. (2017) or Monte Carlo simulations (MC) uses the Hamiltonian in Eq. (1) as a starting point.

(1) |

where denote the exchange interactions between two atomic magnetic moments at site and . We are using the sign convention in which ferromagnetic interactions are characterized by positive values of and antiferromagnetic with negative sign.

Within ASD, the temporal evolution of the atomic moments is governed by the Landau-Lifshitz-Gilbert (LLG) equations (see c.f. Ref. Eriksson et al., 2017 for a detailed account of the method). The inclusion of finite temperature effects is obtained through Langevin dynamics, by adding a stochastic Gaussian shaped field (white noise) to the effective magnetic field , that acts on each atomic moment . The stochastic, or thermal, field is characterized by

(2) | |||||

where and denote lattice sites, and the carteisian components and is the amplitude of the stochastic field given by

(4) |

where is the Gilbert damping parameter. Using the Fokker-Plank equation and searching for stationary solution, the amplitude of can be related to the temperature through the scaling factor Eriksson et al. (2017).

If equilibrium properties are desired, Monte Carlo (MC) simulations based on the Hamiltonian, Eq. (1) are normally more efficient than ASD simulations. The quantum corrections as outlined above can easily be included in Monte Carlo as well, and for simplicity we restrict the discussion here to the Metropolis algorithm Metropolis et al. (1953). These temperature considerations are fully transferable to other importance sampling MC methods, including the heat-bath algorithm Miyatake et al. (1986). The crucial step for advancing in time is the transition probability between two states and in the Markov chain. In the Metropolis algorithm, the following form of is used

(5) |

Thus the temperature dependence in both ASD and MC simulations is governed by the scaling factor .

### ii.2 Quantum statistics

In classical statistics, i.e. using Boltzmann distribution, the scaling factor , is given by

(6) |

where the superscript denotes classical for clarity. When transitioning from classical to quantum statistics, it is then this temperature dependent that is to be modified. Using quantum statistics with the Bose-Einstein distribution, Woo et. al Woo et al. (2015) showed that the relation in Eq. (6) is modified to the more complex form

(7) |

where is the magnon density of states (mDOS) at temperature and the energy of a magnon at a wave-vector . Here, the superscript denotes quantum for clarity.

An alternative to the proper Bose-Einstein expression in Eq. (7), is to perform the temperature rescaling based on empirical data, such as the experimental magnetization curve, which was proposed by Evans et. al Evans et al. (2015) using a simplistic model including free parameters. The big advantage with the present method is, depending on how the mDOS in Eq. (7) is modelled, the number of free parameters can be reduced, or even eliminated.

### ii.3 Magnon density of states

The most important quantity entering the expression for the rescaled field arising from the quantum fluctuations controlling the temperature in Eq. (7) is the temperature dependent magnon density of states (mDOS) and in this section we discuss two different methods for obtaining it.

#### Adiabatic magnon spectra

The adiabatic magnon spectra (AMS)Halilov et al. (1998) is directly connected to the real-space exchange interactions as appear in Eq. (1) through Fourier transformation of the Hamiltonian. A detailed account on how the procedure is performed is found in Appendix A and the result is a set of eigenvalues depending on the wave-vector . By integrating the AMS over all wave-vectors, the resulting mDOS is obtained. In practice, in this work we convolute the mDOS with a Gaussian to ensure numerical stability, and then normalize it so that .

The adiabatic mDOS is temperature independent and valid only at the ground state, i.e T=0 K. In order to create a temperature dependent mDOS from the AMS, we apply the quasiharmonic approximation (QHA) as was outlined in Ref. [Woo et al., 2015]. If is the highest energy that give contribution to the mDOS, which is also put as upper limit in the integration in Eq. (7), then at finite temperature ( the cutoff energy becomes rescaled as

(8) |

where is the critical temperature and the critical exponent associated with the magnetization in the 3D Heisenberg model. As pointed out in Ref. [Woo et al., 2015], this particular model assumes that above the cutoff frequency is zero and thus the quantum statistics are sharply reverted back to classical statistics.

#### Dynamical structure factor

The magnon dispersion can also be obtained by simulating the dynamical structure factor Bergman et al. (2010); Bergqvist et al. (2013); Etz et al. (2015) and identifying the peak values of for each wave-vector . In contrast to the adiabatic treatment, temperature effects from the Gilbert damping processes are included so that a finite temperature description of the magnon dispersion is obtained. It is however worth noting that the present implementation assumes a fixed magnitude of the magnetic moments and does thus not include longitudinal fluctuations which give rise to Stoner excitations and an additional damping mechanism for magnons (Landau damping). The simulated finite temperature is obtained from integrating over all sampled wave-vectors. In order to minimize the numerical noise in the simulated magnon spectra and resulting mDOS, we perform a sampling and post-processing protocol as described in Appendix. B.

### ii.4 Details of calculations

All first-principles calculations in this study was performed using a multiple-scattering (Korringa-Kohn-Rostoker, KKR) implementation of the density functional theory (DFT) in the local spin density approximation (LDA) as implemented in the SPR-KKR software Ebert et al. (2011). The calculations were performed scalar relativistically using full potential and employing a basis set consisting of -orbitals. The coherent potential approximation (CPA) was employed for treating the disordered FeCo alloy. The magnetic exchange interactions that appear in Eq. (1) were obtained from the magnetic force theorem using the LKAG formalism Lichtenstein et al. (1987) in the low temperature magnetic reference state. The exchange interactions are treated as temperature independent in this study, however this approximation can be lifted using either noncollinear reference statesSzilva et al. (2013) or using an extended spin modelPan et al. (2017). The atomistic simulations, either the Monte Carlo or atomistic spin dynamics, were performed using the UppASD software Skubic et al. (2008); Eriksson et al. (2017).

## Iii Results

### iii.1 Quasi-harmonic approximation (QHA)

#### Simple cubic model with nearest neighbour interactions

As a first application of the method, we study a model system consisting of moments distributed on a bulk lattice in simple cubic geometry interacting with only nearest neighbours. It is a very well characterized system that has been served as a benchmark for many Monte Carlo simulations in the past. Although no analytical result exist for its Curie temperature, it has been calculated to very high precision Peczak et al. (1991) to , which we also obtain (within error bars less than 0.005%) using our software. The magnon density of states obtained from AMS is displayed in Fig. 1.

In combination with the QHA in order to obtain the finite temperature mDOS, we show in Fig. 2 the results of the temperature dependence of the magnetization and the specific heat and compare those with the classical results. At low temperatures, it is well known that a classical description breaks down and the magnetization is decreasing too rapidly (linearly) with temperature and does not follow the Bloch law that is expected from linear spin wave theory analysis. Perhaps even more severe is the fact that the specific heat does not go to zero but to a constant value of 1 since each degree of freedom contributes with 1/2 and there are two transversal components. It can be seen as an analogy with the Dulong-Petit law for phonons where the classical value reaches a constant value of 3R, where R is the gas constant but where a quantum mechanical description in the form av either Einstein or Debye model obtains a vanishing specific heat when . Similar, using the quantum statistics in the simulations, the magnon specific heat is approaching zero at low temperatures. In fact, the particular form of the mDOS for this system causes the fluctuations to be very small in a large temperature interval and only become significant in a narrow region close to .

#### FeCo alloys

An important class of materials for technological applications are transition metal alloys such as FeCo. Fe-Co alloys exhibit a large saturation magnetization and Curie temperature, larger than the values for Fe or Co separately and are commonly used in read head devices in magnetic storage applications. In addition, around the composition of FeCo, the alloy has a measured Gilbert damping which is unusually low for a metallic systemSchoen et al. (2016). A great advantage with the present method is the applicability to not only elemental magnetic systems but also magnetic alloys and compounds which we here demonstrate by simulating finite temperature properties of FeCo.

Fe-Co around 50-50 composition can either be synthesized in an ordered magnetic compound forming CsCl (B2) crystal structure or as a completely random alloy in the bcc crystal structure. We have run calculations on both these structures and the results were found very similar between the two. Hence, we decide to present only the results of the random alloy. Calculation of the mDOS from AMS of a random alloy is slightly more involved than for an elemental material due to the random arrangement of the atoms. We therefore calculated the mDOS for 10 000 different disorder configurations and afterwards took the average of all those configuration in order to get the final averaged mDOS.

The results from Monte Carlo simulation are displayed in Fig. 3. The shape of the magnetization is similar to that of elemental Fe, Section III.2 with the key difference however that calculated is higher, around 1395 K which is in good agreement with experimentsKawahara et al. (2003). The quantization effects arising from using quantum statistics and QHA are also very evident when inspecting the internal energy of the system. At low temperatures, classical statistics overestimate magnon excitations causing a rapid increase of the internal energy in contrast to the internal energy from using quantum statistics that are much less affected. In the critical region around , the quantization effects drastically reduces and vanish at .

#### GdFe

A phenomenological correction using renormalized temperatures to the magnetization shape for ferromagnetic materials were used in an earlier studyEvans et al. (2015) but it has the shortcoming that it will not work for other magnetic orderings like antiferromagnets and ferrimagnets. In comparison, a great advantage with the present methodology it does not have this restriction which we will here demonstrate for the ferrimagnetic system GdFe. A ferrimagnetic is characterized by two sublattices with antiferromagnetic coupling but not fully compensating each other such as at low temperatures there is a finite (small) magnetization in the sample. GdFe is the prototype system for so called all-optical switchingStanciu et al. (2007); Chimata et al. (2015) in which a laser pulse is used as a stimulus for magnetic switching where both the magnetic sublattices obtain a moment orientation opposite to the initial state after the pulse. Crystalline GdFe, as treated here, has a trigonal unit cell (spacegroup R-3mH,166) with 3 different Fe sites and 2 different Gd sites. We employed the experimental crystal structure and volume in our calculations using the LDA+U approximation for the Gd -states with a value of the Hubbard parameter set to 6.7 eV and Hunds exchange set to 0.7 eV. Experimentally, the critical and compensation temperatures of GdFe are 721 K and 618 K, respectively Kadziolka-Gawel and Bajorek (2014).

In the upper panel of Fig. 4, calculated Gd and Fe sublattice magnetization from Monte Carlo simulations are displayed. Each of the sublattice magnetization is behaving similar to the ferromagnetic systems in the sense that the magnetization is decreasing too rapidly with temperature when classical statistics are employed. This is corrected when using quantum statistics and applying the QHA to the mDOS which result in a much more stable magnetization with respect to temperature. Experimentally, the rare-earth (Gd) moment has a more rapid decrease of its sublattice magnetization than the Fe, causing a zero total moment at a temperature lower than . This is the compensation temperature of a ferrimagnet. Taking the difference between our calculated sublattice magnetization, we do not obtain any compensation point, suggesting that an even more elaborate approach is needed to capture this behaviour, in particular additional longitudinal spin fluctuations of the moment. Nevertheless, we obtain a of around 650 K, in reasonable close agreement with experiments and keeping in mind that the calculated value is slightly dependent on chosen values of and . The value of is even more pronounced when inspecting the peak of the specific heat as displayed in the lower panel of Fig. 4. Similar to ferromagnets, using quantum statistics recover the correct behaviour at low temperatures with vanishing specific heat and in turn the entropy (not shown).

### iii.2 Simulated mDOS

In order to study the difference between using QHA on the adiabatic magnon spectrum compared to using a simulated mDOS we here consider the finite temperature properties of bcc Fe. Fe in the body centered cubic (bcc) structure is the prototypical ferromagnet and often serves as a benchmark for testing and validating any new features in atomistic simulations. Experimentally, Fe has a Curie temperature of 1043 K and the magnetization curve as function of temperature is well known from experiments. The size of the magnetic moments in Fe is relatively robust with respect to rotations and temperature and therefore rather well described by an Heisenberg Hamiltonian. By calculating exchange parameters up to three lattice constants from the low temperature ferromagnetic configuration and subsequent Monte Carlo simulations, we obtain of around 930 K, in agreement with previous studies Rosengaard and Johansson (1997); Pajda et al. (2001); Ruban et al. (2004). Improvements to can be obtained by including electron entropy effects and/or include longitudinal spin fluctuations Ruban et al. (2007); Pan et al. (2017), but that was not considered in this study.

As emphasized by Eq. (7), the mDOS is a central quality to extracting the quantum statistical description of thermal excitations. How the temperature dependence of the mDOS is modelled is thus directly determining the resulting statistics. Results for the finite temperature mDOS from AMS using the QHA, i.e. the same methodology as was used for the systems in Section III.1, are displayed in the upper panel of Fig. 5. Qualitatively the spectral weight is shifting towards lower frequencies with increasing temperatures and vanishes at , where the distribution become continuous i.e. classical statistics apply. The sharp transition from quantum to classical statistics at comes directly from the QHA model through the rescaling in Eq. (8). This can be motivated by the fact that the spin stiffness vanishes along with the magnetization at . However, there are experimental and theoretical indications of magnons even above Qin et al. (2017); Tao et al. (2005) and thus the validity of the QHA cut-off temperature can be questioned. The scaling in Eq. (8) does thus infer that is a free parameter in the model. In some aspects that can be a warranted feature, for example when using the experimental for an empirical description of the magnetization behaviour. Since the that is used as input for the QHA rescaling can be determined from MC simulations, using calculated exchange interactions, it can be argued that the is not a free parameter in that case. Another step towards a parameter-free description is to use a simulated mDOS instead of a model one. To obtain the simulated finite temperature magnon densities of states, we perform atomistic spin dynamics simulation to sample the excitation spectra represented by the dynamical structure factor , as described in Sec. II.3.1 and Appendix B, which is then integrated to obtain . We will in the following refer to this combination for obtaining the simulated mDOS as the QHB method. The resulting simulated mDOS curves are shown in the lower panel of Fig. 5 . At low temperatures well below , the results from QHA and the simulated QHB results agree very well as expected. For elevated temperatures the trend of the QHB mDOS follows the trend from the QHA based data, with spectral weight transfer to lower frequencies with increasing temperatures. However, at high temperatures close to , the difference of mDOS from QHA and QHB becomes evident. Instead of a vanishing mDOS as enforced by QHA, the QHB simulations yield an mDOS that fits rather well to a classical Boltzmann distribution.

The simulated mDOS is then used as input for the determination of the temperature rescaling factor for which MC simulations give the resulting temperature dependent magnetization as displayed in Fig. 6. For comparison with experiments, an empirical magnetization curve according to Kuz’minKuz’min (2005) is also shown in Fig. 6. The Kuz’min curve follows an analytic expression for the shape of the magnetization that reproduces the experimental curve by matching the low temperature Bloch law with the high temperature critical behaviour resulting in the expression

(9) |

where , , and parameters with values 0.35 and 4.0, respectively.

Using classical statistics, we obtain the usual shortcomings that the magnetization drops too rapidly at low temperatures and has a linear temperature dependence instead of the expected Bloch behaviour. However, using quantum statistics and QHA, apart from small deviations close to which arises from finite size effects from the simulations, we obtain a similar shape as the empirical curve from Eq. (9) (using our calculated and M(0)). This agreement is obtained regardless if the underlying mDOS is taken from AMS or the simplified quartic model mDOS from Ref.[Woo et al., 2015] (not shown). From the observation that the magnetization curve is very similar for the simplified quartic model and for the AMS mDOS it can be noted that the quartic model is an acceptable approximation under QHA since minor differences in the high energy region of the mDOS does not seem crucial for determining the M vs T behaviour. Why the shape of the high energy mDOS is less important can readily be understood by the exponential prefactor to the mDOS in the determination of the scaling factor as expressed in Eq. (7).

The simulated temperature dependent mDOS within the QHB method, is much more general in the sense that it does not require any additional parameter and should thus yield a more realistic description of the magnon properties at finite temperatures. At lower temperatures, the magnetization is similar to the QHA treatment but at elevated temperatures two pronounced differences can be noticed. The curvature of the QHB curve is less pronounced, resulting in a flatter shape compared to both the Kuz’min and the QHA curves. More important, the obtained from the QHB curve differs as well. That means that when is not given as a parameter to the model, as in the QHA case, the use of quantum instead of classical statistics also affects the critical temperature of the simulated system.

In the results above, the QHB mDOS was obtained using classical statistics and then applied to get a quantum statistical description of the following MC simulations. Since the use of quantum statistics affects the thermodynamics of the simulations, good arguments can be made that the QHB mDOS should also be obtained from quantum statistical simulations. This can be achieved in a self-consistent way by performing first a classical simulation, to get the first mDOS, and then successive QHB simulations, always using the mDOS obtained from the previous run, until the M vs T curves are converged. How this iterative approach performs can be seen in Fig. 7 where it is found that a self-consistent result is obtained after 3-4 QHB simulations. Comparing the self-consistent results in Fig. 7 with the single-shot output in Fig. 6 it can be seen that is slightly increased further and also that the shape of the curve becomes more similar to the Kuz’min and QHA results.

In the case of bcc Fe considered here, the obtained increase in when using the self-consistent QHB mDOS could be seen as unfortunate since the classical simulations give a better agreement with the experimental compared to the quantum statistical QHB simulations. However, the increase in is consistent with picture that the critical temperature would have an enhancement factor of , where is the spin quantum number, in a fully quantum mechanical treatmentLichtenstein et al. (1987) (the enhancement factor is equal to 1 in classical statistics). Following the similar arguments as in Ref. [Lichtenstein et al., 1987], the good agreement in classical statistics may be due to mutual error cancellation of various effects, for instance the use of temperature independent exchange interactions. By correcting the statistics, perhaps it is also needed to combine with more sophisticated methods to compensate, such as temperature dependent exchange interactions.

As an example where an enhanced from using quantum statistics is improving the theoretical results with respect to comparison with experiments, we show results for QHA and QHB simulations for fcc Ni in Fig. 8. Here we see that the increase in from 330 K in the classical case to 661 K, as obtained by self-consistent QHB simulations, is indeed an improvement when compared to the experimental value of Kuz’min (2005). Since Ni has a smaller moment and thus lower spin quantum number than Fe, quantum effects are indeed expected to be larger in Ni. However, it is also seen that the steep shape of the simulated M vs T curve is strongly exaggerated compared to the experimental curve. We have not been able to figure out why the simulated M vs T curve takes this form but there are several potential reasons for it. First of all, it is well known that the exchange interactions in Ni are strongly temperature dependent and differs considerably depending on reference state. Related, the finite temperature description of Ni from a pure Heisenberg model is very poor, in particular it is needed to also include longitudinal spin fluctuations that are responsible for stabilizing a finite magnetic moment at high temperatures. Nevertheless, in many cases it is mostly interesting to obtain reliable values of and in those cases, the present method may provide that without too much complications.

### iii.3 Thermodynamic properties

We now turn our attention towards how a quantum statistical description can improve the calculation of the free energy which is crucial in order to correctly determine phase stabilities. Returning to bcc Fe, relevant thermodynamic properties are gathered in Fig. 9. The specific heat , that measures the amount of fluctuations in the system, has the familiar shape with a peak at and large difference between classical and quantum treatment where the latter is required in order to obtain the physical result when K. The magnetic entropy can be obtained from Monte Carlo simulations in two alternative ways, either by using a model or by direct thermodynamic integration. In the paramagnetic limit, the (transversal) magnetic entropy , where is the size of each atomic moment. Without accounting for longitudinal spin fluctuations, using the value , the high temperature magnetic entropy has value of 1.16, indicated by a dashed line in the middle panel of Fig. 9. In the opposite end, from the third law of thermodynamics, it is required that the entropy is zero when approaching zero temperature. In order to have a simplified model between these two limits, Heine and JoyntHeine and Joint (1988) formulated a model interpolating between these two limits, later employed by GrimvallGrimvall (1989) for fcc Fe. In this model, which we denote the Heine-Joynt model, the magnetic entropy has the following expression:

(10) |

where and is the average angle between nearest neighbour magnetic moments. The angle and its temperature dependence is easily obtained in Monte Carlo simulations from the static correlation function and in turn the magnetic entropy using Eq. (10). An alternative and more general method to obtain magnetic entropy is through thermodynamic integration

(11) |

From this expression, it is clear that classical statistics fail at low temperatures since at any non-zero temperature, the entropy will be finite. From our results, we obtain a lower entropy in the whole temperature interval from the thermodynamic integration compare to the Heine-Joynt model. It is however worth stressing that only the transversal fluctuations contribute to the entropy in these calculations, at elevated temperatures the longitudinal fluctuations of the momentsPan et al. (2017) have important contributions to the entropy. Once the magnetic entropy is obtained, it is trivial to obtain the (Helmholtz) free energy , where is the internal energy which is shown in the bottom panel of Fig. 9.

## Iv Summary

A general improvement of finite temperature properties of magnetic systems from atomistic simulations are found using quantum statistics instead of classical statistics. The quantum effects enters through the temperature dependent magnon density of states. A simple but very effective method is the quasi harmonic approximation where the easy obtainable zero temperature mDOS is rescaled at finite temperatures and vanishing above the critical temperature in which the method is reverting back to classical description. The methodology has here been shown to work not only for ferromagnetic systems but also for anti- and ferri- magnetic magnetic systems as well as for disordered alloys. Perhaps the main shortcoming is the required knowledge of the critical temperature which separate the quantum and classical description such that the simulations must be performed in two steps, first to determine the critical temperature using classical statistics and secondly rerun the simulations using quantum statistics. To eliminate this shortcoming, we have constructed a more sophisticated approach where the temperature dependent magnon density of states is self-consistently obtained through magnetodynamical simulations. With this scheme it, the low-temperature behaviour is closely agreeing with the QHA method, where both methods have the important properties that the M vs T curve follows the Bloch and that the magnetic specific heat goes to zero at low temperatures. The latter property is crucial in order to be able to extract a good description of the magnetic contribution of the free energy. At larger temperatures, we notice that our suggested scheme gives an enhanched critical temperature compared to when using classical statistics. The method is already suitable to employ for atomistic simulations and free energy calculations and have the prospect of being improved further if longitudinal fluctuations and/or coupled spin-lattice contributions to the thermodynamics of the simulated systems are taken into account. Such efforts are ongoing.

###### Acknowledgements.

The authors thank Danny Thonig for fruitful discussions. We acknowledge support from the VR (Swedish Research Council), SeRC (Swedish e-Science Research Centre), GGS (Göran Gustafssons Stiftelser), eSSENCE, and the CEA-Enhanced Eurotalents program, co-funded by FP7 Marie Skłodowska-Curie COFUND Programme (Grant Agreement n^{o}600382). The computations were performed on resources provided by SNIC (Swedish National Infrastructure for computing) at NSC (National Supercomputer Centre) in Linköping, Sweden.

## Appendix A Adiabatic magnon spectra and mDOS

Let denote the Fourier transform of the exchange interaction between atom type and with a wave-vector lying in the Brillouin zone (BZ). is calculated as

(12) |

where is a position vector connecting the two sites. If there are atoms per unit cell, for each wave-vector , the magnon energies will then be given by the eigenvalues of a general matrix Halilov et al. (1998); Anderson (1963) here expressed in block form

(13) |

From the magnon frequencies, the adiabatic mDOS is obtained by summing up the energies as

(14) |

## Appendix B Simulated magnon spectra and mDOS

By measuring the trajectories of the simulated atomic moments, the time and space correlation function is obtained

(15) |

The correlation function defined in Eqn. (15) can thus describe how the magnetic order evolves both in space and over time. The perhaps most valuable application of is however obtained by a Fourier transform over space and time to give the dynamic structure factor

(16) |

Due to a combination of finite size effects, finite sampling time, temperature fluctuations and damping in the simulations, the raw output of the magnon density of states from simulations has associated noise that needs to be taken care of. For this purpose we employ a Hann window function for the time-domain transform and also apply a Gaussian convolution of the summed up mDOS combined with a normalization. In addition to the noise, there are also other numerical difficulties with the sampling of . While the correlation function is formulated in cartesian coordinates , the proper correlations stemming from the magnons in the system is only sampled in the spatial directions aligned perpendicularly to the magnetization. That means that in order to get a good mDOS from the dynamic structure factor, the perpendicular components of it must be identified, and more importantly, the longitudinal component should not be sampled. For that reason, we perform a rotation of the measured trajectories towards the quantization axis so that and is sampled. At finite temperatures, there can however be a small drift of the magnetic system due to the temperature dependent stochastic fluctuations present in the simulations so defining a static quantization axis might prove difficult. As a result, the simulated mDOS might show spurious peaks for the low-energy part of the spectrum, in particular at . While these peaks can be partially suppressed by Gaussian convolution, we have found that a better way to remove the spurious peaks is to update the quantization axis over time during the simulations. Physically, this can be motivated by the fact that the Goldstone mode of a ferromagnet should not carry any spectral weight in the mDOS. The quantization axis update is however performed slowly compared to the time-scale of the simulations in order to minimize the effect on the simulated spectra at .

### References

- C. Ortiz, O. Eriksson, and M. Klintenberg, Computational Materials Science 44, 1042 (2009).
- A. Jain, S. P. Ong, G. Hautier, W. Chen, W. D. Richards, S. Dacek, S. Cholia, D. Gunter, D. Skinner, G. Ceder, and K. A. Persson, APL Materials 1, 011002 (2013).
- G. Pizzi, A. Cepellotti, R. Sabatini, N. Marzari, and B. Kozinsky, Computational Materials Science 111, 218 (2016).
- N. M. Rosengaard and B. Johansson, Phys. Rev. B 55, 14975 (1997).
- M. Pajda, J. Kudrnovský, I. Turek, V. Drchal, and P. Bruno, Phys. Rev. B 64, 174402 (2001).
- L. Bergqvist, O. Eriksson, J. Kudrnovský, V. Drchal, P. Korzhavyi, and I. Turek, Phys. Rev. Lett. 93, 137202 (2004).
- E. Şaşıoğlu, L. M. Sandratskii, P. Bruno, and I. Galanakis, Phys. Rev. B 72, 184415 (2005).
- M. Ležaić, P. Mavropoulos, and S. Blügel, Applied Physics Letters 90, 082504 (2007), http://dx.doi.org/10.1063/1.2710181 .
- J. Kudrnovský, V. Drchal, and P. Bruno, Phys. Rev. B 77, 224422 (2008).
- K. Sato, L. Bergqvist, J. Kudrnovský, P. H. Dederichs, O. Eriksson, I. Turek, B. Sanyal, G. Bouzerar, H. Katayama-Yoshida, V. A. Dinh, T. Fukushima, H. Kizaki, and R. Zeller, Rev. Mod. Phys. 82, 1633 (2010).
- A. Bergman, A. Taroni, L. Bergqvist, J. Hellsvik, B. Hjörvarsson, and O. Eriksson, Phys. Rev. B 81, 144416 (2010).
- L. Bergqvist, A. Taroni, A. Bergman, C. Etz, and O. Eriksson, Phys. Rev. B 87, 144401 (2013).
- C. Etz, L. Bergqvist, A. Bergman, A. Taroni, and O. Eriksson, J. Phys: Condens. Matter 27, 243202 (2015).
- O. Eriksson, A. Bergman, L. Bergqvist, and J. Hellsvik, Atomistic Spin Dynamics: Foundations and Applications (Oxford University Press, 2017).
- F. Körmann, A. Dick, T. Hickel, and J. Neugebauer, Phys. Rev. B 83, 165114 (2011).
- C. H. Woo, H. Wen, A. A. Semenov, S. L. Dudarev, and P.-W. Ma, Phys. Rev. B 91, 104306 (2015).
- B. Skubic, J. Hellsvik, L. Nordström, and O. Eriksson, Journal of Physics: Condensed Matter 20, 315203 (2008).
- N. Metropolis, A. W. Rosenbluth, M. N. Rosenbluth, A. H. Teller, and E. Teller, The Journal of Chemical Physics 21, 1087 (1953).
- Y. Miyatake, M. Yamamoto, J. J. Kim, M. Toyonaga, and O. Nagai, Journal of Physics C: Solid State Physics 19, 2539 (1986).
- R. F. L. Evans, U. Atxitia, and R. W. Chantrell, Phys. Rev. B 91, 144425 (2015).
- S. Halilov, H. Eschrig, A. Perlov, and P. Oppeneer, Phys. Rev. B 58, 293 (1998).
- H. Ebert, D. Ködderitzsch, and J. Minâr, Reports on Progress in Physics 74, 096501 (2011).
- A. Lichtenstein, M. Katsnelson, V. Antropov, and V. Gubanov, Journal of Magnetism and Magnetic Materials 67, 65 (1987).
- A. Szilva, M. Costa, A. Bergman, L. Szunyogh, L. Nordström, and O. Eriksson, Phys. Rev. Lett. 111, 127204 (2013).
- F. Pan, J. Chico, A. Delin, A. Bergman, and L. Bergqvist, Phys. Rev. B 95, 184432 (2017).
- P. Peczak, A. M. Ferrenberg, and D. P. Landau, Phys. Rev. B 43, 6087 (1991).
- M. Schoen, D. Thonig, M. Schneider, T. Silva, H. Nembach, O. Eriksson, O. Karis, and J. Shaw, Nat. Phys. 12, 839 (2016).
- K. Kawahara, D. Iemura, S. Tsurekawa, and T. Watanabe, Mater. Trans. 44, 2570 (2003).
- C. D. Stanciu, F. Hansteen, A. V. Kimel, A. Kirilyuk, A. Tsukamoto, A. Itoh, and T. Rasing, Phys. Rev. Lett. 99, 047601 (2007).
- R. Chimata, L. Isaeva, K. Kádas, A. Bergman, B. Sanyal, J. H. Mentink, M. I. Katsnelson, T. Rasing, A. Kirilyuk, A. Kimel, O. Eriksson, and M. Pereiro, Phys. Rev. B 92, 094411 (2015).
- M. Kadziolka-Gawel and A. Bajorek, Hyperfine Interact 226, 321 (2014).
- A. V. Ruban, S. Shallcross, S. I. Simak, and H. L. Skriver, Phys. Rev. B 70, 125115 (2004).
- A. V. Ruban, S. Khmelevskyi, P. Mohn, and B. Johansson, Phys. Rev. B 75, 054402 (2007).
- H. J. Qin, K. Zakeri, A. Ernst, and J. Kirschner, Phys. Rev. Lett. 118, 127203 (2017).
- X. Tao, D. P. Landau, T. C. Schulthess, and G. M. Stocks, Phys. Rev. Lett. 95, 087207 (2005).
- M. D. Kuz’min, Phys. Rev. Lett. 94, 107204 (2005).
- V. Heine and R. Joint, Europhys. Lett. 5, 81 (1988).
- G. Grimvall, Phys. Rev. B 39, 12300 (1989).
- P. W. Anderson, Solid State Physics 14, 99 (1963).