A Comparison of Three Time-stepping Methods for the LLG Equation in

TVE-F 17 005 maj

Examensarbete 15 hp Juni 2017

A Comparison of Three Time-stepping Methods for the LLG Equation in

Dynamic Micromagnetics

Anton Kroner

Tomas Berg

Simon Wredh

Teknisk- naturvetenskaplig fakultet UTH-enheten

Besöksadress:

Ångströmlaboratoriet Lägerhyddsvägen 1 Hus 4, Plan 0

Postadress:

Box 536 751 21 Uppsala

Telefon:

018 – 471 30 03

Telefax:

018 – 471 30 00

Hemsida:

http://www.teknat.uu.se/student

Abstract

A Comparison of Three Time-stepping Methods for the LLG Equation in Dynamic Micromagnetics

Anton Kroner & Tomas Berg & Simon Wredh

Micromagnetism is the study of magnetic materials on the microscopic length scale (of nano to micrometers), this scale does not take quantum mechanical effects into account, but is small enough to neglect certain macroscopic effects of magnetism in a material.

The Landau-Lifshitz-Gilbert (LLG) equation is used within micromagnetism to determine the time evolution of the magnetisation vector field in a ferromagnetic solid. It is a partial differential equation with high non linearity, which makes it very difficult to solve analytically. Thus numerical methods have been developed for approximating the solution using computers.

In this report we compare the performance of three different numerical methods for the LLG equation, the implicit midpoint method (IMP), the midpoint with

extrapolation method (MPE), and the Gauss-Seidel Projection method (GSPM).

It was found that all methods have convergence rates as expected; second order for IMP and MPE, and first order for GSPM. Energy conserving properties of the schemes were analysed and neither MPE or GSPM conserve energy. The computational time required for each method was determined to be very large for the IMP method in comparison to the other two. Suggestions for different areas of use for each method are provided.

Handledare: Doghonay Arjmand & Gunilla Kreiss

Popul¨ arvetenskaplig sammanfattning

Mikromagnetism ¨ ar ett omr˚ ade som studerar magnetism p˚ a mikroskopiska (dock ej kvant- mekaniska) l¨ angdskalor. Omr˚ adet ¨ ar av intresse bland annat inom h˚ arddiskindustrin, d¨ ar magnetiska material anv¨ ands f¨ or att lagra information.

Inom mikromagnetism finns Landau-Lifshitz-Gilberts (LLG) ekvation. L¨ osningar till denna beskriver hur magnetiseringen f¨ or¨ andras inom ett material ¨ over tid. LLG- ekvationen ¨ ar en partiell differentialekvation med variabler i rum och tid, vilket inneb¨ ar att l¨ osningen ¨ ar en funktion som beror p˚ a position och tidpunkt.

Ekvationen har en komplicerad struktur och att hitta exakta l¨ osningar ¨ ar inte l¨ att.

F¨ oljdaktligen anv¨ ander man numeriska metoder och datorer f¨ or att ber¨ akna ett n¨ armev¨ arde till l¨ osningen. Det finns ett antal olika metoder som ¨ ar anpassade f¨ or LLG-ekvationen och syftet med det h¨ ar projektet ¨ ar att j¨ amf¨ ora egenskaper f¨ or tre av dessa: implicita mittpunktsmetoden (IMP), den extrapolerade mittpunktsmetoden (MPE), samt Gauss- Siedel projektionsmetoden (GSPM).

Det uppt¨ acktes att alla metoder uppfyllde sina f¨ orv¨ antade noggranhetsordningar, andra ordningen f¨ or IMP och MPE, samt f¨ orsta ordningen f¨ or GSPM. Energibevar- ingsf¨ orm˚ agan unders¨ oktes och varken GSPM eller MPE klarade av att ber¨ akna energin korrekt. Ber¨ akningstiden som kr¨ avdes f¨ or IMP best¨ ammdes till att vara mycket h¨ ogre ¨ an f¨ or de andra tv˚ a metoderna. GSPM ¨ ar l¨ ampig f¨ or snabba och onoggranna simuleringar.

MPE b¨ or anv¨ andas f¨ or ber¨ akningar som inte kr¨ aver exakt energibevaring och som upp- fyller dess stabilitetsvillkor. IMP ger mest korrekta l¨ osningar oavsett diskretisering men

¨

ar relativt l˚ angsam.

Contents

1 Introduction 1

1.1 Background . . . . 1

1.2 Purpose . . . . 1

1.3 Simplifications . . . . 1

2 Theory 2 2.1 The LLG equation . . . . 2

2.2 Energy . . . . 3

2.3 Numerical properties . . . . 4

2.4 Numerical methods . . . . 4

2.4.1 Implicit midpoint method . . . . 4

2.4.2 Midpoint with explicit extrapolation . . . . 8

2.4.3 Gauss-Seidel projection . . . . 9

3 Method of comparison 11 3.1 Physical requirements . . . . 11

3.2 Numerical viability . . . . 11

4 Results 12 4.1 Magnitude conservation . . . . 13

4.2 Energy conservation . . . . 13

4.3 Stability analysis . . . . 15

4.4 Convergence studies . . . . 17

4.5 Computation time . . . . 18

5 Discussion 18 5.1 Conservative properties . . . . 18

5.2 Numerical performance . . . . 18

6 Conclusion 19

1 Introduction

1.1 Background

Magnetism in a material arises from the electron spin and orbital angular momentum in the atoms of which the material consists of. However, on larger scales, it is possible to analyse magnetism without taking these quantum mechanical effects into account.

The field of micromagnetics studies magnetism at the microscopic size scale of nano to micrometers, and has been of interest since the advent of modern computers. In the year 1960, the field had developed enough to be useful in the development of the magnetic hard drive [1].

A commonly used mathematical model in micromagnetics is the LLG equation. It was introduced by Landau and Lifshitz in 1935 and later modified by Thomas Gilbert.

The equation is a vectorial partial differential equation (PDE) and is used to describe the precessional motion of magnetisation in solids, as seen in figure 1. Analytical solutions can only be found in special cases and consequently numerical schemes have been developed to reliably solve the equation using computers. [1]

1.2 Purpose

The purpose of this project it to implement three different time stepping methods for the LLG equation and numerically compare their properties. A classical implicit midpoint method, known to be a good performer albeit slow is to be put up against two presumably faster methods: A mid point method with explicit extrapolation and a projection based scheme implemented with a Gauss-Seidel correction. Most importantly the magnitude of the solution vector has to be conserved for every time step. Furthermore energy conserva- tion is to be studied together with numerical stability, convergence rates and computation time. These properties are to be weighed against each other to draw conclusions about efficiency. The methods are to be implemented and tested with Matlab.

1.3 Simplifications

Due to the complexity of the equation we have limited ourself to a 1-dimensional setting,

whereas in general the problem is formulated over bounded domains in R

³

. The LL

equation, see equation (2.2), includes parameters, which are assumed equal to 1 in all

the simulations here. Moreover, due to the nature of the LL equation, the length of

the magnetisation m remains constant for all the time. Without loss of generality, the

magnitude is also assumed to be 1. As the interest of this project is the time stepping,

we have as much as possible kept to a constant space discretisation of 100 grid points on

the interval [0,100].

2 Theory

2.1 The LLG equation

The Landau-Lifshitz-Gilbert-equation is used to determine the time-evolution of the mag- netisation field in a ferromagnetic material (for a derivation see: [2]). To simplify equa- tions this paper will only consider its normalised form which can be written as

∂m

∂t = −



m × h

ef f

− am × ∂m

∂t



, (2.1)

where m(x, t) is the normalised magnetisation vector field of the ferromagnetic material, h

_{ef f}

is an effective magnetisation field and a is a damping constant. The equation can also be written on the equivalent [3] Landau-Lifshitz (LL) form

∂m

∂t = −βm × h

_{ef f}

− αm × (m × h

_{ef f}

), (2.2) which is explicit with respect to the time derivative of the magnetisation field. In this formulation β is a gyromagnetic constant and α is a damping constant. The LL version is suitable for explicit numerical schemes whereas the LLG variant can be beneficial for implicit schemes [4]. The effective magnetisation field is given by

h

_{ef f}

= 2A

µ

₀

M

_s²

∇

²

m + h

_ani

+ h

_ext

+ h

_d

, (2.3) where A is an exchange constant, µ

₀

is the permeability and M

_s

is the magnitude of the unnormalised magnetisation field. h

_ani

is the anisotropic magnetisation field, h

_ext

is the external magnetisation field and h

_d

is the demagnetisation field. To calculate h

_d

we would need to couple the LLG equation with one of Maxwell’s equations, however in this study we only want to study methods for the LLG equation so the demagnetisation field has been set to zero.

This study is limited to one space dimension which reduces the laplacian to a single second order derivative. If we have uniaxial anisotropy, along a unit vector, p, we can write the anisotropy field as [5]

h

ani

= hm, pip (2.4)

and as the external field is any known applied field, we can write the effective field as h

_{ef f}

= B ∂

²

m

∂x

²

+ hm, pip + h

_ext

. (2.5)

By taking the inner-product of equation (1.2) with m it is easy to see that |m(t, x)| =

|m(0, x)| for all x. Moreover, it is assumed that |m(0, x)| = 1. This is because the

magnitude of the magnetisation must be equal to the magnetisation saturation of the

material [6]. Figure 1 shows the precessional motion of a single magnetisation vector

under the effect of an effective h-field.

b a

H

eff

m y

z

x

Figure 1: The path of the magnetisation vector, m, as it precesses around the effective H field, where a and b are the first and second terms in the LLG equation, respectively [4].

2.2 Energy

Magnetic materials exhibit very complex behaviour due to the competition between dif- ferent energy terms. There are several contributions describing different phenomena that add up to a total free energy [4]. In this report we are considering three of these contri- butions.

The exchange energy is given by E

_exc

= 1

2 Z

J

_exc

|∇m|

²

, (2.6)

where J

exc

is a material constant.

The anisotropy energy is given by E

_ani

= − 1

2 Z

K

_ani

hm, pi

²

, (2.7)

where K

_ani

is a material constant and p is the direction of an easy axis in a uniaxial material.

The applied field energy is given by E

_ext

= −

Z

µ

₀

hh

_ext

, mi, (2.8)

where h

_ext

is an externally applied field. Important to note is that when the damping term vanishes the LL equation (2.2) conserves energy [4]. This conservative property is typically sought after when designing a numerical method for approximating the solution.

With N space steps the total continuous energy can be discretised as follows: [7]

E = − 1 2

N

X

k

J

exc

(m

k

· m

k−1

+ m

k

· m

k+1

) − 1 2

N

X

k

K

ani

(m

k

· p)

²

−

N

X

k

µ

0

h

ext

· m

k

(2.9)

2.3 Numerical properties

In numerical analysis there are certain core properties that play vital roles for the viability of a method. First of all, an explicit method calculates the solution using only current and previous time states. Comparatively a method of implicit nature expresses the solution at the upcoming time state in terms of the solution itself. Consequently an equation needs to be solved for every time step when using an implicit scheme.

The notion of stability, in general, is used to describe the sensitivity of the solution when the data in a given model problem are perturbed slightly. A necessary condition for the convergence of a numerical solution towards the exact solution is the stability of the numerical method, which may be fulfilled under some constraints on the discretisation parameters, as is the case for most explicit methods. Nevertheless, there are methods (implicit methods) which are unconditionally stable. Naturally, for a conditionally stable scheme to be usable at all the discretisation constraint has to be fulfilled.

To determine how accurate the solution of a method is one can analyse the order of convergence/accuracy. This is a number p that gives information about how the error  relates to a discretisation ∆t as shown in equation (2.10).

 ∝ ∆t

^p

(2.10)

Different methods use different algorithms that require certain amount of computing power. For large problems or small discretisation parameters the computational time can increase to unfeasible scales. Although dependent on hardware and code optimisations, a comparison can be achieved by simply timing the execution time.

2.4 Numerical methods

The nonlinearity of the LL and the LLG equations makes it very difficult to determine a solution analytically, this results in the use of numerical methods to reliably solve the equations. The part that gives rise to interesting problems is the time discretisation and there are different numerical methods that can be used to solve the equations, with their own strengths and weaknesses. This report will focus on analysing the properties of three specific methods.

2.4.1 Implicit midpoint method

Finite difference methods work by replacing the derivatives with finite differences that approximate the derivative at a discrete point,

dy(t

ⁿ⁺¹

)

dt ≈ y

ⁿ⁺¹

− y

ⁿ

∆t (2.11)

where n denotes the time step and y

ⁿ

is an approximation to the solution at y(t

ⁿ

). The differential equation then reduces to a series of simple algebraic equations that can be solved by ”stepping” as long as y

⁰

is known.

The popular Crank-Nicolson method uses an average of the Euler forward and Euler

backward methods to achieve second order accuracy. It may also be referred to as the

implicit midpoint (IMP) scheme and can be written on the form y

ⁿ⁺¹

− y

ⁿ

∆t = f (t

_n+¹

2

, y

ⁿ⁺¹²

), (2.12)

where

y

ⁿ⁺¹²

= 1

2 (y

ⁿ⁺¹

+ y

ⁿ

). (2.13)

Implicit means that unlike equation (2.11), it can not be directly solved because y

ⁿ⁺¹

depends on itself.

However, methods for solving non linear algebraic equations have been known for a long time, and the Newton-Rhapson method is commonly applied for solving problems of this type. Say we have some function g(x), the Newton-Rhapson method can is used to approximate a root, x ,that solves this equation by iteration given a starting ”guess”

x

₀

x

_i+1

= x

_i

− g(x

_i

)

g

⁰

(x

_i

) (2.14)

where g

⁰

(x) is the derivative of g(x). The stopping criterion for the algorithm is kx

_i+1

− x

_i

k <

tol for a given value of tol.

We can apply equation (2.14) to solve (2.12) by rearranging it and giving it a new name, G(y

ⁿ⁺¹

), we would then iterate the solution y

ⁿ⁺¹

in each time step, and use y

₀ⁿ⁺¹

= y

ⁿ

, the solution in the previous time step, as a starting guess.

Even though the LL equation is a vector PDE, the method works exactly the same.

We discretise in time by replacing the time derivative in equation (2.2) with a finite difference

∂m(t

ⁿ

)

∂t ≈ m

ⁿ⁺¹

− m

ⁿ

∆t (2.15)

and using the midpoint average, equation (2.13), on the right hand side we get m

ⁿ⁺¹

− m

ⁿ

∆t = −βm

ⁿ⁺¹²

× h

_{ef f}

− γm

ⁿ⁺¹²

× (m

ⁿ⁺¹²

× h

_{ef f}

). (2.16) Now define G(m

ⁿ⁺¹

) by rearranging equation (2.16) into

G(m

ⁿ⁺¹

) = m

ⁿ⁺¹

− m

ⁿ

∆t + βm

ⁿ⁺¹²

× h

_{ef f}

+ γm

ⁿ⁺¹²

× (m

ⁿ⁺¹²

× h

_{ef f}

). (2.17) Now we want to use Newton-Rhapsons method to solve for m

ⁿ⁺¹

, however, both G and m

ⁿ⁺¹

are vectors. Equation (2.14) can be expanded into solving systems of equations by replacing the scalar quantities by their respective vector equivalents. The derivative in equation (2.14) is replaced by the Jacobian matrix of G, J, which gives us

m

ⁿ⁺¹_i+1

= m

ⁿ⁺¹_i

− J

⁻¹

G(m

ⁿ⁺¹_i

) (2.18)

where as before we use m

ⁿ⁺¹₀

= m

ⁿ

. Notice that at this point, it is important to not

confuse the Newton-Rhapson iteration index, i, with the time stepping index, n. The

time stepping index increases only when the km

ⁿ⁺¹_i+1

− m

ⁿ⁺¹_i

k < tol condition is met, and

m

ⁿ⁺¹_i

”becomes” the solution for that time step after a series of iterations of over i. That

value is then used as an initial guess for calculating the solution at the next time step.

To calculate h

_{ef f}

we must now choose a discretisation in space. Even a one-dimensional domain has three components of the vector m at each point in space, so it is useful to represent the magnetization field by a N × 3 matrix, where N is the number of discrete space points, and the columns contain the components of the vector at each point in space.

m =



















m

_x0

m

_y0

m

_z0

m

_x1

m

_y1

m

_z1

m

_x2

m

_y2

m

_z2

.. . .. . .. . m

x_{N −1}

m

y_{N −1}

m

z_{N −1}

m

_xN

m

_yN

m

_zN



















Figure 2: Space discretisation matrix, where m

_x,y,z

denotes a component of m We replace the derivative in equation (2.5) by a finite difference approximation

∂

²

m(x

_k

)

∂x

²

≈ m

_k−1

− 2m

_k

+ m

_k+1

∆x

²

(2.19)

where k ∈ {0, 1, 2, . . . , N } denotes a point in space. By expanding equation (2.19) into the component form















∂²mx(xk)

∂x²

∂²my(xk)

∂x²

∂²mz(xk)

∂x²















≈ 1

∆x

²





m

_xk−1

− 2m

_xk

+ m

_xk+1

m

y_k−1

− 2m

y_k

+ m

y_k+1

m

_zk−1

− 2m

_zk

+ m

_zk+1



 (2.20)

we can see that this discrete derivative can be represented by a multiplication of each column of the space discretisation matrix, m, by a matrix

A = 1

∆x

²















−2 1

1 −2 1

1 −2 . ..

. .. ... 1 1 −2















N ×N

(2.21)

containing the coefficients of the finite difference approximation. This leads to

∂

²

m

∂x

²

≈ Am (2.22)

which is the second order space derivative. However, equation (2.22) is not complete

without providing boundary conditions. One type of boundary conditions are the periodic

boundary conditions, which means that we imagine the domain as a circle. To implement this, we modify A to contain two additional non zero elements in the corners

A = 1

∆x

²















−2 1 1

1 −2 1

1 −2 . ..

. .. ... 1

1 1 −2















N ×N

(2.23)

which will make it so that m

_{x,y,zN +1}

= m

_x,y,z0

and m

_x,y,z−1

= m

_x,y,zN

, completing equation (2.22).

We can now define a discrete version of the effective field, h

_{ef f}

, from equation (2.5).

In place of the derivative we use equation (2.22) and the other terms follow from their discrete counterparts and we can write

h

_{ef f}

= BAm + hm, pip + h

_ext

(2.24)

where p, h

_ext

and h

_{ef f}

have been extended to be matrices that contain their respective vector components at a point n at each row, analogous to m (see: figure 2.4.1). The inner product is taken at each individual row of m and p to return a scalar that scales each row of p.

Finally, we can combine both the time and space discretisations. We extend equation (2.17) by replacing each m with m and h

ef f

by h

_{ef f}

G(m

ⁿ⁺¹

) = m

ⁿ⁺¹

− m

ⁿ

∆t + βm

ⁿ⁺¹²

× h

_{ef f}

+ γm

ⁿ⁺¹²

× (m

ⁿ⁺¹²

× h

_{ef f}

) (2.25) where G has been extended to a matrix, as before. And

h

_{ef f}

= BAm

ⁿ⁺¹²

+ hm

ⁿ⁺¹²

, pip + h

_ext

(2.26) uses the midpoint average

m

ⁿ⁺¹²

= 1

2 (m

ⁿ⁺¹

+ m

ⁿ

) (2.27)

from equation (2.13) just like previously.

G is now three systems of N equations, but equation (2.18) still applies, except that for each iteration we have to compute a Jacobian for each column of G, multiply them together separately, and finally put them back into a matrix together. We write this as

m

ⁿ⁺¹_i+1

= m

ⁿ⁺¹_i

− J

⁻¹_x

G

x

(m

ⁿ⁺¹_i

) J

⁻¹_y

G

y

(m

ⁿ⁺¹_i

) J

⁻¹_z

G

z

(m

ⁿ⁺¹_i

) 

(2.28) where, just as before, we use m

ⁿ⁺¹₀

= m

ⁿ

as an initial value, J

_x,y,z

are the Jacobians of G

_x,y,z

, which are the columns of G, and i is the iteration index that runs until the tolerance km

ⁿ⁺¹_i+1

− m

ⁿ⁺¹_i

k < tol is fulfilled.

The Jacobians in each iteration are calculated by a central finite difference J

^k_x

= G

^k_x

(m

^k_x

+

^d₂

, m

^k_y

, m

^k_z

) − G

^k_x

(m

^k_x

−

^d₂

, m

^k_y

, m

^k_z

)

d

J

^k_y

= G

^k_y

(m

^k_x

, m

^k_y

+

^d₂

, m

^k_z

) − G

^k_y

(m

^k_x

, m

^k_y

−

^d₂

, m

^k_z

) d

J

^k_z

= G

^k_z

(m

^k_x

, m

^k_y

, m

^k_z

+

^d₂

) − G

^k_z

(m

^k_x

, m

^k_y

, m

^k_z

−

^d₂

) d

(2.29)

where k runs over every row of the spatial discretisation matrix, and d is a small number (≈ 10

⁻⁶

) used to compute the finite difference.

2.4.2 Midpoint with explicit extrapolation

The midpoint with explicit extrapolation (MPE) uses, just like Runge-Kutta methods for ordinary differential equations, a weighted average of several different time-steps when calculating the next time-step [4].

If one considers the LL equation (2.2), we see that it can be written on the following

form ∂m

∂t = −γm × H, (2.30)

where

H = h

ef f

+ αm × h

ef f

, (2.31)

where α is a damping constant. We can see in equation 2.30 that m ·

^∂m_∂t

= 0 and as a consequence that kmk = const [8].

To create a numerical method we introduce the following midpoint approximation

 dm dt



n+¹₂

= m

ⁿ⁺¹

− m

ⁿ

∆t + O(∆t

²

), (2.32)

m

ⁿ⁺¹²

= m

ⁿ⁺¹

+ m

ⁿ

2 + O(∆t

²

), (2.33)

where the subscript n indicates the time-step. The O(∆t

²

) is used to denote the order of accuracy in the time region, this being a order of accuracy of 2.

To derive a numerical method we also need to express H in the midpoint, we can do this by extrapolating H as follows

H

ⁿ⁺¹²

= h

ⁿ⁺¹²

+ O(∆t

²

), (2.34) where

h

ⁿ⁺¹²

= 3

2 H

ⁿ

+ 1

2 H

ⁿ⁻¹

. (2.35)

We can now substitute equations 2.32, 2.33 and 2.34 into the LL equation 2.30. This gives us the following expression

m

ⁿ⁺¹

− m

ⁿ

= γ∆t

2 (m

ⁿ⁺¹

+ m

ⁿ

) × h

ⁿ⁺¹²

, (2.36) which can be expressed as

m

ⁿ⁺¹

+ γ∆t

2 m

ⁿ⁺¹

× h

ⁿ⁺¹²

= m

ⁿ

− γ∆t

2 m

ⁿ

× h

ⁿ⁺¹²

. (2.37)

This gives rise to a system of equations that can be expressed as a product of matrices









1

^γ∆t₂

h

ⁿ⁺

1

z 2

−

^γ∆t₂

h

ⁿ⁺

1

y 2

−

^γ∆t₂

h

ⁿ⁺

1

z 2

1

^γ∆t₂

h

ⁿ⁺

1

x 2

γ∆t 2

h

ⁿ⁺

1

y 2

−

^γ∆t₂

h

ⁿ⁺

1

x 2

1











 m

ⁿ⁺¹_x

m

ⁿ⁺¹_y

m

ⁿ⁺¹_z





=









m

ⁿ_x

−

^γ∆t₂

(m

ⁿ_y

h

ⁿ⁺

1

z 2

− m

ⁿ_z

h

ⁿ⁺

1

y 2

) m

ⁿ_y

−

^γ∆t₂

(m

ⁿ_z

h

ⁿ⁺

1

x 2

− m

ⁿ_x

h

ⁿ⁺

1

z 2

) m

ⁿ_z

−

^γ∆t₂

(m

ⁿ_x

h

ⁿ⁺

1

y 2

− m

ⁿ_y

h

ⁿ⁺

1

x 2

)







 ,

where the subscripts x,y and z represents the components of the corresponding vectors.

Now that we have a numerical method for temporal discretisation we can take a look at the spatial discretisation. The spatial discretisation lies in the effective field h

_{ef f}

or more precisely the exchange field

_µ^2A

0M_s²

∇

²

m. To descretizise this we use a second order finite difference scheme with an order of accuracy of 2:

∇

²

m = ∂

²

m

∂x

²

= m

_k+1

− 2m

_k

+ m

_k−1

∆x

²

+ O(∆x

²

). (2.38)

This is then simply substituted into the effective field and solved for every k 2.4.3 Gauss-Seidel projection

The Gauss-Seidel projection method (GSPM) applies a Gauss-Seidel correction to a frac- tional step scheme for symplectic flow of harmonic maps and couples the solution with a projection method for heat flow of harmonic maps [9]. Considering the undamped LL equation with only the gyromagnetic term present

∂m

∂t = −m × ∆m, (2.39)

it can be linearised using a fractional step scheme [9]

m

^∗

− m

ⁿ

∆t = ∆

_h

m

^∗

, m

ⁿ⁺¹

= m

ⁿ

− m

ⁿ

× m

^∗

,

(2.40)

where ∆t is the time step size and ∆

_h

is the discrete laplacian. Combining the two equations one gets

m

ⁿ⁺¹

= m

ⁿ

− m

ⁿ

× (I − ∆t∆

_h

)

⁻¹

m

ⁿ

. (2.41) This scheme can be shown to be unstable when applied to a PDE such as the LLG equation [10]. Due to its vectorial product structure the stability of the scheme can be significantly improved by implementing an implicit Gauss-Seidel correction where new component values are reused immediately. In its component form the improved scheme looks as follows

g

ⁿ

= (I − ∆t∆

_h

)

⁻¹

m

ⁿ

(2.42)



 m

ⁿ⁺¹_x

m

ⁿ⁺¹_y

m

ⁿ⁺¹_z



 =





m

ⁿ_x

+ (g

ⁿ_y

m

ⁿ_z

− g

_zⁿ

m

ⁿ_y

) m

ⁿ_y

+ (g

_zⁿ

m

ⁿ⁺¹_x

− g

_xⁿ⁺¹

m

ⁿ_z

) m

ⁿ_z

+ (g

_xⁿ⁺¹

m

ⁿ⁺¹_y

− g

_yⁿ⁺¹

m

ⁿ⁺¹_x

)



 . (2.43)

When the gyromagnetic term is omitted and we have a case of infinite damping the LL equation becomes

∂m

∂t = −m × (m × ∆m), (2.44)

or equivalently

∂m

∂t = ∆m + ∇|m|

²

m. (2.45)

(2.45) describes heat flow for harmonic maps. An unconditionally stable two-step pro- jection scheme for the equation is shown in [11]. The term ∇|m|

²

m is considered a point wise constraint: hm,mi = 1. First solve

m

^∗

− m

ⁿ

∆t = ∆

h

m

^∗

(2.46)

for m

^∗

, then project onto the unit sphere

m

ⁿ⁺¹

= m

^∗

km

^∗

k . (2.47)

The two above methods for solving equations (2.39) and (2.45) can be coupled together to solve the full LL equation as the Gauss-Seidel projection method [10]. Considering a splitting of the effective field

h

_{ef f}

= ∆m + f, (2.48)

where

f = h

_ani

+ h

_ext

, (2.49)

the full equation can be written as

∂m

∂t = −βm × (∆m + f) − αm × m × (∆m + f). (2.50) GSPM can be applied to (2.50) in three steps:

Step 1: Implicit Gauss-Seidel for gyromagnetic term g

ⁿ

= (I − ∆t∆

_h

)

⁻¹

(m

ⁿ

− f)

g

^∗

= (I − ∆t∆

_h

)

⁻¹

(m

^∗

− f) (2.51)



 m

^∗_x

m

^∗_y

m

^∗_z



 =





m

^∗_x

+ (g

_yⁿ

m

ⁿ_z

− g

_zⁿ

m

ⁿ_y

) m

^∗_y

+ (g

ⁿ_z

m

^∗_x

− g

_x^∗

m

ⁿ_z

) m

^∗_z

+ (g

^∗_x

m

^∗_y

− g

^∗_y

m

^∗_x

)



 . (2.52)

Step 2: Projection method for heat flow of harmonic maps

f

^∗

= h

^∗_ani

+ h

_ext

(2.53)



 m

^∗∗_x

m

^∗∗_y

m

^∗∗_z



 =





m

^∗_x

+ α∆t(∆

_h

m

^∗_x

∗ +f

_x^∗

) m

^∗_y

+ α∆t(∆

_h

m

^∗_y

∗ +f

_y^∗

) m

^∗_z

+ α∆t(∆

_h

m

^∗_z

∗ +f

_z^∗

)



 . (2.54)

Step 3: Projection onto unit sphere



 m

ⁿ⁺¹_x

m

ⁿ⁺¹_y

m

ⁿ⁺¹_z



 = 1

km

^∗∗

k



 m

^∗∗_x

m

^∗∗_y

m

^∗∗_z



 . (2.55)

3 Method of comparison

To make a comparison several numerical properties and physical constraints have been investigated for each method. The specific initial value problem we have studied is



 

 

∂m

∂t

= −m × h

_{ef f}

− m × (m × h

_{ef f}

), m(0, t) = m(L + ∆x, t),

m(x, 0) = [cos(

^2πx_L

), sin(

^2πx_L

), 0],

(3.1)

with h

_{ef f}

=

^∂_∂x^2m2

+ h

_ani

+ h

_ext

, where h

_ani

= hm, pip , p =

^√1

3

[1, 1, 1] and h

_ext

= [0, 0, 1].

For error estimations an ”exact” solution calculated using IMP with ∆t = 10

⁻⁵

has been used.

3.1 Physical requirements

Essential from a physical viewpoint is that the magnitude is conserved for every time step, as explained in section 2.1. The restrictions of this project forces a constant magnitude of 1. To investigate whether the methods fulfilled this constraint correctly the maximum norm was calculated for every time step. Its deviation from 1 was then plotted to see the evolution of the magnitude error.

The gyromagnetic term of the LL equation is of energy conserving Hamiltonian struc- ture [12] . For a numerical method to correctly model the physics this implies that energy has to be conserved in the case of no damping. The energy conserving properties of the methods were checked by calculating and plotting the time evolution of the discrete en- ergy as seen in equation (2.9). Although a method does not compute energy correctly it may still result in a correct solution. This is a consequence of the interaction between the magnetic moments and the laplacian in the term for the exchange energy.

3.2 Numerical viability

To determine how viable a numerical method is, certain core properties should be anal- ysed, as detailed in section 2.3. To determine which discretisations satisfy correct solu- tions for a method its stability needs to be analysed. Due to the non linearity of the LL equation, performing an analytical stability analysis is difficult. Aside for IMP, which is implicit, we have performed numerical tests where the error is looked at for increasing time steps with different space steps. If there is a major error increase for a certain dis- cretisation it can be concluded that stability has been lost. Stability results were double checked by comparing energy values.

How accurate a method is can be concluded from its order of accuracy. This prop- erty determines how fast the error decreases with smaller step sizes and is consequently important for the choice of discretisation. The orders of accuracies were calculated for each scheme by looking at the errors for different time steps and inserting these values in equation (2.10).

The last property analysed was the computation time. It is obviously preferable for

a method to be as fast as possible and for large simulations the feasibility might even be

constrained by long execution times. The computation time was investigated by timing

the three methods for different amounts of time steps. Linear regression was performed on the results to get an approximate function for the computation time.

4 Results

General results are compiled in table 1 and a time evolution of the solution can be seen in figure 3.

Table 1: General results for the three schemes

Scheme Mag. conserv. Energy conserv. Stability Conv. rate Comp. time

IMP Yes Yes Uncond. 2 Long

MPE Yes No Cond. 2 Short

GSPM Yes No Uncond. 1 Short

t = 0

0 50 100

x -1

0 1

y -1

0 1

z

-1 1 0

100

z

t = 0.5

y 0

x 1

-1 0 50

-1 1 0

100

z

t = 1

y 0

x 1

50 -1 0

-1 1 0

100

z

t = 2

y 0

x 1

50 -1 0

-1 1 0

100

z

t = 3

y 0

x 1

-1 0 50

-1 1 0

100

z

t = 5

y 0

x 1

-1 0 50

Figure 3: Solution at different times t. Computed using GSPM with ∆t = 10

⁻⁴

and 100

grid points on [0,100]. The solution converges towards the effective h-field.

4.1 Magnitude conservation

The magnitude conservation was confirmed by finding the largest magnitude error |m|−1 for every time step. Errors for ∆t = 10

⁻³

solutions are plotted in figure 4 and as can be seen the fluctuations are very small and negligible. The errors for GSPM are on the scale of machine epsilon.

0 500 1000 1500 2000 2500 3000 3500 4000 4500 5000 0

2 4

6 10

^-14

Implicit midpoint

0 500 1000 1500 2000 2500 3000 3500 4000 4500 5000 0

1 2

|m|-1

10

^-14

Midpoint extrapolation

0 500 1000 1500 2000 2500 3000 3500 4000 4500 5000 Time steps

0 1 2

3 10

^-16

Gauss-Seidel projection

Figure 4: Magnitude fluctuations around 1 for the three methods with ∆t = 10

⁻³

and 100 grid points on [0,100]

4.2 Energy conservation

With zero damping the equation reduces to

∂m

∂t = −m × h

_{ef f}

, (4.1)

which is of Hamiltonian structure and conservative with respect to energy. To observe

how well the methods conserve this property their discrete energies were calculated using

equation (2.9). The energy evolution over time for the three methods are shown in figure 5. It is clear from the plots that implicit midpoint conserves energy whereas the other two schemes do not. However, the projection method diverges a lot faster than midpoint with extrapolation which is very clear from figure 6 that shows the energies in the same plot.

0 10 20 30 40 50 60 70 80 90 100

-116.5 -116.4

-116.3 Implicit midpoint

0 10 20 30 40 50 60 70 80 90 100

-116.5 -116.4 -116.3

Energy

Midpoint extrapolation

0 10 20 30 40 50 60 70 80 90 100

Time -180

-160 -140 -120

-100 Gauss-Seidel projection

Figure 5: Energy evolution without damping for the three methods with ∆t = 10

⁻²

,

T = 100 and 100 grid points on [0,100].

0 20 40 60 80 100 Time

-119 -118.5 -118 -117.5 -117 -116.5 -116

Energy

Energy evolution without damping

Implicit midpoint Midpoint extrapolation Gauss-Seidel projection

Figure 6: Energy evolution without damping for the three methods on the same energy scale with ∆t = 10

⁻²

, T = 100 and 100 grid points on [0,100].

4.3 Stability analysis

Due to the complexity of the LL equation performing analytical stability analyses is rather difficult. Instead one can analyse the methods numerically by observing how the error behaves at different time steps. If the error suddenly spikes and starts to behave erratically the stability criteria is no longer fulfilled and the method has become unstable.

The IMP method is of implicit nature and therefore unconditionally stable. According to [9] the GSP method is also unconditionally stable. Test runs with large time steps have not been able to disprove this thesis and it is thus assumed true.

MPE is an explicit scheme and should lose stability for large step sizes. Figure 7 shows three plots for different space steps around time steps where there is a marginal increase in maximum error. The error was approximated by comparing the last time step with solutions of finer discretisations. From this information an approximation of the stability relation was calculated. The plots in figure 7 implies a stability relation that follows the trend

_∆x^∆t2

≤ C, where C is some constant. From table 2 it is possible to approximate the value of the constant C to: 0.1539 ± 3.22%.

This result can also be confirmed by analysing the energy evolution with damping

around time steps close to the stability limit. As shown in figure 8 the energy starts

increasing after supposed convergence for time steps larger than the stability condition.

Table 2: ∆t for loss of stability for different ∆x

∆t ∆x

0.143 1 0.087 0.75 0.041 0.5

0.02 0 0.04 0.06 0.08 0.1 0.12 0.14 0.16 0.18 0.2 1

2

ERROR

Stability dx=0.75

0.02 0.04 0.06 0.08 0.1 0.12 0.14 0.16 0.18 0.2 Time Step

0 2

4 Stability dx=0.5

0.02 0.04 0.06 0.08 0.1 0.12 0.14 0.16 0.18 0.2 0

0.5

1 Stability dx=1

Figure 7: MPE: Error as a function of ∆t for three different space steps.

0 2 4 6 8 10 Time

-240 -220 -200 -180 -160 -140 -120 -100

Energy

Energy evolution with damping for different t

t = 0.143 t = 0.149 t = 0.159 t = 0.189

Figure 8: MPE: Energy over time with damping for different time steps with 100 grid points on [0,100].

4.4 Convergence studies

Convergence rates were tested by using an IMP solution of (3.1) with ∆t = 10

⁻⁵

and 100 grid points on [0,100] as ”exact”. Absolute errors for decreasing time steps of each of the methods were then computed. The errors were calculated as the maximum error norms of all space steps at the end time T = 5. These results confirm an order of accuracy of 2 for IMP and MPE and 1 for GSPM.

Table 3: Errors at time T = 5 for different ∆t

∆t IMP MPE GSPM

0.1 2.9842E-4 1.0808E-3 8.7727E-3

0.1/2 7.4119E-5 2.6964E-4 4.0920E-3

0.1/4 1.8570E-5 6.7320E-5 1.9713E-3

0.1/8 4.7222E-6 1.6996E-5 9.6701E-4

0.1/16 1.2874E-6 4.4605E-6 4.7894E-4

Order of accuracy 1.9862 1.9829 1.0471

4.5 Computation time

A major difference between the IMP method and the other methods is the computational time required to complete a time step. Computational time was measured for each of the methods to find the dependency on the amount of time steps. In table 4 we can see that the GSPM is the fastest, with a computational time of about 2 seconds per 1000 time steps. The MPE method is about half as fast, but the IMP method has a time cost that is more than two orders of magnitude greater than the GSPM.

Table 4: Computation time, in seconds, for each method using K time steps with ∆x = 1.

K IMP MPE GSPM

100 18.43 0.59 0.29

400 60.96 1.54 0.75

800 105.17 2.93 1.42

1000 132.17 3.54 1.76

4000 423.34 13.67 7.36

8000 834.59 27.33 14.54

10000 1017.66 34.09 18.38

T

_comp

≈ 0.1K + 22 0.0034K + 0.2 0.0018K

5 Discussion

5.1 Conservative properties

We know from the theory that kmk is supposed to be constant over time. Looking at figure 4 that shows how the magnitude error evolves we can see that all three methods conserve this property well. The errors are very small and should be negligible in almost all cases. The IMP plot has a rather odd increase but the error is still very small and in a run with more time steps the error stayed at the same order of magnitude, similar to the MPE case. For GSPM m is normalised in every time step and consequently the error is constantly on the scale of machine epsilon.

Figure 6 shows the energy over time and we can see that the IMP method conserves energy very well. This agrees with [13] which shows analytically that the method should conserve energy. Furthermore it is clear that both GSPM and MPE fails to keep the energy constant. Although, whereas the projection method quickly diverges the midpoint scheme does fairly well with a relatively low margin of error. In situations where the energy quantity is of importance the GSPM is consequently not ideal to use but MPE might be suitable in some cases.

5.2 Numerical performance

IMP is an implicit scheme and therefore unconditionally stable. The lack of stability

constraint is one of the main benefits of this method. Combined with an order of con-

vergence of two and its energy conserving properties the IMP method is a popular LLG

solver. However, due to the implicit nature its downsides are implementation difficulties

as well as long computation times. (Although not considered in this project, it should be possible to reduce the computational time of the IMP method by disregarding a few elements when computing the approximation to the Jacobian matrix [13].)

Comparatively MPE is explicit and as figure 7 shows stability is lost for large time steps in relation to the space step. The constraint for the studied problem was approx- imated to be ∆t 6 0.154∆x

²

which allows for some discretisation variety but not as freely as with IMP. As the order of convergence also is two what differs the most from IMP is computation time where MPE performs a lot better. In cases where the stability constraint is not limiting and energy does not have to be very accurate MPE seems to be the preferred mid point scheme.

From table 4 we can see that GSPM is the fastest scheme, although similar to MPE.

Seemingly unconditionally stable as well, it is suited for quick computations. However, it is unable to conserve energy and only has an order of convergence of one. It is thus debatable whether GSPM can hold up to the two mid point methods when more accurate solutions are required.

6 Conclusion

We implemented three different methods for solving the LLG equation. Experiments to confirm convergences rates were performed and the results match expected rates from cited sources. The energy conservative properties of the solutions to the undamped LLG equation were investigated and it was found that the GSPM and the MPE method failed to compute this correctly. The stability constraint for the MPE method was approx- imated. The computational time for the IMP method was found to be two orders of magnitude greater than the other two methods.

We can therefore conclude that the IMP should be used in cases where energy com-

putation, or the use of large time steps is of importance. The GSPM should be used for

fast calculations where the energy is of negligible importance. And the MPE could be

used in cases where some energy inaccuracy is allowed and necessary discretisations fulfil

the stability criterion as a faster replacement for the IMP.

References

[1] D. Wei, Micromagnetics and recording materials. Springer Science & Business Media, 2012.

[2] B. Gurney, M. Carey, C. Tsang, M. Williams, S. Parkin, R. Fontana Jr, E. Gro- chowski, M. Pinarbasi, T. Lin, and D. Mauri, “Ultrathin magnetic structures iv, edited by b. heinrich and jac bland,” 2005.

[3] G. Bertotti, Hysteresis in magnetism: for physicists, materials scientists, and engi- neers. Academic press, 1998.

[4] I. Cimr´ ak, “A survey on the numerics and computations for the landau-lifshitz equation of micromagnetism,” Archives of Computational Methods in Engineering, vol. 15, no. 3, pp. 1–37, 2007.

[5] L. Baˇ nas, “Numerical methods for the landau-lifshitz-gilbert equation,” Numerical Analysis and Its Applications, pp. 158–165, 2005.

[6] A. Aharoni, Introduction to the Theory of Ferromagnetism, vol. 109. Clarendon Press, 2000.

[7] D. Arjmand, M. Poluektov, and G. Kreiss, “Atomistic-continuum multiscale modelling of magnetisation dynamics at non-zero temperature,” arXiv preprint arXiv:1702.05173, 2017.

[8] C. Serpico, I. Mayergoyz, and G. Bertotti, “Numerical technique for integration of the landau–lifshitz equation,” Journal of Applied Physics, vol. 89, no. 11, pp. 6991–

6993, 2001.

[9] X.-P. Wang, C. J. Garcıa-Cervera, and E. Weinan, “A gauss–seidel projec- tion method for micromagnetics simulations,” Journal of Computational Physics, vol. 171, no. 1, pp. 357–372, 2001.

[10] C. J. Garcia-Cervera, “Numerical micromagnetics: A review,” 2007.

[11] E. Weinan and X.-P. Wang, “Numerical methods for the landau–lifshitz equation,”

SIAM journal on numerical analysis, 2006.

[12] E. Kim and K. Lipnikov, “The mimetic finite difference method for the landau–

lifshitz equation,” Journal of Computational Physics, vol. 328, pp. 109–130, 2017.

[13] M. d’Aquino, C. Serpico, and G. Miano, “Geometrical integration of landau–lifshitz–

gilbert equation based on the mid-point rule,” Journal of Computational Physics,

vol. 209, no. 2, pp. 730–753, 2005.