Adaptively time stepping the stochastic Landau-Lifshitz-Gilbert equation at nonzero temperature: implementation and validation in MuMax3

Thermal fluctuations play an increasingly important role in micromagnetic research relevant for various biomedical and other technological applications. Until now, it was deemed necessary to use a time stepping algorithm with a fixed time step in order to perform micromagnetic simulations at nonzero temperatures. However, Berkov and Gorn have shown that the drift term which generally appears when solving stochastic differential equations can only influence the length of the magnetization. This quantity is however fixed in the case of the stochastic Landau-Lifshitz-Gilbert equation. In this paper, we exploit this fact to straightforwardly extend existing high order solvers with an adaptive time stepping algorithm. We implemented the presented methods in the freely available GPU-accelerated micromagnetic software package MuMax3 and used it to extensively validate the presented methods. Next to the advantage of having control over the error tolerance, we report a twenty fold speedup without a loss of accuracy, when using the presented methods as compared to the hereto best practice of using Heun's solver with a small fixed time step.


I. INTRODUCTION
Micromagnetic simulations of systems at nonzero temperatures are an increasingly important tool to numerically investigate magnetic systems relevant for technological applications. Historically, the foundations for a description of thermal fluctuations in micromagnetic systems were laid by Brown when he investigated the thermal switching of single-domain particles 2,3 . Today, these particles are used in promising biomedical applications such as disease detection and tumor treatment [4][5][6][7] . In order for these applications to be successful, a full understanding of the particles' thermal switching is important. For example, many characterization procedures [8][9][10][11][12][13][14] require this knowledge to accurately determine the particles' properties. Also diagnostic particle imaging [15][16][17][18][19][20][21][22][23] , and therapeutic applications 24,25 rely on models of the particles' thermal switching. Currently, these models are often based on approximations that not always take into account that, e.g. the magnetization state in large particles can deviate from a uniform magnetization 26 , or the particles might interact with each other via the magnetostatic interaction 27 . In such cases, the analytical models do not accurately reflect the true magnetization dynamics of the particles, and one has to rely on numerical models, the most accurate of which are based on a micromagnetic approach 28-34 . Next to their relevance in magnetic nanoparticle research, thermal fluctuations also play an important role in (exchange-coupled) continuous magnetic systems. One a) Electronic mail: jonathan.leliaert@ugent.be technologically relevant example is domain wall motion through a magnetic nanostrip, proposed as the operating principle for the racetrack memory [35][36][37] and for logic devices [38][39][40][41][42] . Recently, even smaller magnetization structures, i.e. skyrmions, have been proposed in both memory 43 and logic devices 44 . As the information carriers in these devices become smaller, the influence of thermal fluctuations further grows in importance: at such small spatial scales, the thermal stability of the bits themselves starts to become a relevant research question 45 . At the same time, their thermal depinning becomes an inherently stochastic process 46 . When numerically investigating domain wall motion at low driving forces, the dynamics can only be captured by considering the interplay between thermal fluctuations, the disorder energy landscape of the material and the driving forces. The resulting motion is then called domain wall creep 47 . Until now, full micromagnetic simulations are still prohibitively expensive in all but the smallest of such systems 48 .
Thermal fluctuations are also critical to the design of magnetic storage elements. They do not only determine the data retention limit of any magnetic storage system, but can also influence the read and write process of a MRAM cell. So estimation of read and write errors requires stochastic micromagnetic modeling of the spin valve during the application of the spin-torque current. This very challenging because of large time scales that are involved 49 .
There exist different theoretical approaches, each with their respective advantages and disadvantages, to study thermally induced magnetization dynamics 32,50 . Following Brown 3 , it is possible to derive the Fokker-Planck equation describing the time-dependent probability distribution of the magnetization directions of an ensem-ble of uniformly magnetized magnetic nanoparticles 3,32 . However, only in the simplest cases, e.g. when the particles' anisotropy axes are aligned with the applied field, an analytical solution can be found. In more complex cases approximations have to be introduced and when considering continuous systems consisting of several exchangecoupled finite difference cells, this approach becomes intractable.
Alternatively, thermal fluctuations can be included as a stochastic term in the Landau-Lifshitz-Gilbert (LLG) equation, henceforth named stochastic LLG (sLLG) equation by adding a thermal field term to the effective field. This approach was presented by Lyberatos 51 and is based on the fact that finite difference cells can be considered as dipoles comparable to single-domain particles. This method has the advantage that it is a straightforward extension of the LLG equation. On the other hand, in contrast to the Fokker-Planck approach, the resulting sLLG equation has to be solved many times in order to gather enough data to draw conclusions about averaged quantities.
Integrating stochastic differential equations (SDE) requires the use of non-Riemann calculus to be able to deal with the discontinuous thermal field term. A full discussion of this topic lies beyond the scope of this article, and for an excellent introduction we refer to Ref. 32 . In general, SDE's are not trivial to numerically integrate, and require specialized methods 52 , which are only suited to integrate SDE's written in either their Ito or Stratonovich form, as otherwise a drift term might appear in the solution. When considering variable time stepping algorithms, the complexity further increases 53,54 , sometimes even canceling the relative advantage obtained by using such methods. However, Berkov and Gorn have shown that the drift term in the sLLG equation manifests itself only in the length of the magnetization vector 1 which, in contrast to the Landau-Lifshitz-Bloch equation 55 , is held constant. As shown in Ref. 1 , this can be seen more clearly when writing the sLLG equation in spherical coordinates. In Cartesian coordinates, numerical noise would build up in this direction if it weren't accounted for by renormalizing the magnetization after each time step to ensure that the drift term does not influence the magnetization dynamics. Consequently, the Ito and Stratonovich interpretation are equivalent for integrating the sLLG equation, enabling the use of higher order solvers to integrate the sLLG eqation 56 . Despite this result, in literature one often still finds the recommendation to use the second order Heun's solver (here denoted with "RK12", because it is a second order Runge-Kutta type solver with embedded first order solution) with a very small fixed time step of the order of femtoseconds to simulate micromagnetic systems at finite temperatures [57][58][59] . In this paper, we present the use of higher order solvers with adaptive time stepping for such simulations and show that this method offers significant advantages.

II. METHODS
The Landau-Lifshitz-Gilbert (LLG) equation 60,61 contains a precession and a damping term, and describes the magnetization dynamics at the nanometer length scale and the picosecond timescale.
In this equation, γ denotes the gyromagnetic ratio, α the dimensionless Gilbert damping parameter and B eff the effective field. There is more than one way to add thermal fluctuations to the LLG equation. We choose to add a stochastic thermal field B therm as a contribution to the effective field term in both the precession and damping term, although it has been shown that the field contribution could also be omitted in the damping term if the size of the thermal field is rescaled adequately 58 . A second common option 32 is to add thermal fluctuation directly as an extra thermal torque term to the LLG equation.
The properties of the thermal field B therm were determined by Brown when he investigated the thermal switching of single-domain particles 3 . Later, it was realized that this theory was also applicable to micromagnetic simulations as each finite difference cell can be considered as such a particle 51 . The thermal field is given by Here, the operator · denotes a time average, ·· a correlation, δ the Dirac delta function and the indices i and j run over the x, y and z axes in a Cartesian coordinate system. The thermal field has zero average [Eq.
(3)] and its size q is given by Eq. (4). In this equation, k B denotes the Boltzmann constant, T the temperature, M s the saturation magnetization, and V the volume on which the thermal fluctuations act, i.e. the volume of a single finite difference cell. Equations (2) to (4) are determined such that the effect of the thermal fluctuations is independent of the spatial discretization used: when splitting up a volume into subvolumes and averaging the thermal fluctuations within those, one will recover the same resulting dynamics as in the undivided volume. The same is also true for the time step ∆t: when averaged out over a larger time, thermal fluctuations decrease in strength and again, this proportionality is determined such that the average dynamics do not depend on the time discretization.

A. Implementation in MuMax 3
MuMax 3 is a GPU-accelerated micromagnetic software package which numerically solves the LLG equation using a finite-difference discretization 59 . The thermal field is included in the effective field as where ∆t denotes the time step and η is a random vector drawn from a standard normal distribution whose value is redetermined after every time step 57,59 .
MuMax 3 provides several explicit Runge-Kutta methods to time step the LLG equation, the details of which can be found in Ref. 59 . Here, we will only mention the ones relevant for this work. Previously, simulations at nonzero temperatures in MuMax 3 were performed with the widely used Heun's method (RK12), using a very small time step of the order of 5 fs. In Section III, we will validate our results by comparing them to the solution obtained with this solver when there are no analytical solutions available.
For dynamical simulations at zero temperature, the default solver is the Dormand-Prince method (RK45). This solver offers 5th order error convergence and contains an embedded 4th order method to estimate the error. Generally speaking (i.e. when not investigating very fast dynamics, or in the absence of thermal fluctuations), it is not advantageous to implement even higherorder solvers. The reason for this is threefold: 1) For the moderately small torques encountered in typical micromagnetic simulations, the performance of the solver is limited by its stability regime. This means that using even slightly larger time steps will result in much larger errors no matter how small the exerted torques are. 2) Higher order methods typically need more intermediate torque evaluations per time step, thus reducing the advantage obtained by taking a larger time step. 3) Due to the memory required to store the results of the intermediate torque evaluations, higher-order solvers disproportionately increase the memory consumption compared to the obtained gain in performance.

B. Sixth order Runge-Kutta-Fehlberg solver
Due to the large size of the stochastic thermal field, simulations at nonzero temperatures require much smaller time steps, so that the solver performance is not longer limited by its stability regime, and large enough time steps can be taken to justify the additional intermediate torque evaluations. Therefore, we also implemented the 6th order Runge-Kutta-Fehlberg (RKF56) method with 5th order embedded solution shown in Table I. Unlike the RKF56 solver, some of the solvers used in MuMax 3 like the RK45 method, benefit from the firstsame-as-last (FSAL) property. In these solvers, the last torque evaluation of the current time step corresponds to the first evaluation of the next time step, thus effectively reducing the number of evaluations per step by 1.  However, as the stochastic thermal field is not constant in between time steps, the torque continuity requirement is not longer fulfilled, and the first and last torque evaluations have to be performed separately. Because the RKF56 solver never has the FSAL property, its performance can only compete with these other methods at nonzero temperatures, where the other methods do not benefit from the FSAL property either.
The implementation of the Sixth order Runge-Kutta-Fehlberg solver is subject to the same tests used for the other solvers implemented in MuMax 359 , and a full report is beyond the scope of this article. Nonetheless, Fig. 1 shows that our solution to standard problem 4, proposed by the µMag modeling group 63 , agrees with the solution obtained with OOMMF 64 (not using the RK56 solver).
One of the most sensitive checks one can perform to test the implementation of a solver is to investigate the scaling of its error convergence. Figure 2 shows the error after a single precession without damping of a single spin in a field of 0.1 T as function of the time step ∆t. The solver shows the expected sixth order convergence up to the limit of the single precision implementation 59 ( ≈ 10 −7 ).

C. Time stepping with adaptive time steps
When performing simulations at nonzero temperatures, it is important to note that the size of the thermal field is determined by 1/ √ ∆t. This implies that, when a large thermal field is generated leading to a bad step (defined as a step where the torque was too large for the used time step), the step will be undone, and the adaptive time step algorithm will decrease the time step, thus further increasing the size of the field. Luckily, this 1/ √ ∆t dependency makes the time step smaller at a slower rate than that the error is reduced (∆t N with N the order of the solver 65 ). However, it is important to use higher order solvers like the RK45 or RKF56 method in order to maintain a large time step.  To give the correct solution, the statistical properties of the random numbers η in Eq.(5) should correspond to the ones determined in Eqs.
(2) to (4). However, if one would redraw these random numbers after a bad step, the small thermal fields (small η) would be applied during longer time steps and large thermal fields (large η) during shorter time steps, thus virtually changing the distribution of the random numbers, and eventually giving rise to incorrect solutions. In our implementation we avoid this by keeping the previously drawn random numbers and rescaling the thermal field with a factor ∆t new /∆t badstep in case a bad step is encountered to ensure that the correct statistical properties of the thermal field are maintained.

III. VALIDATION
The adaptive time stepping at nonzero temperatures will be tested in several test cases, focusing on different aspects, i.e. static vs. dynamic properties of uncoupled spins or continuous magnets. Each time, the simulation results obtained with adaptive time stepping will be compared either to analytical solutions or to the solutions obtained with the RK12 method with fixed time step.
A. Spectra of a single spin It will be verified whether the thermal field and the resulting magnetization dynamics of a single spin in the absence of an external field shows the theoretically expected behavior. The thermal field should display a white spectrum S H construction, as its size is given by Gaussian random numbers. Figure 3 proves that this is indeed the case. The random thermal fields acting on a single isotropic finite-difference cell gives rise to a random walk on the unit sphere 11,32 . The shape of the spectral density S M (f ) of such a random walk is described by the square root of a Lorentzian 66 , i.e. white noise with a 1/f cutoff at a cutoff frequency f 0 given by 33,58 Figure 4 shows the obtained magnetization spectra (gray dots) indeed coincide with the red lines determined by Eqs. (6) and (7). Because all spectra coincide, they were rescaled with an arbitrary factor for clarity. The magnetization spectra of a single spin (rescaled with an arbitrary factor for clarity reasons) obtained with the RK45 solver with fixed time steps, and with the RK45 and RKF56 solver with adaptive time steps. All spectra display the theoretically expected shape depicted by the red lines.

B. Equilibrium magnetization
This validation problem checks whether the magnetization of an ensemble of uncoupled spins in thermal equilibrium in an externally applied field is described by the Langevin function L(ξ), where the argument ξ stands for Figure 5 proves that this is indeed the case for 4 different ξ (realized by 2 different cell sizes and 2 different temperatures) simulated with the RK45 and RKF56 solver with adaptive time steps over a large range of applied fields.

C. Thermal switching
After the equilibrium magnetization addressed in the previous problem, we will now concern ourselves with a dynamical problem consisting of the thermal switching rate of a single (macro-)spin particle with uniaxial anisotropy. In the limit of a high energy barrier (compared to the thermal energy), the switching rate ν is given by 67 For an ensemble of uncoupled spins, initialized with all spins pointing in the same direction, this switching gives rise to an exponentially decaying magnetization, with a decay constant 1/2ν. In this test problem, we simulate this decay to determine the numerical switching rate, and compare these values to their theoretical prediction. Figure 6 shows Arrhenius plots for the temperaturedependent switching rate ν of uncoupled finite difference cells with volume V =(10 nm) 3 and uniaxial anisotropy constant K=1×10 4 or 2×10 4 J/m 3 . Again, a quantitative agreement is seen between the MuMax 3 simulations and the theoretically predicted behavior. or 20 kJ/m 3 . Simulations were performed using the RK12 solver with fixed time steps (∆t=5 fs) or the RK45 or RKF56 solver with adaptive time steps on an ensemble of 2 18 noninteracting cells for 1 µs or until the ensemble magnetization crossed 0. All results agree with the red solid lines depicting the analytically expected switching rates (Eq. 10).

D. Thermally excited magnetization spectrum
In this problem we look at the thermally excited magnetization spectrum of a 10 nm thick disk with a diameter of 512 nm, M s =1 MA/m, exchange constant A ex =10 pJ/m, and α = 1, discretized in cells measuring 4 by 4 by 10 nm 3 . The equilibrium magnetization structure in such a disk is a vortex structure, as depicted in the inset of Fig. 7. We apply a thermal field corresponding to 300 K, thereby thermally exciting the sample, resulting in the spectra shown in Fig. 7. The spectrum depicted in red was obtained with the RK12 solver with fixed time step (∆t=5 fs), and serves as a reference solution. The gray lines, which agree almost perfectly with the benchmark solution, correspond to the three spectra obtained with the RK45 solver with the time step fixed to ∆t=300 fs and with the RK45 and RKF56s solvers with adaptive time step and = 10 −5 . Note that the spectra overlap even at high frequencies, indicating that the adaptive time stepping does not lead to spectral leakages. The red line corresponds to the spectrum obtained with the Heun solver with ∆t=5 fs. The three gray lines, which overlap almost perfectly with the red one, correspond to the RK45 solver with fixed time step ∆t=300 fs, and with the solution obtained with the RK45 and RKF56 solver with adaptive time stepping with = 10 −5 . All spectra were averaged out over 25 realizations with a different random seed for the thermal field. The inset shows the equilibrium vortex magnetization structure in the system under consideration a a The minimum in the error is a result between a direct match between the used time step and the gyration period used for the evaluation of the precision 59 .

E. Thermal diffusion of a domain wall
In a last validation problem we investigate the thermal driven diffusion of transverse domain walls in a non-disordered permalloy nanowire 68 . We simulate a nanowire with cross-sectional dimensions of 100×10 nm 2 discretized in cells of 3.125 × 3.125 × 10nm 3 . We use the material parameters of permalloy: M s = 860 kA/m, Aex = 13 pJ/m, α=0.01 and simulate the domain wall in the center of a moving window, as shown in Fig. 8. The thermally driven motion can be described by a random walk characterized by a diffusion constant D which scales linearly with temperature. Similarly as in Ref. 68 we simulate a transverse domain wall for 100 ns and repeat this simulation with a large number of different random seeds. In Tab. II, we compare the results obtained from the full micromagnetic simulations with the diffusion constant D of approximately 310 nm 2 /ns predicted by the model introduced, and numerically validated in Ref. 68 using the RK12 solver with fixed time step. The results show that the standard errors s are larger than the difference between the obtained diffusion constants and the expected value. This also indicates that, for this problem, an error tolerance = 10 −3 suffices for all practical purposes, as the variance between different simulations gives rise to an uncertainty that is larger than the errors due to the use of this relatively large . As this is problem dependent, the default value in MuMax 3 remains 10 −5 , but we suggest that, depending on the system under consideration, larger values might be suitable in simulations at nonzero temperatures.

IV. PERFORMANCE
To assess the performance of the presented methods we will consider the problems detailed in III D and III E, i.e the thermal spectrum of a disk and thermally driven domain wall diffusion, respectively, as benchmarks. The benchmark results are shown in Fig. 9 a) and b). For each of these problems, the simulation time was first determined when solving the problem with the RK12 solver with a fixed time step of 5 fs. As this is a second order solver with embedded first order solver, the difference between both solutions serves as an error estimate, allowing us to estimate with which error tolerance in the adaptive time stepping the results should be compared. This data point is indicated by a black cross in the figure. Next, the problems were solved using adaptive time stepping, once with the RK45, and once with the RKF56 solver, with ranging from 10 −3 to 10 −7 , shown in gray and red, respectively. The simulation runtime [Figs. 9 a) and b)] show the time it took to simulate 10 ns of magnetization dynamics, while Figs. 9 c) and d) indicate the average time step ∆t used by the solver, for each . When comparing the results from the adaptive time stepping methods with the result of the RK12 solver at the same estimated , one sees that the adaptive time stepping methods use a considerably larger time step without a loss of accuracy. For these two benchmark problems, this results in a twenty fold speedup of the simulation.
We also investigated the performance of the system described in Section III C, i.e. thermal switching of uncoupled spins. There, an even higher speedup was achieved, but this was attributed to the fact that such systems do not require the calculation of the demagnetizing field so that other factors, like the generation of the random numbers for the thermal field, become the limiting factor. Because the random numbers have to be generated only once per time step independently of the used solver, a very large performance gain can be achieved by using higher order solvers which allow time steps that are over a 1000 times larger than the ones necessary for the RK12 method with the same accuracy. However, as this highly depends on the used hardware, the performance gained by using the adaptive time stepping methods can lie anywhere between a factor 20 to 10 000, depending on accuracy. These observations indicate that the presented methods are particularly suited for magnetic nanoparticle research, where micromagnetic simulations are becoming increasingly important 34 .
The relative performance of the RK45 and RKF56 solver is comparable and show the expected trends that the RK45 solver is faster at higher and simulated temperatures while the RKF56 solver is faster at low and temperatures. Generally, only when simulating systems at very high temperatures, or with very small , it does pay off to use the RKF56 solver.

V. CONCLUSIONS
In this paper, we have exploited the fact that the drift term in the stochastic Landau-Lifshitz-Gilbert equation is only able to manifests itself in the direction of the magnetization length, which is fixed. Therefore, we were able to straightforwardly extend existing high order solvers with adaptive time stepping at nonzero temperatures. In an effort to further increase the performance, we have implemented the sixth order Runge-Kutta-Fehlberg solver, The top row shows the runtime required to simulate 10 ns of magnetization dynamics, while the bottom row shows the average time step used. The black crosses indicate the performance of the fixed time step RK12 method with estimated from the difference between the second order and embedded first order solution. All results were obtained an NVIDIA Titan Xp GPU in a system running on a 7th generation i5 CPU.
and we extensively validated both the correctness of this newly implemented solver and the adaptive time stepping method used at nonzero temperature. All presented methods are included in the open-source micromagnetic software package MuMax 3 and are thus freely available online.
The main advantages of the presented adaptive time stepping methods at nonzero temperatures are that they offer an inherent error control, which is unavailable with fixed time stepping methods, and without a loss of accuracy one can obtain a twenty fold speedup compared to the commonly best practice of using the RK12 solver with small fixed time step. This enables simulations which previously took too long to be considered feasible and will be useful for micromagnetic research of continuous (exchange coupled) systems like spin valves, or domain wall motion in nanowires, and for uncoupled spins, e.g. in magnetic nanoparticle research.