Semi-implicit Structure Preserving Method for The Landau-Lifshitz Equation

1. Problem Statement

The paper addresses the numerical simulation of the Landau-Lifshitz-Gilbert (LLG) equation, which governs the temporal evolution of magnetization in ferromagnetic materials. The LLG equation is a vector-valued nonlinear partial differential equation subject to a strict pointwise norm constraint ($|m| = 1$), representing the constant magnitude of the magnetization vector.

Key Challenges Identified:

• Projection Deficiencies: Existing numerical methods often rely on a "projection step" to enforce the norm constraint after each time step. While this preserves the magnitude, the projection is inherently nonlinear. This introduces significant difficulties in rigorous theoretical analysis regarding the stability and convergence of the scheme.
• Stability vs. Complexity: Fully implicit methods are stable but computationally expensive due to nonlinear solvers. Explicit methods are cheap but suffer from strict CFL (Courant-Friedrichs-Lewy) stability constraints.
• Lack of Rigorous Justification: Many structure-preserving schemes lack rigorous theoretical proofs for the stability of their normalization steps.

2. Methodology

The authors propose a semi-implicit numerical scheme that avoids artificial projection steps by exploiting the intrinsic structure of the LLG equation. The method combines two specific iterative procedures:

A. Governing Equation Reformulation

The study utilizes the standard LLG equation form:
$$m_t = -m \times h_{eff} - \alpha m \times (m \times h_{eff})$$
where $h_{eff} = \Delta m$ (under the simplified exchange energy functional). The authors explicitly choose this form over alternative formulations (like the heat flow of harmonic maps) to retain the original physical structure.

B. The Proposed Two-Step Scheme

The method is a first-order accurate scheme in time and second-order in space, consisting of two distinct steps:

1. Stable Evolution Step (Semi-implicit BDF1):
   A predictor step calculates an intermediate solution $\tilde{m}^{n+1}$ using a semi-implicit Backward Differentiation Formula (BDF1). This step handles the damping and gyromagnetic terms stably without enforcing the norm constraint yet.
   $$\frac{\tilde{m}^{n+1}_h - m^n_h}{k} = -m^n_h \times \Delta_h \tilde{m}^{n+1}_h - \alpha m^n_h \times (m^n_h \times \Delta_h \tilde{m}^{n+1}_h)$$

2. Norm-Preserving Step (Crank-Nicolson Type):
   A corrector step updates the solution to $m^{n+1}_h$ using a Crank-Nicolson-type formulation involving the cross product with the average of the current and new states. This structure naturally preserves the norm without an explicit projection:
   $$\frac{m^{n+1}_h - m^n_h}{k} = -\left(\frac{m^{n+1}_h + m^n_h}{2}\right) \times \Delta_h \tilde{m}^{n+1}_h - \alpha \left(\frac{m^{n+1}_h + m^n_h}{2}\right) \times (m^n_h \times \Delta_h \tilde{m}^{n+1}_h)$$

Key Theoretical Insight:
By applying the Crank-Nicolson structure to the cross-product term, the scheme ensures that the inner product of the update with the vector sum $(m^{n+1} + m^n)$ vanishes, mathematically guaranteeing $\|m^{n+1}\| = \|m^n\|$. This eliminates the need for a nonlinear projection operator.

3. Key Contributions

• Structure-Preserving without Projection: The primary contribution is a scheme that enforces the unit-norm constraint naturally through the discretization structure, rather than via an external, nonlinear projection step.
• Theoretical Advantages: By removing the projection step, the method avoids the nonlinearity-induced complications in stability and convergence analysis that plague traditional projection methods. It guarantees the uniqueness of the numerical solution.
• Computational Efficiency: The scheme results in a linear system at each time step (unlike fully implicit nonlinear solvers), making it computationally efficient while maintaining unconditional stability (in the sense of norm preservation) and avoiding strict CFL constraints associated with explicit schemes.
• Accuracy: The method achieves first-order accuracy in time and second-order accuracy in space.

4. Numerical Results

The authors conducted extensive numerical experiments in 1D and 3D settings to validate the method:

• Accuracy Verification:
  • Temporal Accuracy: Using exact analytical solutions, the scheme demonstrated first-order convergence in time ($O(k)$) for both 1D and 3D cases across $L_\infty$, $L_2$, and $H_1$ norms.
  • Spatial Accuracy: With a fixed small time step, the scheme demonstrated second-order convergence in space ($O(h^2)$).
• Norm Preservation:
  • Numerical tests confirmed that the deviation $\| |m_h| - 1 \|$ remained at the level of machine precision ($\approx 10^{-15}$) across various time steps and grid sizes, validating the theoretical norm-preserving property.
• Stability and Robustness:
  • The method was tested against various initial conditions (smooth and complex modulations) and compared with the standard BDF1 projection method.
  • The proposed method showed superior stability and consistency, maintaining the physical structure of the magnetization vector without the oscillations or drifts sometimes seen in projection-based methods.

5. Significance and Future Work

Significance:
This work provides a robust alternative to projection-based methods for micromagnetic simulations. By deriving the scheme directly from the LLG equation's structure, it offers a mathematically cleaner framework for theoretical analysis. It is particularly valuable for simulations requiring long-time stability and strict adherence to physical constraints (unit magnetization) without the computational overhead of nonlinear solvers.

Future Directions:
The authors outline several avenues for future research:

1. Rigorous Analysis: Conducting a full theoretical proof of stability and convergence for the proposed scheme.
2. Spatial Discretization: Extending the method to use the Finite Element Method (FEM) to handle complex geometric domains.
3. Higher-Order Schemes: Generalizing the approach to higher-order temporal accuracy using BDF-$k$ schemes ($k=2,3,4,5$).
4. Practical Application: Validating the method in complex, real-world micromagnetic modeling scenarios.