# Numerical integration of the stochastic Landau-Lifshitz-Gilbert equation in generic time-discretization schemes

###### Abstract

We introduce a numerical method to integrate the stochastic Landau-Lifshitz-Gilbert equation in spherical coordinates for generic discretization schemes. This method conserves the magnetization modulus and ensures the approach to equilibrium under the expected conditions. We test the algorithm on a benchmark problem: the dynamics of a uniformly magnetized ellipsoid. We investigate the influence of various parameters, and in particular, we analyze the efficiency of the numerical integration, in terms of the number of steps needed to reach a chosen long-time with a given accuracy.

## I Introduction

The design of magnetic devices used to store and process information crucially relies on a detailed understanding of how the magnetization dynamics are influenced not only by external magnetic fields but also by dissipation and thermal fluctuations Hillebrands2002 (); Bertotti09 (). In the simplest scenario, the time evolution of the magnetization is governed by the stochastic generalization of the Landau-Lifshitz-Gilbert (LLG) equation introduced by Brown to study the relaxation of ferromagnetic nanoparticles Brown (). In recent years, much attention has been directed to the theoretical and experimental understanding of how the magnetization can be manipulated with spin polarized currents via the spin torque effect originally discussed by Slonczewski Slonczewski1996 () and Berger Berger1996 (), an effect that can be described by a simple generalization of this equation.

Explicit analytical solutions to the stochastic LLG equation are available in very few cases; in more general circumstances information has to be obtained by direct numerical simulation of the stochastic equation, the study of the associated Fokker-Planck equation (see Ref. review, for a recent review), or via functional methods Aron14 ().

The stochastic LLG equation is a stochastic equation with multiplicative noise. It is a well-known fact that in these cases, a careful analysis of the stochastic integration prescriptions is needed to preserve the physical properties of the model. In the stochastic LLG case one should force the modulus of the magnetization to stay constant during evolution, and different schemes (Ito, Stratonovich or the generic ‘alpha’ prescription) require the addition of different drift terms to preserve this property (for a recent discussion see Ref. Aron14, ). All these issues are by now well-understood and they are also easy to implement in the continuous time treatment of the problem. Nevertheless, this problem has not been analyzed in as much detail in the numerical formulation of the equation.

Indeed, most of the works focusing on the numerical analysis of the stochastic equation use Cartesian coordinates Garcia98 (); Nowak2000 (); Li2004 (); Cheng2005 (); Cheng2006 (); dAquino (); Usadel2006 (); Kazantseva2008 (); Titov2010 (); Weiler2011 (); Haase2012 (). Although there is nothing fundamentally wrong with this coordinate system, most algorithms based on it do not preserve, in an automatic way, the norm of the magnetization during time evolution. These algorithms require the explicit magnetization normalization after every time step, a trick that is often hidden behind other technical difficulties Martinez2004 (); Cimrak2007 (). This problem can be avoided only if the specific midpoint prescription (Stratonovich) is used dAquino ().

Given that the modulus of the magnetization should be constant by construction, a more convenient way to describe the time evolution should be to use the spherical coordinate system. Despite its naturalness, no detailed analysis of this case exists in the literature. The aim of this work is to present a numerical algorithm to solve the LLG equation in the spherical coordinates system and to discuss in detail how different discretization prescriptions are related, an issue which is not trivial due to the multiplicative character of the thermal noise.

In order to make precise statements, we focus on the study of the low-temperature dynamics of an ellipsoidal Cobalt nanoparticle, a system that has been previously studied in great detail by other groups dAquino (). Our goal is to introduce the numerical method in the simplest possible setting and the uniaxial symmetric potential involved in this problem seems to us a very good choice.

The paper is organized as follows. In Sec. II we present the problem. We first recall the stochastic LLG equation in Cartesian and spherical coordinates. In both cases we discuss the drift term needed to ensure the conservation of the magnetization modulus as well as the approach to Boltzmann equilibrium. We then describe the concrete problem that we solve numerically. In Sec. III we present the numerical analysis. We first introduce the algorithm and then discuss the results. Section IV is devoted to the conclusions.

## Ii The problem

### ii.1 The stochastic Landau-Lifshitz-Gilbert equation

The stochastic Landau-Lifshitz-Gilbert (sLLG) equation in the Landau formulation of dissipation Landau () reads

(1) | |||||

where . is the product of , the gyromagnetic ratio relating the magnetization to the angular momentum, and , the vacuum permeability constant. The gyromagnetic factor is given by and in our convention with Bohr’s magneton and Lande’s -factor. The symbol denotes a vector product. For the first term in the right-hand-side describes the magnetization precession around the local effective magnetic field . The term proportional to is responsible for dissipation. Thermal effects are introduced à la Brown via the random field Brown () which is assumed to be Gaussian distributed with average and correlations

(2) |

for all . The parameter is, for the moment, free and is determined below. is the dissipation coefficient and in most relevant physical applications, . An equivalent way of introducing dissipation was proposed by Gilbert Gilbert () but we have chosen to work with the Landau formalism in this work.

This equation conserves the modulus of and takes to Boltzmann equilibrium only if the Stratonovich, mid-point prescription, stochastic calculus is used. Otherwise, for other stochastic discretization prescriptions, none of these physically expected properties are ensured. The addition of a carefully chosen drift term is needed to recover the validity of these properties when other stochastic calculi are used. The generic modified sLLG equation Aron14 ()

(3) | |||||

where the time-derivative has been replaced by the -covariant derivative

(4) |

ensures the conservation of the magnetization modulus and convergence to Boltzmann equilibrium for any value of . The reason for the need of an extra term in the covariant derivative is that the chain-rule for time-derivatives of functions of the stochastic variable is not the usual one when generic stochastic calculus is used. It involves an additional term (for a detailed explanation see Ref. Aron14, ). In addition, having modified the stochastic equation in this way, one easily proves that the associated Fokker-Planck equation is independent of and takes the magnetization to its equilibrium Boltzmann distribution at temperature provided the parameter is given by

(5) |

where is the volume of the sample that behaves as a single macrospin, the Boltzmann constant, and the saturation magnetization. The parameter is constrained to vary in . The most popular conventions are the Ito one that corresponds to and the Stratonovich calculus which is defined by . Note that this is not in contradiction with the claim by García-Palacios Garcia-Palacios () that Ito calculus does not take his LLG equation to Boltzmann equilibrium as he keeps, as the starting point, the same form for the Ito and Stratonovich calculations and, consequently, he obtains the Boltzmann result only for Stratonovich rules.

As the modulus of the magnetization is conserved, this problem admits a more natural representation in spherical coordinates. The vector defines the usual local basis () with

(6) |

and

(7) | |||||

The sLLG equation in this system of coordinates becomes Aron14 ()

(8) | |||||

(9) | |||||

(10) | |||||

where the and components of the stochastic field are defined as

(11) | |||||

(12) |

and similarly for .

We introduce an adimensional time, , and the adimensional damping constant , and we normalize the field and the magnetization by defining, , , , to write the equations as

(13) | |||||

(14) |

The random field statistics is now modified to and .

### ii.2 The benchmark

We focus here on the dynamics of a uniformly magnetized ellipsoid with energy per unit volume

(15) |

is the external magnetic field and are the anisotropy parameters. This case has been analized in detail in Ref. dAquino, and is used as a benchmark with which to compare our results. We normalize the energy density by , and write

(16) |

The effective magnetic field is . Once normalised by , it reads

(17) |

We study the dynamics of a Cobalt nanoparticle of prolate spherical form with radii [nm] (in the easy-axis direction) and [nm] (in the and directions, respectively), yielding a volume [m]. There is no external applied field, the saturation magnetization is [A/m], the uniaxial anisotropy constant in the direction is [J/m], and the temperature is [K]. In the following we work with the adimensional damping constant and the physical value for it is . For this nanoparticle one has (where are the demagnetization factors Cullity ()), and since . The constant takes the value [m/(As)]. We recall that in Ref. dAquino, the time-step used in the numerical integration is [ps], that is equivalent to .

For future reference, we mention here that in the absence of an external field, , the energy density can be written in terms of the component of the magnetization as

(18) |

Note also that in this system the anisotropy-energy barrier is , and therefore the ratio , that indicates that the dynamics take place in the low-temperature regime.

## Iii Numerical analysis

In this section we first give some details on the way in which we implemented the numerical code that integrates the equations, and we next present our results.

### iii.1 Method

First, we stress an important fact explained in Ref. Aron14, : the random fields and are not Gaussian white noises but acquire, due to the prefactors that depend on the angles, a more complex distribution function. Therefore, we do not draw these random numbers but the original Cartesian components of the random field which are uncorrelated Gaussian white noises. We then recover the field and by using Eqs. (11) and (12) and the time-discretization of the product explained below. Most methods used to integrate the sLLG equation rely on explicit schemes. Such are the cases of the Euler and Heun methods. While the former converges to the Ito solution, the latter leads to the Stratonovich limit Garcia98 (). To preserve the module of , in these algorithms it is necessary to normalize the magnetization in each step, a nonlinear modification of the original sLLG dynamics Cimrak2007 (). Implicit schemes, on the other hand, are very stable and, for example, the mid-point method (Stratonovich stochastic calculus) provides a simple way to automatically preserve the module under discretization dAquino (). In what follows, we describe our numerical-implicit scheme which keeps the module length constant and, unlike previous approaches, is valid for any discretization prescription.

Next, we define the -prescription angular variables according to

(19) | |||

(20) |

with . In the following we use the short-hand notation , , and so on. The discretized dynamic equations now read and with

(21) | |||||

(22) | |||||

where , the effective fields at the -point are and , and and . As we said above, we first draw the Cartesian components of the fields () as

(23) |

where the are Gaussian random numbers with mean zero and variance one, and we then calculate and using Eqs. (11) and (12).

The numerical integration of the discretised dynamics consists in finding the roots of the coupled system of equations and with the left-hand sides given in Eqs. (21) and (22). We used a Newton-Raphson routine Press () and we imposed that the quantity be smaller than . To avoid singular behavior when the magnetization gets too close to the axis, or , we apply in these cases a rotation of the coordinate system around the axis.

All the results we present below, averages and distributions, have been computed using independent runs.

### iii.2 Results

#### iii.2.1 Stratonovich calculus

We start by using the Stratonovich discretization scheme, , to numerically integrate the stochastic equation using the parameters listed in Sec. II.2 which are the same as the ones used in Ref. dAquino, . We simply stress here that these are typical parameters (in particular, note the small value of the damping coefficient ). Although we solved the problem in spherical coordinates, we illustrate our results in Cartesian coordinates [using Eqs. (7) to transform back to these coordinates] to allow for easier comparison with the existing literature.

##### Trajectories.

Figure 1 displays the three Cartesian components of the magnetization, , , and , as a function of time for a single run starting from an initial condition that is perfectly polarized along the axis, . The data show that while the and components fluctuate around zero, the component has telegraphic noise, due to the very fast magnetization reversal from the ‘up’ to the ‘down’ position and vice versa. Indeed, the working temperature we are using is rather low, but sufficient to drive such transitions.

##### Equilibrium criteria.

In Fig. 2 we show the relaxation of the thermal average of the component, , evolving from the totally polarized initial condition, , during a maximum adimensional time . In the inset one can see temporal fluctuations around zero in the averages of the other two components, and . The error bars in these and other plots are estimated as one standard deviation from the data average, and when these are smaller than the data points we do not include them in the plots. The data in Fig. 2 demonstrate that for times shorter than the system is still out of equilibrium while for longer times this average is very close to the equilibrium expectation, .

A more stringent test of equilibration is given by the analysis of the probability distribution function (pdf) of the three Cartesian components , , and . In Fig. 3 (a) we present the numerical pdf’s of , , and we compare the numerical data to the theoretical distribution function in equilibrium. We computed the former by sampling over the second half of the temporal window, that is, by constructing the histogram with data collected over , and then averaging the histograms over independent runs. For the equilibrium we note that the equilibrium probability density of the spherical angles is

(24) |

where , which implies

(25) | |||||

Here is the partition function and, for the parameters used in the simulation, and .

It is quite clear from Fig. 3 that the numerical curves for the two shortest are still far from the equilibrium one, having excessive weight on positive values of . The last curve, obtained for the longest running time, is, on the contrary, indistinguishable from the equilibrium one in this presentation. A more quantitative comparison between numerical and analytic pdf’s is given in the inset in Fig. 3 (a), where the probability distribution ‘H-function’ Kubo92 ()

(26) |

with given in Eq. (25), is plotted as a function of the inverse time . The two sets of data in the inset correspond to computed with the continuous analytic form (25), data falling above, and with a discretized version of it, where the same number of bins as in the numerical simulation is used (specifically, 51), and data falling below and getting very close to zero for the longest used. The latter is the correct way of comparing analytic and numerical data and yields, indeed, a better agreement with what was expected. Finally, in Fig. 3 (b) we show the pdf of for the same three used in Fig. 3 (a), and we observe a faster convergence to an equilibrium distribution with a form that is very close to a Gaussian. For symmetry reasons the behavior of is the same.

Figure 4 shows the pdf’s for the same set of parameters but starting from the initial condition . The approach to equilibrium is faster in this case: all curves fall on top of the theoretical one. Insets show the time dependence of and which still fluctuate around zero with larger temporal fluctuations for the latter than the former. We reckon here that the fluctuations of and are quite different. The oscillations of around zero are due to the telegraphic noise of this component and to the fact that the average is done over a finite number of runs. The amplitudes of these oscillations tend to zero with an increasing number of averages.

We conclude this analysis by stating that the dynamics in the spherical coordinate system for the Stratonovich discretization scheme behave correctly, with the advantage of keeping the norm of the magnetization fixed by definition.

#### iii.2.2 Generic calculus

Although it was shown in Ref. Aron14, that in the limit every discretization of the stochastic equation leads (at equilibrium) to the Boltzmann distribution, the numerical integration of the equations is done at finite and then both, the time-dependent and the equilibrium averaged observables may depend on . With this in mind, we investigated which discretization scheme is more efficient in terms of computational effort. The aim of this section is to study the dependence of the numerical results for different values of and to determine for which one can get closer to the continuous-time limit () for larger values of .

Figure 5 shows the temporal dependence of the for Stratonovich calculus, i.e. for , for a small window of time and from the initial condition . A very fast decay followed by a slow relaxation is observed. The phenomena can be well fitted by a sum of two exponential functions, one (of small amplitude) describing the rapid intra-well processes and the another one the dominant slow over-barrier thermo-activation Garanin1996 (). Concretely, we used

and we found that the best description of data is given by

that is and , consistently with statements in Ref. Garanin1996, . In addition, we compared our result for to the one arising from Eq. (3.6) in Ref. Garanin1996, for our parameters and we found

(27) |

which is less than our numerical estimate, a very reasonable agreement, in our opinion.

Most importantly, we reckon that the time-dependent results do not depend strongly on for (see Fig. 5, where data for , , and prove this claim) and we can assert that the master curve is as close as we can get, for the numerical accuracy we are interested in, to the one for , that is, to the correct relaxation. Instead, for other discretization prescriptions, the dependence on is stronger. For example, for (Ito calculus) the curves for and are still notably different from each other (see the insert to Fig. 5), and they have not yet converged to the physical time-dependent average. Even smaller values of are needed to get close to the asymptotically correct relaxation, shown by the dotted black line. We do not show the pdf’s here, but consistently, they are far away from the equilibrium one for these values of .

In Fig. 6 we use the initial condition to see whether the efficiency of the Ito calculus improves in this case. Although the values of are very close to the expected vanishing value both distributions, (main panel) and (inset), are still far from equilibrium. We conclude that also for this set of initial conditions smaller are needed to reach the continuous-time limit. We have investigated other values of and in all cases we have found that convergence is slower than for the case.

We conclude that the Stratonovich calculus is ‘more efficient’ than all other -prescriptions in the sense that one can safely use larger values of (and therefore reach longer times) in the simulation. This does not mean that other discretization schemes yield incorrect results. For one must use smaller values of the time-step to obtain the physical behavior.

#### iii.2.3 Effect of the damping coefficient

Previously, we studied the magnetization relaxation for a physically small damping coefficient, . Under these conditions relaxation is very slow and it is difficult to reach convergence for generic values of . To overcome this problem, we increased slightly the damping up to . As a consequence, and because we are still in the low damping regime, relaxation to equilibrium is expected to be faster [in Fig. 7 (a) we show below that this statement is correct]. Note that more subtle issues can arise in the non axially symmetric case if initial conditions are not properly chosen Kalmykov2010 (). Then, in this subsection we check whether the dependence found for the calculus improves under these new dissipation conditions.

In Fig. 7 (a) we test the dependence of for and five values of ranging from to and increasing by a factor of two. Filled and open data points of the same color correspond to and , respectively. The agreement between the two data sets is very good for all . Indeed the agreement is so good that the data are superimposed and it is hard to distinguish the different cases. The curves also show that the dynamics are faster for increasing . Figure 7 (b) displays the decay of as a function of time for and a rather large value of the damping coefficient, , for different time increments, . The curves tend to approach the reference one shown by the solid black line and corresponding to for decreasing values of . A quantitative measure of the convergence rate is given by another parameter, defined as

(28) |

and shown in the inset for . Here, is the average of for and , while is the curve corresponding to other values of and . For all -schemes tends to zero for . Note that for this parameter is very close to zero for all the values shown in Fig. 7 (b), confirmation of the fact that this prescription yields very good results for relatively large values of and it is therefore ‘more efficient’ computationally.

## Iv Conclusions

In this paper we have introduced a numerical algorithm that solves the the sLLG dynamic equation in the spherical coordinate system with no need for artificial normalization of the magnetization. We checked that the algorithm yields the correct evolution of a simple and well-documented problem dAquino (), the dynamics of an ellipsoidal magnetic nanoparticle. We applied the algorithm in the generic ‘alpha’-discretization prescription. We showed explicitly how the finite dynamics depend on , despite the fact that the final equilibrium distribution is -independent. We showed that at least for the case reported here, the Stratonovich mid-point prescription is the ‘more efficient one’ in the sense that the dependence of the dynamics on the finite value of is less pronounced so, larger values of can be used to explore the long time dynamics.

We think it would be worthwhile to explore, both analytically and numerically, if this is a generic result of the sLLG dynamics. A priori, it is not clear what will be the optimal prescription to deal with other problems such as a system under a non-zero longitudinal external magnetic field or for a more general non-axially symmetric potential Coffey1995 (); Garanin1996 (); Kalmykov2010 ().

Finally, we mention that it is well known that for a particle on a line with multiplicative noise, the addition of an inertial term acts as a regularization scheme that after the zero mass limit “selects” the Stratonovich prescription (see Ref. Aron10, ). For the case of magnetization dynamics non Markovian generalizations of the LLG equation have been considered in Refs. Miyasaki1998, ; Atxitia09, . It could be interesting to analyze in this case, how the markovian limit relates to any specific stochastic prescription. We plan to report on this issue in the near future.

###### Acknowledgements.

We thank C. Aron, D. Barci and Z. González-Arenas for very helpful discussions on this topic. F.R. acknowledges financial support from CONICET (Grand No. PIP 114-201001-00172) and Universidad Nacional de San Luis, Argentina (Grand No. PROIPRO 31712) and thanks the LPTHE for hospitality during the preparation of this work. L.F.C. and G.S.L. acknowledge financial support from FONCyT, Argentina (Grand No. PICT-2008-0516).## References

- (1) B. Hillebrands and K. Ounadjela, eds. Spin dynamics in confined magnetic structures (Springer, Berlin, 2002).
- (2) G. Bertotti, I. Mayergoyz, and C. Serpico, Nonlinear magnetization dynamics in nanosystems (Elsevier, Amsterdam, 2009).
- (3) W. F. Brown, Phys. Rev. 130, 1677 (1963).
- (4) J. C. Slonczewski, J. Magn. Magn. Mat. 159, L1 (1996).
- (5) L. Berger, Phys. Rev. B 54, 9353 (1996).
- (6) Thermal fluctuations of magnetic nanoparticles: Fifty years after Brown W. T. Coffey and Y. P. Kalmykov, J. App. Phys. 112, 121301 (2012).
- (7) C. Aron, D. G. Barci, L. F. Cugliandolo, Z. González Arenas, and G. S. Lozano, arXiv:1402.1200. In press in J. Stat. Mech. (2014).
- (8) J. L. García-Palacios and F. J. Lázaro, Phys. Rev. B 58, 14937 (1998).
- (9) U. Nowak, R. W. Chantrell, and E. C. Kennedy, Phys. Rev. Lett. 84, 163 (2000).
- (10) Z. Li and S. Zhang, Phys. Rev. B 69, 134416 (2004).
- (11) X. Z. Cheng, M. B. A. Jalil, H. K. Lee, and Y. Okabe, Phys. Rev. B 72, 094420 (2005).
- (12) X. Z. Cheng, M. B. A. Jalil, H. K. Lee, and Y. Okabe, Phys. Rev. Lett. 96, 067208 (2006).
- (13) M. d’Aquino, C. Serpico, G. Coppola, I. D. Mayergoyz, and G. Bertotti, J. of Appl. Phys. 99, 08B905 (2006).
- (14) K. D. Usadel, Phys. Rev. B 73, 212405 (2006).
- (15) N. Kazantseva, D. Hinzke, U. Nowak, R. W. Chantrell, U. Atxitia, and O. Chubykalo-Fesenko, Phys. Rev. B 77, 184428 (2008).
- (16) S. V. Titov and P. M. Déjardin, H. El Mrabti, and Y. P. Kalmykov, Phys. Rev. B 82, 100413 (2010).
- (17) M. Weiler, L. Dreher, C. Heeg, H. Huebl, R. Gross, M. S. Brandt, and S. T. B. Goennenwein, Phys. Rev. Lett. 106, 117601 (2011).
- (18) C. Haase and U. Nowak, Phys. Rev. B 85, 045435 (2012).
- (19) E. Martínez, L. López-Díaz, L. Torres, and O. Alejos, Physica B 343, 252 (2004).
- (20) I. Cimrák, Arch. Comput. Meth. Eng. 15, 1 (2007).
- (21) L. D. Landau and E. M. Lifshitz, Phys. Z. Sowjetunion 8, 153 (1935).
- (22) T. L. Gilbert, Phys. Rev. 100, 1243 (1955). T. L. Gilbert, IEEE Trans. Mag. 40, 3443 (2004).
- (23) J. L. García-Palacios, Adv. Chem. Phys. 112, 1 (2000).
- (24) B. D. Cullity and C. D. Graham, Introduction to Magnetic Materials (Wiley, New Jersey, 2009).
- (25) W. H. Press, S. A. Teukolsky, W. T. Vetterling, and B. P. Flannery, Numerical Recipes in C: the art of scientific computing 2nd ed. (Cambridge University Press, New York, 1992).
- (26) R. Kubo, M. Toda, and N. Hashitume, Statistical Physics II. Nonequilibrium Statistical Mechanics (Springer-Verlag, Berlin, 1992).
- (27) D. A. Garanin, Phys. Rev. E 54, 3250 (1996).
- (28) Y. P. Kalmykov, W. T. Coffey, U. Atxitia, O. Chubykalo-Fesenko, P.-M. Déjardin, and R. W. Chantrell, Phys. Rev. B 82, 024412 (2010).
- (29) W. T. Coffey, D. S. F. Crothers, Y. P. Kalmykov, and J. T. Waldron, Phys. Rev. B 51, 15947 (1995).
- (30) C. Aron, G. Biroli, and L. F. Cugliandolo, J. Stat. Mech. P11018 (2010).
- (31) K. Miyasaki and K. Seki, J. Chem Phys 108, 7052 (1998).
- (32) U. Atxitia, O. Chubykalo-Fesenko, R. W. Chantrell, U. Nowak, and A. Rebei, Phys. Rev. Lett. 102, 057203 (2009).