Time integration schemes for ODEs are commonly classified as explicit or implicit, and as one-step or multistep methods. This section reviews the classes most relevant for integrating the LLG equation.
5.2.1 Explicit Runge–Kutta Methods
The LLG equation can be written as a general initial value problem \(\dot {\mathbf {m}} = f(t, \mathbf {m})\), where \(f\) denotes the right-hand side (RHS) of Equation (3.12) including all effective field contributions and, where applicable, spin-transfer torque terms. Explicit Runge–Kutta methods advance this system from \(t_n\) to
\(t_{n+1} = t_n + \Delta t\) via:
where \(A = (a_{ij}) \in \mathbb {R}^{s \times s}\) is the Runge–Kutta matrix, \(\mathbf {b} = (b_1, \ldots , b_s)^{\mathsf {T}}\) the quadrature weights, and \(\mathbf {c} = (c_1, \ldots , c_s)^{\mathsf
{T}}\) the stage nodes at which the right-hand side is evaluated. For an explicit method, \(A\) is strictly lower triangular.
Classic examples include the forward Euler method (\(s = 1\), first-order), Heun’s method (\(s = 2\), second-order), and the classical fourth-order Runge–Kutta (RK4) method (\(s = 4\)). For stiff problems, explicit methods
are subject to severe stability restrictions. The stability region of the forward Euler method, for instance, is the unit disk centered at \((-1, 0)\) in the complex plane. Eigenvalues of the linearized LLG operator corresponding to
exchange modes lie far to the left of this region, necessitating time steps satisfying Equation (5.4) (see Figure 5.2). The classical RK4 method has an extended stability region, extending to approximately \(|\lambda \Delta t| < 2.8\) along the negative real axis.
For the exchange-dominated eigenvalues \(|\lambda _\text {max}| \sim 10^{13}\,\si {\radian \per \second }\) typical of fine micromagnetic meshes, stability requires \(\Delta t \lesssim 10^{-13}\,\si {\second
}\). While this constraint necessitates many time steps (e.g., 20000 steps for a \(2\,\si {\nano \second }\) simulation), RK4 remains a viable option when the computational cost per step is low and fourth-order accuracy in
constraint preservation is desired.
Figure 5.2: Stability regions in the complex plane \(z = \lambda \Delta t\) with physical eigenvalues. The Physics of Interest (Precession, Damping) lies near the origin and could be stable for explicit methods.
The exchange eigenvalue position depends on the time step: for \(\Delta t = 10^{-13}\,\si {\second }\), \(\lambda _\mathrm {exch} \Delta t \approx -1\) lies within the RK4 stability region (black dot). Increasing
\(\Delta t\) moves it left (gray dot), eventually exiting explicit stability regions. Explicit solvers (Red/Blue) are forced to take tiny time steps (\(\Delta t \to 0\)) just to keep the black dot inside their stability boundary. The
Implicit solver (Green) remains stable for any \(\Delta t\), allowing us to step over the stiff exchange modes. IMEX methods combine both approaches: exchange modes are treated implicitly (Green), while precession and damping
remain explicit (Red/Blue).
5.2.2 Implicit Methods and A-Stability
Implicit methods evaluate the RHS at the unknown solution value \(\mathbf {m}_{n+1}\), thereby enlarging the stability region. The simplest example is the backward Euler method:
This method is A-stable: its stability region includes the entire left half of the complex plane, meaning that all stable modes of the continuous problem remain stable regardless of the time-step size [156].
It is also A-stable and additionally preserves the unit-sphere constraint for the LLG equation when combined with appropriate spatial discretization [90].
The computational cost of implicit methods arises from solving the nonlinear system that occurs at each time step. For the LLG equation, Newton iteration requires the Jacobian \(\partial f / \partial \mathbf {m}\), which
is dense due to the demagnetization field. However, the tangent-plane reformulation (Section 5.3) transforms the problem into a linear system for each
time step, dramatically reducing the computational overhead.
5.2.3 Multistep Methods: BDF Schemes
Backward differentiation formulas (BDFs) are implicit linear multistep methods that approximate the time derivative using values from multiple previous time steps. The \(k\)-step BDF method takes the form:
This achieves second-order accuracy. BDF methods up to order 6 exist, but only BDF1 and BDF2 are A-stable. The higher-order methods BDF3–BDF5 are \(A(\alpha )\)-stable, meaning their stability region contains a wedge
\(|\mathrm {arg}(-\lambda )| < \alpha \) for some angle \(\alpha < \pi /2\) rather than the full left half-plane [156]. This weaker condition, also referred to as stiff stability, is sufficient for problems where all
stiff eigenvalues lie near the negative real axis, as is the case for exchange-dominated micromagnetic systems. BDF6 is not even \(A(\alpha )\)-stable and is therefore not used in practice.
The C-language Variable-coefficient Ordinary Differential Equation solver (CVODE) solver from the SUNDIALS library [161] implements variable-order, variable-step BDF methods with automatic step-size and
order selection based on local error estimates. This adaptive capability is particularly valuable for micromagnetic simulations, where the dynamics may alternate between rapidly varying transients and slowly evolving quasi-static
phases.
5.2.4 Implicit-Explicit Method
When the RHS naturally splits into stiff and non-stiff components, it can be written as:
IMEX methods treat \(f_I\) implicitly for stability while evaluating \(f_E\) explicitly for efficiency [162]. For the LLG equation, we assign the stiff exchange-damping contribution to the implicit partition:
\(\seteqnumber{0}{5.}{18}\)
\begin{align}
f_I &= -\gamma _L \alpha \, \mathbf {m} \times (\mathbf {m} \times \mathbf {H}_\text {exch}), \label {eq::TI_fI}\\ f_E &= -\gamma _L \mathbf {m} \times \mathbf {H}_\text {eff} - \gamma
_L \alpha \, \mathbf {m} \times (\mathbf {m} \times \mathbf {H}_\text {eff}^\text {(non-exch)}), \label {eq::TI_fE}
\end{align}
where \(\mathbf {H}_\text {eff}^\text {(non-exch)} = \mathbf {H}_\text {eff} - \mathbf {H}_\text {exch}\). Crucially, only the damping component of exchange enters \(f_I\), while the precessional
term \(-\gamma _L \mathbf {m} \times \mathbf {H}_\text {exch}\) remains in \(f_E\). This choice ensures that the implicit Jacobian \(\partial f_I / \partial \mathbf {m}\) is symmetric positive-definite, enabling
efficient algebraic multigrid (AMG) or preconditioned conjugate gradient (PCG) solvers (cf. Section 5.4.2).
The ARKStep solver from SUNDIALS implements additive Runge–Kutta (ARK) methods for IMEX problems [161]. This approach is investigated in Section 5.5.
