Chapter 5 Numerical Time Integration
of the LLG Equation
The temporal evolution of magnetization in micromagnetic systems is governed by the LLG equation introduced in Chapter 3. This nonlinear partial differential equation
couples precessional dynamics with phenomenological damping. After spatial discretization by finite-element or finite-difference methods (Chapter 4), the LLG equation
reduces to a large system of coupled ordinary differential equations (ODEs). The efficient and accurate numerical integration of this system is essential for practical micromagnetic simulations, yet poses several fundamental
challenges that distinguish it from generic ODE problems.
This chapter develops the theoretical framework for time integrating the spatially discretized LLG equation and presents a systematic comparison of the integration schemes implemented in this work. Section 5.1 identifies the key numerical challenges as stiffness, the unit-sphere constraint, and geometric structure, that guide the choice of time integrator. Section 5.2 reviews the classical taxonomy of explicit, implicit, and implicit-explicit (IMEX) methods. Section 5.3
introduces the tangent-plane scheme, a reformulation of the LLG equation that transforms the constrained nonlinear problem into a sequence of linear saddle-point systems. Section 5.4 discusses adaptive multistep methods provided by the SUNDIALS library. Finally, Section 5.5
presents verification studies using the \(\mu \)MAG Standard Problem 4 and compares the accuracy of the implemented integrators.
5.1 Challenges in LLG Time Integration
Numerical integration of the LLG equation presents three interrelated challenges that must be addressed by any practical time-stepping scheme: stiffness arising from disparate time scales, preservation of the unit-sphere constraint,
and retention of the underlying geometric structure.
5.1.1 Stiffness and Disparate Time Scales
A differential equation is called stiff when its solution contains components evolving on vastly different time scales, and when explicit methods require impractically small time steps to maintain numerical stability rather
than accuracy [156]. The LLG equation is inherently stiff due to the exchange interaction.
The exchange field contribution to the effective field takes the form:
where \(A\) is the exchange stiffness constant. Upon spatial discretization with mesh size \(h\), the Laplacian operator introduces eigenvalues scaling as \(h^{-2}\), leading to an effective exchange field magnitude:
For typical ferromagnetic materials (\(A \approx 10^{-11}\) J/m, \(M_s \approx 10^6\) A/m) and fine meshes (\(h \approx 1\) nm), this yields exchange fields of order \(10^4\)–\(10^5\) A/m,
corresponding to precessional frequencies:
The associated time scale \(\tau _\text {exch} = 2\pi /\omega _\text {exch} \sim 10^{-13}\)–\(10^{-12}\) \(\si {\second }\) is orders of magnitude shorter than the physical time scales of interest in device
applications, such as magnetization switching times of order \(10^{-9}\) \(\si {\second }\).
For explicit time-stepping schemes, the Courant–Friedrichs–Lewy (CFL) stability condition requires:
where \(C\) is an order-unity constant depending on the specific method. This constraint becomes increasingly restrictive as the mesh is refined, rendering explicit schemes computationally impractical for high-resolution
simulations [157, 158, 159].
While implicit time-stepping schemes bypass this strict CFL stability limit, allowing significantly larger time steps, stiffness re-emerges as a bottleneck in linear algebra. The large eigenvalues of the exchange operator result in a
system matrix \(\mathbf {A}\) with a larger condition number \(\kappa _2(\mathbf {A})\), defined in terms of the spectral matrix norm \(\|\cdot \|_2\) as:
Figure 5.1: Hierarchy of time scales in micromagnetic simulations. The exchange interaction imposes the CFL stability limit (\(\tau _\mathrm {exch} \sim 10^{-12}\) s for \(h = 2\) nm mesh). Demagne-
tization dynamics occur at \(\tau _\mathrm {demag} \sim 10^{-11}\) s. Precession due to external fields operates at \(\tau _\mathrm {prec} \sim 10^{-9}\) s. Damping and switching occur at \(10^{-9}\) –
\(10^{-8}\) s, far above the explicit stability limit. Material: Permalloy (\(A = 1.3 \times 10^{-11}\) J/m, \(M_s = 8 \times 10^5\) A/m, \(\alpha = 0.02\)).
Since the maximum eigenvalue of the discretized exchange operator scales with the inverse square of the mesh size (\(\lambda _\text {max} \propto h^{-2}\)), refining the mesh drastically increases \(\kappa _2(\mathbf
{A})\). This ill-conditioning causes standard iterative solvers (such as Generalized Minimal Residual (GMRES)) to stagnate, necessitating the use of effective preconditioners.
Figure 5.1 illustrates the separation of time scales in a typical micromagnetic simulation. The exchange-driven dynamics impose a stability limit that is far
below the time scales associated with damping, demagnetization, and variations in the external field.
5.1.2 The Unit-Sphere Constraint
A fundamental property of the LLG equation is that the magnetization magnitude is conserved:
This constraint follows directly from the structure of the LLG equation. Taking the inner product of Equation (3.9) with \(\mathbf
{m}\) and using the orthogonality of cross products yields:
This approach is simple but introduces \(\mathcal {O}(\Delta t^p)\) perturbations, where \(p\) is the order of the underlying scheme.
(ii)Constraint-preserving schemes: Certain schemes, such as the implicit midpoint rule applied to the LLG equation in Gilbert form, preserve the constraint exactly at the discrete
level [94].
(iii)Tangent-plane formulation: The equation is reformulated to solve for the velocity \(\mathbf {v} = \partial _t \mathbf {m}\) in the tangent space \(T_\mathbf {m} S^2
= \{\mathbf {v} : \mathbf {v} \cdot \mathbf {m} = 0\}\), enforcing orthogonality either exactly via Lagrange multipliers or approximately via projection [135, 160].
The tangent-plane approach, developed in detail in Section 5.3, is employed throughout this work.
5.1.3 Structural Properties
The LLG equation possesses a rich geometric structure that certain numerical schemes can exploit or preserve. In the absence of damping (\(\alpha = 0\)), the equation reduces to:
which is a Hamiltonian system with the total micromagnetic energy as the Hamiltonian. The magnetization precesses on surfaces of constant energy, and the dynamics preserve both the energy and the phase-space volume
(Liouville’s theorem).
When damping is present (\(\alpha > 0\)), the system becomes dissipative. As derived in Chapter 3 (cf. Equation (3.15)), the energy decreases monotonically [94]:
A necessary property of numerical schemes is to preserve this dissipative character at the discrete level, ensuring that artificial energy growth does not occur.
