Full-spin-wave-scaled stochastic micromagnetism for mesh-independent simulations of ferromagnetic resonance and reversal

In this paper, we address the problem that standard stochastic Landau-Lifshitz-Gilbert (sLLG) simulations typically produce results that show unphysical mesh-size dependence. The root cause of this problem is that the effects of spin-wave fluctuations are ignored in sLLG. We propose to represent the effect of these fluctuations by a full-spin-wave-scaled stochastic LLG, or FUSSS LLG method. In FUSSS LLG, the intrinsic parameters of the sLLG simulations are first scaled by scaling factors that integrate out the spin-wave fluctuations up to the mesh size, and the sLLG simulation is then performed with these scaled parameters. We developed FUSSS LLG by studying the Ferromagnetic Resonance (FMR) in Nd2Fe14B cubes. The nominal scaling greatly reduced the mesh size dependence relative to sLLG. We then performed three tests and validations of our FUSSS LLG with this modified scaling. (1) We studied the same FMR but with magnetostatic fields included. (2) We simulated the total magnetization of the Nd2Fe14B cube. (3) We studied the effective, temperature- and sweeping rate-dependent coercive field of the cubes. In all three cases, we found that FUSSS LLG delivered essentially mesh-size-independent results, which tracked the theoretical expectations better than unscaled sLLG. Motivated by these successful validations, we propose that FUSSS LLG provides marked, qualitative progress towards accurate, high precision modeling of micromagnetics in hard, permanent magnets.


INTRODUCTION
Finite element micromagnetic modeling has been proven to be a reliable tool to describe many magnetic phenomena at finite temperatures. Usually, the micromagnetic model utilizes a form of the Landau-Lifshitz-Gilbert (LLG) equation. If the computation requires varying the temperature, as for example in simulations of heat-assisted magnetic recording or permanent magnets in electric motors and generators, thermal excitations and fluctuations need to be represented in the LLG method.
The most popular approach to deal with thermal excitations in micromagnetics is to use a Landau-Lifshitz-Bloch (LLB) based equation 1,2 which combines the LLG equations for low temperatures and the Bloch equations for high temperatures. In contrast to the LLG equations, in LLB the magnetization magnitude is no longer conserved, moreover, the transversal and longitudinal components have different damping parameters. The LLB gives good results for temperatures close to and higher than the Curie temperature, but it achieves this success by introducing several additional temperature-dependent parameters. At high temperatures, the magnetization reverses linearly by changing its length and orientation, a process which can be perfectly described using the LLB equation 3 . Constructing these parameters involves a considerable amount of effort, though, as it requires a multi-scale simulation approach including ab initio methods and atomistic simulations, typically necessitating additional assumptions and phenomenologies 4,5 .
In this paper, we focus on the stochastic Landau-Lifshitz-Gilbert equation (sLLG). For temperatures far smaller than the Curie temperature, as in the case of permanent magnet applications, the magnetization will fluctuate at the surface of its unit sphere. Switching will not occur by linear reversal. Hence, the sLLG is a good choice. However, unless we do a scaling of the parameters, the result will strongly depend on the mesh size. The here presented FUSSS version of the LLG fixes this problem. A key advantage of sLLG for finite-temperature micromagnetism is that it requires less phenomenological considerations 6,7 . In the sLLG approach, the thermal perturbations are represented by adding white noise to the LLG equation, turning it into a Langevin-type stochastic differential equation 8 . The white noise is added to the effective field of the equation in the form of a stochastic field.
However, simulating magnetization dynamics with sLLG still requires either an accurate computation of atomistic-level spin models 9,10 , or an approximate finite element (FE) calculation with mesh size of atomistic length scales. The intrinsic properties serving as input parameters for such finite element calculations are fixed on the single spin level at 0 K and are usually taken from ab initio calculations. Sadly, due to the high demand of such calculations on computing resources, simulations at the atomistic scale are limited to a sample size of few nanometers.
To calculate magnetic behavior of samples on the micromagnetic scale, usually ranging from nanometers up to a few micrometers, the mesh size has to be increased. The term mesh size is used here to describe the average edge length of the elements in the mesh. Since micromagnetism is a continuum theory, its parameters represent thermal averages over the finite elements that are much larger than the material's unit cell. Simulations that use a coarser FE mesh but retain the atomistic parameters fail to include fluctuations on length scales between the interatomic spacing and the mesh size. A class of such fluctuations is the short wavelength spin waves that reduce the magnetization of the mesh-elements, a key parameter of the sLLG simulation. Therefore, the saturation magnetization M s computed with the sLLG equation using the atomistic magnetization on a coarse mesh is too high. Such discrepancies increase with increasing temperatures. Many works address this problem by adopting phenomenological, effective parameters to match the sLLG simulations on a finite difference or finite element mesh with experiments. Without such adaptation, the sLLG results are reliable only at low temperatures.
Another undesirable consequence of the naive sLLG methods that do not account for the fluctuations on intra-mesh length scales is that the results of the simulations are dependent on the mesh size: larger meshes ignore more fluctuations. This is clearly an unphysical state of affairs.
An inspired, quantitative method to account for intra-mesh fluctuations is to construct and to use renormalization group equations to integrate out intra-mesh fluctuations and represent them with scale-dependent, renormalized mesh parameters in sLLG FE simulations [11][12][13][14] . Grinstein and Koch 11 demonstrated the power of this method by modeling the temperature dependence of the magnetization of iron with a Heisenberg model. They used the Heisenberg renormalization group equations to generate scale-dependent, renormalized mesh parameters, which then they used in their sLLG simulations of the magnetization of large samples. Subsequently, Kirschner et al. 15 performed Monte Carlo simulations on an atomistic level and calculated the average spontaneous magnetization for much coarser cells at different temperatures. These cell/mesh size dependent macro-spins were then used in sLLG simulations on length scales two orders of magnitude larger than the unit cell. Finally, some papers use simplified versions of the renormalization group and call the procedure coarse graining, such as Behbahani et al. 14 for their sLLG calculation of hysteresis loops for magnetite nanorods.
However, these renormalization group approaches still have limitations. For example, renormalization group works only in the leading logarithmic approximation, which holds only if the spin waves have a gapless, purely quadratic dispersion. The need to use this idealized dispersion forced the incorporation of the anisotropy K u and the magnetic field h only as perturbative corrections. Further, the effect of the actual crystal symmetries on the spinwave dispersion were also ignored when the momentum integrals were performed in idealized, isotropic spheres. When used for asymptotic calculations, such as near the critical temperature and at very long length scales, these approximations are justified. However, in most sLLG simulations the system is far from criticality, and the length scales of the FE mesh are intermediate, only a few nanometers. For both of these reasons, the quadratic gapless spinwave dispersion approximation is quite poor. It is especially poor for hard, permanent magnets because of their large anisotropy K u . It gets worse still when the magnetic fields are high. And finally, reducing the analysis to the leading logarithmic approximation is justified if the renormalization group equations are integrated out to a very long correlation length, which is the case only very close to criticality. In contrast, the justification and therefore the accuracy of this leading logarithmic approximation is limited when only a limited range of wavelengths/length scales need to be integrated out, such as in the present case, when the integration proceeds from the atomistic scale only to the mesh size. For all the above reasons, for hard magnets in high fields with typical mesh sizes, the accuracy of sLLG simulations with renormalization group-corrected parameters is limited. This is consistent with the fact that the original verification of these renormalization group ideas was performed only for the soft magnet iron, in zero field, and at temperatures approaching the critical temperature T c 11 , where the approximations are the most defensible, as properly discussed in that paper.
To address the above limitations and challenges, in the present paper we report the development of an improved renormalization group, or scaling, method. (1) Instead of the exchange-only, gapless, quadratic spin-wave dispersion, forced by the leading logarithmic approximation, we retain the full-spin-wave dispersion with the anisotropy K u and finite magnetic field h in it. This way we avoid treating these latter terms only perturbatively. Specifically, we use the spin-wave dispersion of the hard magnetic compound Nd 2 Fe 14 B at 300 K temperature, determined and confirmed by experiments 16 . (2) Using these experimental dispersions also retains the proper, real spin-wave dispersion that represents the crytalline symmetry of the magnetic material. With these two modifications, we integrate out the intra-mesh spinwave fluctuations from the atomistic scale to the mesh size. Please note that the proposed method works with any suitable spin-wave dispersion, either measured experimentally or derived theoretically. This integral gives rise to a scaling factor for the magnetization parameter M s appropriate for the sLLG mesh size. We do not use the classical relation to scale the anisotropy constant K u . Instead, we perform simulations of ferromagnetic resonance in order to derive a scaling law that is self-consistent with the simulations. This is done by fitting K u to shift all FMR curves to the same, shared bias field. For the exchange constant Ax we saw very little influence on the FMR results, hence we stuck with the classical relation of A x / M 2 s . With such scaled parameters, micromagnetic simulations of magnetization reversal are nearly mesh size independent. However, we compare the results of our improved scaling method with the results obtained by the classical scaling relation.
We will demonstrate the efficiency of our method in the context of analyzing ferromagnetic resonance (FMR). FMR measurements are a widely used technique to investigate dynamic magnetic behavior, especially to determine damping effects in magnetization dynamics 17 . Micromagnetic FMR simulations have been performed for granular perpendicular media for magnetic recording to extract the Gilbert damping constant 18 . In this paper, we use our spin-wave renormalized sLLG method to simulate the magnetic field dependence of the FMR curves and show that the simulations are in good agreement with experiments, and reassuringly, are mesh size independent. Although our method is applicable for other spatial discretization methods like finite differences, in this paper we focus on finite elements to proof feasibility.

Development of FUSSS LLG by simulating ferromagnetic resonance
Often ferromagnetic resonance (FMR) experiments or simulations are used to determine the effective damping constant in materials of interest 12,18 . In this work, we carry out FMR simulations to critically test the developed full-spin-wave-scaled sLLG theory by checking whether it indeed delivers results independent of the finite element mesh size. In the simulations, we have adopted the often-used Gilbert damping value of α = 0.01.
In the FMR-simulation, the previously prepared cube is exposed to an oscillating field with a maximal amplitude of μ 0 H AC = 5 mT applied orthogonally to the cube's anisotropic easy axis in xdirection (see Fig. 1).
The frequency is chosen to be f AC = 216 GHz. Neglecting the demagnetizing field, the magnetic moments are in resonance at a bias field of μ 0 H bias = 1 T, resulting in an FMR peak at this field: and γ ¼ 1:7608596 10 11 T À1 s À1 : γ is the gyromagnetic ratio. H ani is the theoretical anisotropy field calculated from the macroscopic intrinsic properties denoted by the superscript L.
In order to distinguish between the different length scales and their respective magnetic properties, we use superscripts a for atomistic length scale, l for minimum scale of micromagnetic simulations, and L for macroscopic length scale as in measurements.
The bias field is applied parallel to the cube's easy axis. Initially, the cube's magnetization is saturated parallel to H bias . At a temperature of 300 K the time evolution of the magnetization configuration is calculated by solving the stochastic LLG. The simulations are repeated with different mesh sizes and with different values of the bias field from 0.6 T to 1.9 T. Each simulation was performed for 1 ns. The transients settled after 0.4 ns-the data recorded after this transient time were accepted as the results. Sequentially the x-component of the magnetization, i.e., the component in direction of H AC , was extracted from the remaining signal. The FMR curves were calculated by taking the maximal magnitude,m x , in the frequency spectrum of the x-component of m.
The results of various sLLG simulations for different mesh sizes are shown in Fig. 2. As shown in Fig. 2a, when the sLLG did not use parameters scaled by the spin-wave fluctuations, the simulations yielded strongly mesh size-dependent results. In contrast, Fig. 2b shows that when the full-spin-wave-scaled sLLG used parameters scaled by the spin-wave fluctuations, then the simulations yielded results essentially independent of the mesh size, as discussed below.
The FMR simulations were performed for a uniformly meshed hard magnetic Nd 2 Fe 14 B cube, wherein the mesh size l was varied in the 1-10 nm range. Please note that for FMR experiments as depicted in Fig. 1 we expect an almost uniform magnetization and only slight deviations of the magnetization from the easy axis. Therefore the mesh size can be larger than the domain wall width.
As Fig. 2a shows, without scaling of the parameters, the magnetic (bias) field defining the center of the FMR peak shifted substantially from the external bias field of H bias = 1 T to the much higher bias fields of H bias = 1.6 T, as l varied from 10 nm to 1 nm. This large, 60% shift demonstrates one more time the unphysical dependence of measurable quantities on the mesh size. For completeness, we note that the Gilbert-damping-related FMR line-width (the FWHM of the peaks) stayed essentially invariant at 255mT (σ 2 = 0.095) as the mesh size l was varied.
Feng and Visscher 12 introduced a mesh size dependent damping constant for a system including only exchange and thermal energy. The effective rescaled damping constant was found to increase with temperature and computational cell size. In this work, we focus on permanent magnets, where in addition to the exchange the magnetocrystalline anisotropy is the dominating energy term. For the cell size of about 4 nm and 300 K the ratios of thermal energy to anisotropy energy and to the exchange energy are 0.015 and 0.134, respectively. Following the results of Feng and Visscher, our computational experiments are in the lowtemperature regime where no or only minor scaling of the damping parameter is required.
When the parameters were scaled with the mesh size according to the full-spin-wave-scaling derived in the Methods section of this paper with scaling factors s M (l), s A (l) = s M (l) 2 , and s K (l) = s M (l) 3 , the FMR peaks shifted back to the vicinity of the bias field of H bias = 1 T. The creep of H bias with mesh size was reduced from above 60% to below 20%, a great improvement. The fact that there was some residual dependence still left is most likely due to the fact that our scaling factors are only the leading terms of the overall spin-wave fluctuation reduction of M l s . To eliminate even this residual mesh size dependence, we modified the scaling of K u , so that the FMR peak remains at the mesh size independent H bias = 1 T. To formulate the most natural scaling function for K u , we retained the power-law form and treated its exponent as the adjustable parameter. We found that the choice of the anisotropy scaling exponent K u / M 2:72 s kept the bias field of the FMR peak location fully mesh size independent, as shown in Fig. 2b. This is only a 9% adjustment of the scaling exponent from its Callen value of 3, a very respectable achievement from our first order approximate scaling factors.
The scaling functions we used for the full-spin-wave-scaled stochastic LLG are shown in Fig. 3. Additionally, these same scaling factors and the corresponding intrinsic property values for the used mesh sizes are listed in Table 1.
From now on, we will refer to the above developed and described full-spin-wave-scaled stochastic LLG as FUSSS LLG. Fig. 1 FMR spectra calculation setup. Scheme of setup to calculate the FMR spectra for a permanent magnetic cube with various mesh sizes. The bias field is applied in direction of the z-axis, i.e., parallel to the magnetocrystalline anisotropy axis, and the oscillating field acts in x-direction. The x-component of the magnetization is used to calculate the response of the system for each bias field.

Tests and validation of FUSSS LLG
We performed several tests and validation of our FUSSS LLG. First, we repeated the FMR simulations by including the magnetostatic field. The FMR curves show two peaks (see Fig. 4).
This curve was obtained using the Fourier transform of the average magnetization of the entire cube. To analyze the spectrum we also computed the resonance of the local magnetization at different probe points. The analysis of local magnetization dynamics shows that the peak at the lower H bias corresponds to the local response at the center of the cube. In this bulk mode, the major part of the cube is in resonance. The second peak at higher H bias values is caused by resonance at points near the top and bottom face (near the center of the face or in the middle of an edge). The edge mode arises from a different local demagnetizing field. In addition, there might be magnetostatic interaction between the modes. However, in the following we show the scaling parameters can be derived neglecting magnetostatics.
The existence of these two modes was confirmed experimentally in FMR investigations of thin films 19 . The two peaks were centered at H bias ≈ 0.88 T and H bias ≈ 1.31 T, for the bulk mode and for the edge mode, respectively. As before, when the intrinsic parameters were not scaled, the sLLG curves shifted to higher bias fields as seen in Fig. 4a. For l = 1 nm the lower field peak shifted from 0.88 T by 0.52 T, again a 60% shift. To demonstrate the predictive power of our FUSSS LLG, we then repeated the simulations with the intrinsic parameters scaled by the scaling factors determined in the no-magnetostatic-field simulations earlier (see Table 1). Figure 4b shows that the FUSSS LLG also found the two resonance peaks, and their center bias fields were almost independent of the finite element mesh size l.
Second, we simulated the equilibrium value of the magnetization of a Nd 2 Fe 14 B cube with an edge length of 40 nm at 300 K temperature. For each chosen mesh size the magnetic state was relaxed for 0.5 ns after initial saturation in the easy-axis direction. The magnetostatic field was not taken into account in these simulations. Due to the stochastic field in the sLLG, the magnetic moments fluctuated around the easy axis. In Fig. 5 the mean value over 2 ns of the magnetization component in easy direction    μ 0 〈M z 〉 is plotted against the mesh size. The macroscopic magnetization was measured to be μ 0 M L z ¼ 1:61 T. Without scaling, the sLLG produces reduced values with decreasing mesh size. For l = 1 nm the mean magnetization is reduced by 6% to 1.51 T, as for smaller meshes the restoring force is smaller. Then we repeated the same calculation with FUSSS LLG, using scaled intrinsic parameters. As shown in the same figure, the magnetization reduction has been properly compensated, and the unphysical mesh-size dependence essentially eliminated. The deviations have been reduced from 6% to <1%.
Third, we simulated the reversal of the magnetization by an applied field. Here we took the magnetostatic field into account. A cubic sample was magnetically saturated in the easy-axis direction and then reversed by changing the direction of the external field from the saturation direction to the opposite with a field sweep rate of v = 250 mT ns For every mesh size, the coercive field H c was extracted by averaging over 10 stochastic simulations, with and without scaling the input parameters. In Fig. 6 the mean values of H c are shown and compared to the reversal field calculated by conventional, non-thermal LLG. Without thermal fluctuations, a mesh size of about half the exchange length is required for accurate prediction of the interplay between exchange interaction and the demagnetizing field at corners and edges 20 . This value is 1.4 nm for the intrinsic material properties used in this work. Figure 6 clearly shows, that for this mesh size scaling is required in order to compute coercive fields that are within the range predicted by sweep rate-dependent theory 21 . The anisotropy field for the exchange-interaction-only model has been calculated in (2). With the material parameters of our cube, this comes to μ 0 H ani = 6.71 T. When the magnetostatic field is incorporated, non-thermal LLG simulations show that the coercive field is reduced by 0.35 T, to μ 0 H c = 6.36 T.
By adding thermal fluctuations, the energy barrier against switching the cube can be overcome with lower external fields. This physics can be captured as a temperature-and sweep-ratedependent effective coercive field H c . Various analytical expressions have been derived to describe this temperature-and sweeprate-dependent H c , where an applied field is swept at the constant rate 22 v. El-Hilo et al. 21 proposed The energy barrier of the cube at zero field E 0 = 522 k B 300 K was determined using the nudged elastic band method 23 and μ 0 H 0 = 6.36 T is the switching field without thermal fluctuations, but with magnetostatic fields. f 0 is the attempt frequency in zero external field, that strongly depends on the reversal mode of the sample.
The approximation of f 0 6 for homogeneous reversal of the cube with volume V, gives an attempt frequency of f 0 = 198 GHz. However, recent studies of inhomogeneous reversal 24 have suggested attempt frequencies exceeding 6 THz for hard magnetic single particles. For these two values for f 0 , we obtain the theoretical coercive fields of μ 0 H c = 6.02 T and 5.74 T, respectively. We use these two fields to define a range of expected coercive fields, signaled by the gray band in Fig. 6. We first calculated the effective coercive field with the unscaled sLLG. First, the coercive field had a substantial, unphysical variation with changing mesh size. Second, in the entire range of mesh sizes, the coercive field H c was well outside the gray band of theoretical expectations. Third, in fact, for the smallest mesh size of l = 0.5 nm, H c was 4.4 T, 23% less than the lowest edge of the expected band, a substantial discrepancy.
Finally, we simulated the same reversal with FUSSS LLG, using the same scaling as in the earlier test examples. As visible in Fig. 6, first, the mesh size dependence of H c was largely eliminated. Second, H c remained inside the gray band of theoretical expectations nearly the entire mesh size range. Both these facts are encouraging signs that FUSSS LLG introduces a marked, quantitative improvement over existing sLLG methods.

DISCUSSION
In this paper, we addressed the problem that standard stochastic Landau-Lifshitz-Gilbert (sLLG) simulations typically produce results that show unphysical mesh-size dependence. We identified the root cause of this problem: the effects of spin-wave fluctuations are ignored in sLLG. We proposed to represent the effect of these spin-wave fluctuations by a full-spin-wave-scaled stochastic LLG, or FUSSS LLG method. In FUSSS LLG that uses a mesh size l, the intrinsic parameters of the sLLG simulations M s , A x , and K u are first scaled by scaling factors s M (l), s A (l) = s M (l) 2 , and s K (l) = s M (l) b (with b = 3) that integrate out the spin-wave fluctuations up to the mesh size l, and the sLLG simulation is then performed with these scaled parameters. The scaling can be naturally anchored in the microscopic, or atomistic parameters. However, given that for hard magnets there is no consensus about these values, we chose to anchor our FUSSS LLG on the macroscopic scale, in experimentally measured quantities. In this sense, our scaling factors s M (l), s A (l), and s K (l) were integrating spinwave fluctuations back in to get the effective parameters right on the length scale of the mesh size l. Fig. 6 Coercive field test. Coercive field of a hard magnetic 10 nm cube calculated with a field sweep rate of v = 250 mT ns −1 . The anisotropy field is μ 0 H ani = 6.71 T. Conventional, non-thermal LLG gives a mesh independent coercive field of 6.36 T, reduced from the anisotropy field by the magnetostatic fields. Unscaled sLLG produces a coercive field heavily dependent on the mesh size, and outside the gray band of theoretical expectations in the entire mesh-size range. In contrast, FUSSS LLG greatly reduces the meshsize dependence and produces a coercive field that is inside the gray band of expectations for most mesh sizes. The error bars show the standard deviation, and the symbols show the mean of ten independent stochastic calculations.
We developed FUSSS LLG by studying the Ferromagnetic Resonance (FMR) in Nd 2 Fe 14 B cubes. We found that while the scaling greatly reduced the mesh size dependence relative to sLLG, full mesh size independence was not achieved with the nominal anisotropy scaling exponent of b = 3. However, we discovered that adjusting b from b = 3 to b = 2.72, a <10% adjustment, delivered fully mesh size independent results for the FMR peak. We note, that the universality of 2.72 as exponent needs to be further explored. To apply our method to other problems, we recommend conducting introductory pilot studies with small samples to determine the exponent and using this value for the actual large-scale study of interest.
We then performed three tests and validations of our FUSSS LLG with this modified scaling. First, we studied the same FMR but with magnetostatic fields included. This model exhibited two FMR peaks instead of one. Second, we simulated the total magnetization of the Nd 2 Fe 14 B cube. Third, we studied the effective, temperature-and sweeping rate-dependent coercive field of the cubes. In all three cases, we found that FUSSS LLG delivered essentially mesh-size-independent results, which tracked the theoretical expectations better than unscaled sLLG.
Before closing, we remark that a sister version of the present FUSSS LLG can be developed for magnets where the microscopic parameters are subject of agreement by the community, in which case the microscopically anchored FUSSS LLG can be an equally powerful method.
In sum, motivated by the success of our tests and validations, we propose that FUSSS LLG provides marked, qualitative progress towards accurate, high precision modeling of micromagnetics in hard, permanent magnets, magnetic recording media, and magnetic storage elements. In magnetic data storage, magnetic switching is required in subnanosecond time scale. Therefore, for the design and development of magnetic storage systems, the method proposed in this work can be directly applied. Proper scaling of the intrinsic parameters is also important for the simulation of magnetization reversal at longer time scales in which the solution of the sLLG equation is the main building block. For example in the forward flux sampling method applied by Vogler et al. 24,25 , many dynamic trajectories are computed by the solution of the sLLG equation, in order to simulate magnetization reversal on a time scale of many years.

METHODS
To calculate magnetization dynamics at finite temperatures we use the Langevin-type stochastic LLG equation. We chose to compute the magnetic state evolution of a single cube of Nd 2 Fe 14 B. We develop and analyze the spin-wave renormalized sLLG method by simulating two test scenarios, and using a series of mesh sizes for each scenario to test the mesh size independence of the results.
1. The equilibrium state of the cube was simulated by initially saturating the magnet in the easy-axis direction and let the cube relax without any external influence. The thermal perturbation in form of the stochastic field caused increasing deviations from the easy direction. The mean value of the magnetization's z-component over time 〈M z 〉 was calculated for various mesh sizes to quantify the mesh dependency. 2. The cube was again saturated in the easy-axis direction. Then it was exposed to an increasing opposing field which eventually reversed the magnet's magnetization at the coercive field H c . The switching simulations were performed with different mesh sizes to evaluate the possible mesh dependency of H c . In order to distinguish between the different length scales and their respective magnetic properties, we use superscripts a for atomistic length scale, l for the minimum scale of micromagnetic simulations, and L for macroscopic length scale as in measurements.

Stochastic LLG for finite-temperature micromagnetics
After discretization with the finite element method the Landau-Lifshitz-Gilbert (LLG) equation can be treated like the dynamics of interacting magnetic moments, whereby each moment is associated with one node of the finite element mesh 26 . The dynamics of the magnetic moment is given by the stochastic Landau-Lifshitz-Gilbert (sLLG) equation. The sLLG method starts with constructing the usual LLG equations and then adding to the effective field h eff,i a field representing thermal fluctuations h th,i acting at node i: Here m i is the magnetic moment at node i, the prefactor γ 0 ¼ γ j j=ð1 þ α 2 Þ contains the gyromagnetic ratio γ, and α is the Gilbert damping constant. Through finite element discretization the magnetic moment at node i can be expressed as m i = M l V i . Here, V i is the volume associated with node i and jM l j ¼ M l s . M l s is the spontaneous magnetization at the minimum scale of the micromagnetic simulation, the finite element mesh size l, typically a few nanometers. The total field is the sum of the effective field h eff and the stochastic field h th . The effective field is the sum of the exchange field h ex , the anisotropy field h ani , the demagnetizing field h demag and the applied field h ext .
The exchange field at node i can be derived from the continuous expression h ex ¼ 2A l x =ðM l s Þ 2 ∇ 2 M l using the Galerkin method.
The anisotropy field h ani; Here, u is the unit vector pointing in the easy-axis direction. The demagnetizing field h demag can be computed from a magnetic scalar potential and evaluated at the nodes by averaging over the neighboring elements.
The stochastic field h th introduces thermal fluctuations without correlation between the spatial components, in time or space. Therefore, it is a Gaussian random process that satisfies 27 h th;i ðtÞ ¼ 0; (8) and has the variance h th;i ðtÞ h th;j ðt þ ΔtÞ ¼ 2Dδ ij δðΔtÞ: This equation relates the strength of the thermal fluctuations to the dissipation of the system depending on the node volume V i . Here, k B = 1.38 × 10 −23 JK −1 is the Boltzmann constant, T is the temperature, and μ 0 = 4π × 10 −7 Hm −1 is the magnetic vacuum permeability.
To integrate the system of stochastic differential equations (7), we apply the midpoint scheme 28,29 . The approach outlined above is similar to that introduced by Ragusa and coworkers 8 who showed that the numerical solution converges to the analytic solution obtained by the Fokker-Planck equation.
A discretization cell represents a large number of fluctuating magnetic spins. The larger the node volume V i , the more spins are contributing, hence averaging effects reduce the total perturbation for this node. In the vicinity of atomistically small cells, no averaging over the thermal fluctuation is needed. In order to calculate magnetic states we interpret the stochastic equation in the Stratonovich sense and integrate (7) using a midpoint scheme 28 with a time step Δt = 1 fs. We apply a fixed-point scheme to solve the nonlinear system of equations in each time step 30 . The discretized stochastic LLG model is similar to an atomistic model. In an atomistic model, the magnetic moments at the atom positions interact, while in the discretized model the magnetic moments at the node points are considered. Whereas in the atomistic model the intrinsic parameters can be clearly taken from first-principle calculations, it is not fully obvious how to derive A l x , M l s , and K l u , because the effective magnetic moment at a computational cell arises from the thermal average of all spins contained in the cell. By considering the spin-wave spectrum at the different length scales, we propose a method for determining the intrinsic magnetic properties at the micromagnetic length scale l.
FE-models of a cube with 40 nm edge length were prepared with various mesh sizes from 1 nm to 10 nm using Salome 31 and NETGEN 32 . A representation of a meshed model is shown in the inset of Fig. 5. The mesh is composed of uniformly sized tetrahedral elements, which is important to obtain reliable renormalization factors. Previously tested stochastic simulations proved to be very sensitive to small changes in element size of irregular meshes. To make progress, the sLLG simulations need the magnet's parameters, the magnetization M s , the exchange stiffness constant A x , the magnetocrystalline anisotropy constant K u , and the Gilbert damping α. Typically, sLLG approaches adopt parameters calculated by ab initio methods on the atomic scale and use them directly as the parameters of the finite elements. In simplistic terms, present sLLG approaches assume that the magnetization of the individual finite elements is fully saturated. As such, present sLLG simulations have not taken into account how the intra-mesh spin-wave fluctuations reduce the magnetization from its saturated value.

Spin-wave renormalization of parameters
In order to incorporate the effects of the intra-mesh, short wavelength spinwave fluctuations up to the mesh size l, Grinstein and Koch proposed to integrate out the fluctuations caused by spin waves with wavenumbers higher than π/l by renormalization group theory. Doing so yielded scaledependent effective parameters for the sLLG simulation 11 . The limitations of this approach were discussed earlier, including (1) the use of an idealized, gapless spin-wave dispersion; (2) the approximate, perturbative expansion in both the magnetic field h and in the anisotropy constant K u to only leading logarithmic order; and (3) the disregard of the crystalline symmetries of the actual material, typically resulting in a non-isotropic spin-wave dispersion. As also discussed, these three approximations can be justified in soft magnets in zero fields, close to the critical temperature, where long wavelength critical fluctuations dominate the physics. However, typical sLLG simulations are not performed close to criticality, in zero fields, and exclusively in soft magnets. For such typical simulations, the above approximations become questionable, especially for hard magnets with large K u values.
Motivated by these limitations, in this paper we propose a similarlyinspired, but a distinct method to incorporate spin-wave fluctuations. In this full-spin-wave-scaled sLLG method, we propose to integrate out the spin-wave fluctuations by using the full, preferably experimentally verified spin-wave spectrum that can have a gap, is not expanded in h and in K u , and reflects the symmetries of the crystal. It is natural to expect that our method will achieve superior accuracy of sLLG simulations for magnets with stronger crystalline anisotropy, in higher fields, having non-negligible crystalline structures, at temperatures well below criticality. Taken on face value, this program starts with the atomistic magnetization M a s , determined by ab initio calculation on the length scale of the unit cell a. The reduction of the magnetization by spin waves of wavenumber k, ΔM(k), is then integrated out with wavelengths sweeping between the atomistic scale a and the mesh size l to yield the length-scale dependent magnetization ΔMðkÞ dk: Using the value of magnetization change caused by a spin-wave leads to where the Bose-Einstein occupation factor n EðkÞ; T ð Þis and E(k) is the full-spin-wave spectrum that includes the anisotropy constant K u and the magnetic field h fully, not only in leading perturbative order. E(k) also reflects the discrete symmetries of the crystal, and thus can include the wavevector k in a trigonometric function instead of simply as k 2 . μ B is the Bohr magneton. The ratio of the renormalized magnetization to the non-renormalized magnetization will be referred to as the scaling function, or scaling factor, s M ðlÞ ¼ M l s =M a s . Extending the definition of the scaling function this way increases the precision with which the spin-wave fluctuations are accounted for because our method is not perturbative in the magnetic field h and in the anisotropy constant K u , and furthermore, it also incorporates the discrete symmetries of the crystal. On the other hand, the justification for keeping only this term becomes less compelling because strictly speaking, we are not keeping only the leading logarithmic terms of the standard renormalization group theory. However, the farther the simulated system is from criticality, the justification to keep only the leading logarithmic terms itself becomes less compelling anyway. Therefore, for systems away from criticality the calculation of the scaling function by retaining the non-perturbative spin-wave energy dispersion with the explicit crystalline symmetries becomes a net positive improvement.
The scaling of the other parameters can be constructed by using well-known relationships: the exchange stiffness constant scales 33 as A l x / ðM l s Þ 2 , and so A l x =A a x ¼ s M ðlÞ 2 . Further, the magnetocrystalline anisotropy constant scales according to Callen and Callen's power law 34,35 for uniaxial anisotropy K l u / ðM l s Þ 3 , and thus K l u =K a u ¼ s M ðlÞ 3 . Once all three atomistic/microscopic parameters, M a s , A a x , and K a u have been scaled to their effective values M l s , A l x , and K l u to incorporate the intramesh spin-wave fluctuations up to wavelength l by their scaling factors, the sLLG simulation can be performed with the mesh size l. This microscale anchored full-spin-wave-scaled sLLG should yield mesh size-independent results.
To implement the above steps, we started by consulting the literature for the ab initio parameters of permanent magnets of interest. For hard Nd 2 Fe 14 B magnets, Herbst summarized the results of several ab initio calculations 36 . Quite remarkably, the calculated values showed a substantial variation, often differing by a factor of 2 or more. For this reason, it was quite difficult to establish consensus values of the ab initio calculations for hard Nd 2 Fe 14 B magnets.
Forced by this situation, we looked for data supported by widespread agreement. We found this among the experimentally determined macroscopic intrinsic properties of Nd 2 Fe 14 B 35 , taken at T = 300 K. The intrinsic properties are usually measured on bulk samples with a methodology averaging over a huge number of atomistic spins. In the following, we denote such measured intrinsic magnetic properties with the superscript L. The saturation magnetization is widely agreed to be μ 0 M L s ¼ 1:61 T, the exchange stiffness constant A L x ¼ 7:7 pJm À1 , and the magnetocrystalline anisotropy constant K L u ¼ 4:3 ; MJm À3 . To build on widely accepted data, we propose that the full-spin-wavescaled sLLG method can be applied in the reverse direction as well. When the sample magnetization is measured experimentally far away from reversal processes that take place close to H c and involve nucleation and domain wall propagation, it is reasonable to assume that the entire difference between the experimentally measured magnetization and the magnetization on the scale of the mesh size is caused by spin-wave fluctuations. Therefore, it is possible to determine the effective magnetization on the length scale of the mesh size l by adding back the magnetization reduction caused by spin waves to the experimentally measured magnetization, by integrating the contribution of spin waves with wavenumbers between π/L and π/l, where L is the macroscopic system size. We call this approach the macroscale-anchored full-spinwave-scaled sLLG method. We start from the experimentally measured macroscopic magnetization M L s measured at length scale L, and integrate the spin-wave corrections with wavelengths larger than the mesh size l back in: ΔMðkÞ dk (14) The integral on the right-hand side of (14) can be approximated by building the limit for infinite system length, therefore the lower limit of the integral is taken as zero. An approximation for the microscopic magnetization is Here we approximated M 1 s with M L s . E(k) is the spin-wave spectrum of an anisotropic Heisenberg ferromagnet 37 in an external field h: where Z is the coordination number of the lattice and a i are the nearest neighbor vectors. Experimental measurements of the spin-wave spectrum are consistent with this form 35 . As outlined in the foundational parts, E(k) retained the magnetic field and the crystalline anisotropy in full, instead of perturbatively expanding in them. These factors induced a gap in the spectrum, which would have been incompatible with the standard renormalization group formulation. Based on the above, the scaling function for the magnetization in this macroscale-anchored full-spin-wave-scaled sLLG takes the form: Examples of the scaling factors, s M (l), s A (l), and s K (l) as functions of the mesh size l are shown in Fig. 3 and discussed there.
Once the full-spin-wave-scaled magnetization M l s , exchange constant A l x , and anisotropy K l u have been constructed, the sLLG method with mesh size l can be used to simulate experiments that involve not only spin waves but nucleation, domain wall propagation, ferromagnetic resonance (FMR), and other, non-trivial phenomena. This full-spin-wave-scaled sLLG method should deliver mesh size independent results, and thus should introduce a major step forward in the accuracy and utility of sLLG methods for simulating complex and challenging magnetization dynamics.

DATA AVAILABILITY
All data generated during and/or analyzed during this study are available from the corresponding author on reasonable request.

CODE AVAILABILITY
All data were computed using the commercial finite element software femme.