4.1 Introduction

The purpose of this chapter is to address the general problem of the numerical integration of Landau-Lifshitz-Gilbert (LLG) equation. In fact, we have observed more than once, that due to the nonlinear nature of this equation, analytical solutions can be derived in very few particular cases, or by using linearization techniques (see chapter 2 and references therein). Consequently, the only general (and mostly used) method to study magnetization dynamics is to solve LLG equation by suitable numerical methods. In this respect, as mentioned in chapter 3, the most common procedure is to use a semi-discretization approach. First, the equation is only discretized in space by using finite difference or finite elements methods [52]. This leads to a discretized version of the micromagnetic free energy and a corresponding system of ordinary differential equations (ODEs). Second, this system of ODEs is numerically integrated by using appropriate time-stepping techniques. It is interesting to underline that, while the spatial discretization is generally carried out trying to preserve the main properties of the free energy functional $G(\textbf{M};\mathbf{H}_a)$ introduced in chapter 1, little attention is generally paid to the preservation, after the time discretization, of the peculiar structure of LLG temporal evolution. This is probably due to the fact that, in the past, the main emphasis was on static micromagnetics and on reproducing accurate approximation of the free energy landscape associated to a magnetic system subject to quasi-static external fields. This goal has been generally achieved by using accurate spatial discretization of the free energy $G(\textbf{M};\mathbf{H}_a)$. On the other hand, when dynamic magnetization processes have to be investigated, the issue of using appropriate numerical time integration technique becomes rather crucial. Nevertheless, this problem seems to have been substantially overlooked, and most workers in LLG numerical simulation use `off-the-shelf' algorithms such as Euler, linear multi-step methods (e.g. Adams-Bashforth, Adams-Moulton, Crank-Nicholson, Backward Differentiation Formulas (BDF)) or Runke-Kutta methods[79,80]. We must underline here that these standard methods do not preserve structural properties of LLG time evolution. This equation has indeed peculiar dynamic properties, mentioned in section 1.3.5, which it is convenient to recall below.

a)

First, the magnetization has constant magnitude in time at each spatial location, as indicated by Eq. (1.97). Equation (1.97) is a fundamental constraint on the LLG time evolution that should be respected in the time discretized version of LLG equation. Since, usual time stepping methods do not preserve this property, most researchers follows the naive approach of renormalizing the magnetization vector field at each time step or after a prescribed tolerance has been exceeded. This naive approach is actually a nonlinear numerical modification of the time evolution which in principle can have also relatively strong effect on the subsequent computation of magnetostatic field [86] and for this reason is not recommended, especially when long time regime have to be studied.

b)

Second, for constant external field the LLG evolution has a Lyapunov structure [82], namely the free energy functional is a decreasing function of time along the trajectories of LLG equation, according to Eq. (1.102). This property is fundamental because it guarantees that the system tends toward stable equilibrium points, which are minima of free energy. Usual time-stepping techniques preserve this property only for sufficiently small time-step. Indeed, when the time-step is too large, instability phenomena can produce transient or even steady increase of energy. The stability constraint on the time step is usually rather severe and this generally leads to unnecessary long computational time.

c)

Third, the LLG equation is obtained by adding a phenomenological damping term to an otherwise hamiltonian (conservative) equation and therefore one should expect that in the limit of $\alpha \rightarrow 0$, the numerical integration should preserve energy and, if possible, the hamiltonian structure. This is not only a mathematical requirement. In fact, in most experimental situations LLG evolution is not strongly dissipative and the damping effects can be considered as a perturbation of the conservative motion (see chapter 2, sections 2.4.3, 2.5, 2.6). In this respect, it is quite reasonable from the physical point of view, that the numerical integration scheme is able to reproduce accurately the conservative motion. This is definitely the most challenging part in the numerical simulation since the conservative precession is generally much faster than the slow motion associated to dissipative processes. As it is well-known in hamiltonian dynamics studies, most standard numerical schemes do not preserve energy and/or hamiltonian structure, and particular care must by devoted to develop appropriate time stepping technique.

As matter of fact, it is generally very difficult to obtain the preservation of the above properties in the time discretization by using explicit methods (e.g. Euler, Adams-Bashforth). Implicit methods, on the other hand, have good performances in terms of stability, but do not preserve the amplitude of magnetization or the energy in the limit $\alpha \rightarrow 0$. However, the use of implicit methods generally makes it necessary to solve large system of coupled nonlinear equations at each time step which may lead to unacceptable computational cost. In this respect, most researchers generally try to avoid implicit methods by using appropriate semi-implicit techniques [81]. This has of course the drawback that accurate numerical time integration require stability upper bound for the time step. This in turn can be quite problematic since LLG dynamics, in many relevant cases, may exhibit dynamic processes with very different time scales. In fact, the issue of developing time integrators for LLG equation that preserve relevant properties of the equation under discretization, has received lately some attention [25,84,85,86]. The general point of view presented in these recent works is to use suitable geometric integrators[87] which are techniques designed to preserve geometrical properties of dynamics, namely symmetry, conservation of quantities, hamiltonian structure etc. In particular, the possibility of developing integrators for LLG equation based on Lie-group methods and Cayley transform have been investigated in Refs. [85,86]. These methods preserve the magnetization amplitude, but they do not generally preserve the LLG Lyapunov structure and the energy in the limit of zero damping. The basic idea is to take into account the conservation of magnetization magnitude by an appropriate change of variable (lift of the problem in the Lie-algebra associated to the Lie-group of rotations). The problem is then solved with usual RK time-stepping algorithms. These methods are conditionally stable and the stability requirements are certainly affected by the choice of the new set of variables. This could lead to an increase of the temporal stiffness and, consequently, to an increase of the computational cost. In this chapter, we will apply the (implicit) mid-point rule to the time integration of LLG equation. We shall demonstrate that the use of mid-point rule leads to a numerical time stepping that preserves the fundamental properties of LLG dynamics. This algorithm has been known for a long time and it has been applied extensively in the area of hamiltonian dynamics for its interesting preservation properties [83]. However, in its pure form, it has never been applied directly to the full 3D micromagnetics dynamical problem. A partial use of mid-point rule has been proposed in Ref. [25]. In this work, the mid point rule has been applied along with an appropriate explicit extrapolation formula (second order Adams-Bashforth) for the effective field. This method has the property of preserving magnetization magnitude and, due to the explicit extrapolation formula, does not require the inversion of a large system of coupled nonlinear equations (but just three by three linear system of equations at each location in space). However, the method does not in general preserve the Lyapunov structure of LLG equation neither the energy for zero damping. In addition, the semi-implicit nature of the scheme imposes stability restrictions to the time step. Here we apply the implicit mid-point rule directly to the LLG equation in its pure form. With the use of the mid-point rule we can overcome the drawbacks of the standard methods. The methods is unconditionally stable, preserves exactly, independently from the time step, magnetization magnitude and, in the case of zero damping, free energy of the system. In addition, mid-point rule preserves unconditionally Lyapunov structure of LLG dynamics for constant applied field, namely in the discrete dynamics, the free energy is always decreasing regardless of the time step. The price we have to pay is that now we have to solve a large (generally full) system of nonlinear algebraic equations. As we will discuss in the following, this problem has been dealt with by using quasi-Newton algorithm which allows one to deal with sparse banded matrix inversions only.