An unconditionally convergent BDF2-type finite element integrator for the inertial Landau-Lifshitz-Gilbert equation

1. Introduction

• The iLLG equation extends classical LLG with a second-time-derivative term associated with rapid nutation superimposed on slower damped precession.The magnetization has a unit-length constraint, making the problem nonconvex and highly nonlinear.Source: §1 · p. 1

1.1. Literature review.

• Existing iLLG finite-element methods include a linear TPS with damping and mesh restrictions, and nonlinear AMM and midpoint methods requiring nonlinear solves.The TPS and AMM converge to weak solutions; the midpoint scheme is not analyzed in.Source: §1.1 · pp. 1–2
• Unique solvability of nonlinear systems in cited methods requires CFL-type conditions, ranging from k = O(h) for AMM to k = O(h^2) for the midpoint scheme.The cited methods use fixed-point iterations to ensure solvability.Source: §1.1 · p. 2
• Projection-free TPS variants avoid weakly acute mesh requirements but incur unit-length error proportional to the time step; cited gradient-flow methods reach second-order error under discrete time regularity.These methods concern LLG or harmonic maps rather than the proposed iLLG integrator.Source: §1.1 · p. 2

1.2. Contributions.

• The proposed integrator combines a BDF2 magnetization update and tangent-space extrapolation with a β-weighted BDF2/backward Euler treatment of inertia.It is well posed, has a mimetic discrete energy law, and is unconditionally stable without mesh/time-step coupling; for β ∈ (0, 1], it gives unconditional second-order unit-length accuracy, while third-order accuracy requires discrete time regularity.Source: §1.2 · p. 2
• The discrete energy inequality requires no CFL restriction, and experiments report unconditional second-order unit-length accuracy, conditional third-order accuracy, and numerical evidence of second-order time accuracy in H1 for β = 0.The third-order constraint estimate requires a discrete time-regularity condition.Source: §1.2 · p. 3

1.3. Outline.

• The paper proceeds from the iLLG model and preliminaries to the algorithm, numerical experiments, and proofs.The experiments address constraint order, time-accuracy order, and the effect of β.Source: §1.3 · p. 3

1.4. General notation.

• The notation section fixes conventions for norms and generic constants, with constants independent of discretization parameters.The shorthand A ≲ B denotes an inequality up to such a constant, and A ≃ B denotes two-way inequality.Source: §1.4 · p. 3

• The iLLG model governs unit-length magnetization and adds an inertial second-time-derivative term that produces nutation alongside damped precession.The angular momentum relaxation time is τ, and α is the Gilbert damping parameter.Source: §2 · pp. 3–4
• The initial data impose |m0| = 1 and m0 · v0 = 0, and the magnetization magnitude is conserved in time.The initial velocity is specified by ∂tm(0) = v0.Source: §2 · pp. 3–4
• Unlike LLG, iLLG includes kinetic energy τ ∥∂tm∥2 /2; total magnetization energy comprises potential energy E[m] plus this term, and sufficiently smooth solutions satisfy a dissipative energy law.The effective field is the negative Gâteaux derivative of the energy functional and drives magnetization dynamics, responsible for precession.Source: §2 · p. 4
• A global weak solution belongs to L∞(0, ∞; H1(Ω)) ∩ W1,∞(0, ∞; L2(Ω)) and satisfies the unit-length constraint and initial data.Its definition also includes a variational formulation and weak energy law; integration by parts in time reduces regularity requirements.Source: §2 · pp. 4–5
• The weak formulation uses ∂t(m × ∂tm) = m × ∂ttm to express the inertial contribution after integration by parts in time.The identity follows because ∂tm × ∂tm vanishes.Source: §2 · p. 5

3.1. Time discretization.

• Time discretization uses uniform steps and difference quotients based on backward Euler and BDF2, with extrapolation buj+1 = 2uj − uj−1.The BDF2-adapted G-norm is defined for pairs of consecutive approximations using a positive-definite matrix.Source: §3.1 · pp. 5–6
• The G-norm is equivalent to a norm on consecutive solution pairs and satisfies identities used in the time-discrete analysis.The supplied passages specify the matrix entries and its positive eigenvalues.Source: §3.1 · p. 6

3.2. Space discretization.

• The spatial mesh is a quasi-uniform family of triangles or tetrahedra, with continuous piecewise-affine finite-element spaces and nodal interpolation.The vector-valued discrete magnetization is sought in the corresponding vector finite-element space.Source: §3.2 · p. 6
• Discrete magnetizations may violate unit length, with the stated L1 error controlled by δ; δ = 0 instead enforces unit length at mesh vertices.The cited projection-free algorithms have δ = O(k), while nodal projection or nonlinear constraint-preserving discretizations can produce the δ = 0 set.Source: §3.2 · pp. 6–7
• The tangent-plane construction targets the iLLG orthogonality condition m · ∂tm = 0, but imposes the discrete orthogonality property only at mesh vertices.The approximate time derivative is sought in the discrete tangent space.Source: §3.2 · p. 7

Algorithm and main results

• Algorithm A uses a BDF2 magnetization update and splits the inertial-term discretization between backward Euler and BDF2 with weight β.The first-order component supplies a bound needed in the convergence analysis; the BDF2 velocity discretization reduces dissipation but makes that analysis more involved.Source: §4 · pp. 7–8
• Each algorithm step is well posed, and its discrete energy law combines α-modulated energy dissipation with artificial dissipation and unsigned time-discretization perturbations.Source: §4 · p. 8
• For β in (0, 1], the scheme has an unconditional second-order unit-length constraint bound and convergence analysis; for β = 0, only the first-order bound is proved unconditionally.The third-order constraint bound requires a discrete regularity condition not guaranteed by the scheme; experiments report second-order H1-norm convergence for β = 0, which the analysis does not cover.Source: §4 · pp. 9–10
• For β in (0, 1], subsequences of the discrete functions converge unconditionally to a weak solution as h and k tend to zero, without a CFL-type condition.Source: §4 · p. 10

5.1. Unit length constraint.

• Second- and third-order unit-length constraint accuracy is observed for Algorithm A, including at β = 0; pfTPS1 is only first-order, while the τ = 0 LLG algorithm does not attain third-order accuracy.The third-order bound is linked to uniformly bounded auxiliary quantities A and B; the experiments show A and B remain bounded, while B approaches about 1.44×10^5 and β = 0 converges faster than larger β values.Source: §5.1 · pp. 11–12, Figure 1 · p. 11

5.2. Effects of the β-parameter.

• Larger β increases dissipation; Algorithm A's energy evolution need not decrease monotonically, while pfTPS1 and β near 1 produce rapidly decreasing magnetization amplitudes.The additional dissipation is most pronounced for pfTPS1, which has a distinct, fully monotone-decreasing energy evolution and the fastest amplitude decrease.Source: §5.2 · p. 12, Figure 3 · p. 13, Figure 4 · p. 13

5.3. Convergence order in time.

• Second-order H1-error convergence persists only for small β, with the drop toward first order depending on the initial condition.It persists through β = 10^-3 for the constant initial condition and through β = 10^-2 for the high-energy condition; the authors characterize any β > 0 as formally first-order.Source: §5.3 · pp. 13–14, Figure 5 · p. 14

5.4. Finite-time blow-up of weak solutions.

• For the pfTPS1/full BDF2 (β=0) comparison, with k=hℓ/10 and ℓ∈{4,5,6,7}, at t≈0.3 the gradient norm reaches values close to 2^(ℓ+1), and apparent blow-up timing is mesh-independent.The displayed value can exceed 2^(ℓ+1) because the schemes do not enforce unit length at the nodes.Source: §5.4 · p. 15

5.5. Self-consistency.

• After initialization, β = 0 quickly brings the effective Gilbert parameter near α = 0.1, whereas pfTPS1 has larger overall damping and approaches α more gradually.Increasing β raises the effective parameter, though it remains below the pfTPS1 value; the initialization value hovers slightly above 160.1.Source: §5.5 · p. 19, Figure 10 · p. 19

6.1. Well-posedness, discrete energy law, and stability.

• β = 1 is analytically excellent because it guarantees the needed boundedness of k2 Pn−1, while β = 0, the BDF2 discretization of the inertial term, performs better numerically, especially for large time steps in energy conservation.Source: §6.1 · p. 22
• Algorithm A is unconditionally stable, and its H1(Ω) estimate holds for all β ∈ [0,1] with a hidden constant depending only on α and τ, uniformly in β.The analysis gives explicit unit-length constraint-error bounds: the first-order estimate holds for all β ∈ [0,1] with C2 independent of β, while the sharper estimate for β ∈ (0,1] has C3 depending on β^-1.Source: §6.1 · pp. 23, 25, 26

6.2. Proof of Theorem 4.5.

• Uniform bounds for the time reconstructions yield a subsequential limit m with |m| = 1 almost everywhere as (h, k) → (0, 0).The limit has regularity m ∈ L∞(0, ∞; H1(Ω)) ∩ W1,∞(0, ∞; L2(Ω)); the constraint follows by passing to the limit in the estimates for β = 0 and β > 0.Source: §6.2 · p. 26
• Passing the summed discrete variational equations to the limit establishes the weak formulation, with the inertial-term contributions requiring separate analysis.The proof distinguishes the three LLG terms from the first three terms arising from the iLLG inertial term.Source: §6.2 · p. 27
• The discrete variational formulation converges to its continuous counterpart, so the limit satisfies Definition 2.1(iii).This completes the proof that the limit satisfies that part of the weak-solution definition.Source: §6.2 · p. 32
• The limit satisfies the energy inequality for every t ∈ [0, T], extending the almost-everywhere result through weak time continuity.The proof establishes weak continuity of ∂tm in L2(Ω) before extending the inequality to all times.Source: §6.2 · pp. 34–35
• The initial conditions hold continuously as t → 0, with ∇m(t) → ∇m0 and ∂tm(t) → v0 strongly in L2(Ω).Together with m(t) → m0 in L2(Ω), this gives continuity in H1(Ω) × L2(Ω).Source: §6.2 · p. 36

Appendix A.

• The lemma assumes Banach spaces X and Y with continuous injection X ⊂ Y, reflexivity of X, and u ∈ L∞(0,T;X) weakly continuous in Y.Weak continuity in Y means every g ∈ Y* yields a continuous map t ↦ ⟨g,u(t)⟩ on [0,T].Source: Appendix A · p. 37
• Under these assumptions, u(t) belongs to X at every t ∈ [0,T] and satisfies ∥u(t)∥X ≤ ∥u∥L∞(0,T;X).Source: Appendix A · p. 37
• The lemma also establishes that u is weakly continuous in time with values in X.That is, every f ∈ X* yields a continuous map t ↦ ⟨f,u(t)⟩ on [0,T].Source: Appendix A · p. 37