An Efficient Energy Stable Structure Preserving Method for The Landau-Lifshitz Equation

One of the main difficulties in micromagnetics simulation is the norm preserving constraints $|\mathbf{m}|=1$ at the continuous or the discrete level. Another difficulty is the stability with the time step constraint. Using standard explicit integrators leads to a physical time step of sub-pico seconds, which is often two orders of magnitude smaller than the fastest physical time scales. Direct implicit integrators require solving complicated, coupled systems. Another major difficulty with the projection method in this field is the lack of rigorous theoretical guarantees regarding its stability of the projection step. In this paper, we introduce a first order method. Such a method is structure preserving based on a combination of a Gauss-Seidel iteration, a double diffusion iteration and a Crank-Nicolson iteration to preserve the norm constraints.

💡 Research Summary

The paper addresses two fundamental challenges in micromagnetic simulations of the Landau‑Lifshitz‑Gilbert (LLG) equation: the pointwise unit‑norm constraint ‖m‖ = 1 and the severe time‑step restriction that forces explicit schemes to use sub‑picosecond steps. Existing strategies—projection, Lagrange‑multiplier, and fully implicit methods—either lack rigorous stability guarantees for the projection step, introduce extra variables, or require solving large coupled nonlinear systems.

The authors propose a first‑order, structure‑preserving algorithm that simultaneously enforces the norm constraint and guarantees energy stability without imposing a restrictive CFL condition. The core idea is to combine three ingredients: a Gauss‑Seidel sweep for the nonlinear cross‑product terms, a double‑diffusion iteration that effectively inverts the Laplacian operator, and a Crank‑Nicolson‑type predictor‑corrector that yields an orthogonal update matrix.

Starting from the simplified, undamped LLG form m_t = − m × Δm, they apply a Crank‑Nicolson discretisation. The resulting linear system can be written as m^{n+1} = A m^{n} with A = (I + βC)^{-1}(I − βC), where C is a skew‑symmetric matrix and β = Δt^{2}. Because C^{T}=−C, A is orthogonal, which directly implies ‖m^{n+1}‖₂ = ‖m^{n}‖₂; thus the unit‑norm is preserved exactly at the discrete level.

To handle the spatial diffusion term Δm, three variants are introduced:

1. Scheme I (explicit) updates the Laplacian using the known solution at time n, leading to a classic CFL restriction Δt ≤ C h^{2}.

2. Scheme II (fully implicit) evaluates Δm at the new time level, requiring the solution of a nonlinear system; this is computationally expensive.

3. Scheme III (semi‑implicit) evaluates Δm at an intermediate state \tilde{m}^{n+1}. The intermediate state is obtained by a double‑diffusion iteration \tilde{m}^{n+1} = (I − k Δ)^{-1} m^{n}, where k is a small parameter. This step can be performed efficiently with FFT‑based solvers, reducing the cost to O(N log N) for N grid points. The semi‑implicit formulation retains the orthogonal property of the Crank‑Nicolson update, thereby preserving the norm without a restrictive time‑step condition.

The Gauss‑Seidel sweep updates the three components of m sequentially, which keeps the algorithm linear in each substep and maintains O(N) complexity. The double‑diffusion iteration supplies a stable approximation of the Laplacian inverse, mitigating the stiffness that typically forces tiny Δt in explicit schemes.

The authors extend the method to the full damped LLG equation m_t = − m × Δm − α m × (m × Δm). The damping term is treated explicitly in the predictor‑corrector stage, preserving the overall first‑order accuracy while still guaranteeing ‖m^{n+1}‖ = ‖m^{n}‖.

A comprehensive set of numerical experiments validates the theory. In 1‑D temporal convergence tests, halving Δt reduces the L₂ error by roughly a factor of two, confirming first‑order accuracy. Spatial convergence tests show that halving the mesh size improves the error by more than a factor of two, indicating that the double‑diffusion step enhances spatial accuracy beyond the nominal first order. Norm‑preservation tests report deviations on the order of 10^{-15}, i.e., machine precision. Energy decay is consistently larger for the proposed semi‑implicit scheme than for the classical Gauss‑Seidel projection method, and spurious oscillations are markedly reduced.

When a modest damping coefficient (α = 0.01) is introduced, the algorithm remains stable over long integration times, and the error growth remains consistent with first‑order behavior.

In summary, the paper delivers a mathematically rigorous, computationally efficient, and practically robust first‑order scheme for the LLG equation. By embedding the unit‑norm constraint into the time‑integration formula via an orthogonal update matrix, and by handling the stiff diffusion term through a cheap double‑diffusion iteration, the method overcomes the principal limitations of existing projection‑based approaches. The approach is well‑suited for large‑scale three‑dimensional micromagnetic simulations, spin‑tronic device modeling, and any application where long‑time dynamics of magnetization must be captured without sacrificing stability or physical fidelity.