BDF2-type integrator for Landau-Lifshitz-Gilbert equation in micromagnetics: IMEX time-stepping and full effective field

State of the art and related works

• Prior higher-order tangent-plane methods split weak-convergence and error analyses, while some IMEX approaches lacked the general field and torque coverage considered here.Existing methods could be first-order, formally almost second-order, require nonlinear solves, or fail to accommodate both contributions together.Source: §1 · p. 2
• The paper’s broader motivation is that much existing analysis uses exchange-only or exchange-plus-applied-field models, whereas micromagnetic energies can also include anisotropy, stray field, and bulk DMI.The introduction notes DMI’s role in favoring twisted textures and stabilizing chiral structures such as skyrmions.Source: §1 · pp. 1–2

Contributions of the present work

• Earlier BDF2 work [AFP26a] proved weak convergence for exchange and Zeeman fields, but its energy inequality used an endpoint restriction such as τ = o(h^2).The restriction entered only in estimates of initial and terminal correction terms; numerical experiments suggested it was an artifact, but its removal remained open.Source: §1 · p. 3
• The analysis removes the τ = o(h^2) restriction at the first and last time steps, proving unconditional weak subsequential convergence for arbitrary (h, τ) → (0, 0).The limit satisfies the modulus constraint and variational formulation of LLG; under strong torque consistency (B5), it also satisfies the physical energy inequality.Source: §1 · p. 3
• The proof uses the BDF2 G-stability quadratic energy and weak lower semicontinuity, avoiding control of the terminal exchange correction and a spatial inverse estimate.This replaces the endpoint-term estimates used in [AFP26a].Source: §1 · p. 3
• For the reduced exchange-plus-external-field model, the integrator coincides with [AFP26a], and the present analysis upgrades its convergence result to require no coupling between h and τ.The paper states that this case follows without changes to the proof.Source: §1 · p. 4
• A-priori analysis for the full effective field is deferred to subsequent work.That work is described as complementary to the present weak-convergence analysis.Source: §1 · p. 4

Energy contributions

• The total magnetic Gibbs free energy is a functional of magnetization and a possible applied external field, and the analysis covers contributions (a)–(e).The earlier works [AFP26a] restrict their energy to contributions (a) and (b).Source: §2 · p. 4
• The energy contributions include exchange, which penalizes magnetization non-uniformity, Zeeman energy, which penalizes deviation from an applied field, and anisotropy.The stated uniaxial anisotropy energy penalizes deviation from a preferred direction; other forms include planar and cubic anisotropy.Source: §2 · pp. 4–5

Magnetostatic contribution

• The stray-field energy is the magnetostatic self-energy from the interaction of magnetization with its own stray or demagnetizing field.The stray field is characterized through a full-space transmission problem for the magnetostatic potential.Source: §2 · p. 5

Dzyaloshinskii-Moriya interaction

• DMI favors twisting magnetization, unlike exchange, which favors alignment, and is described as fundamental to forming magnetic skyrmions.The analysis restricts to bulk DMI, with its coefficient depending on the material’s crystallographic class.Source: §2 · pp. 5–6

3. Numerical BDF2 integrator and unconditional weak convergence

• The work is a complementary extension of [AFP26a], providing a more complete description of the model.The results and proofs extend counterparts in those papers.Source: §3 · p. 6
• Setting λDM = 0 recovers the case without DMI, while setting Q = 0 recovers the case without non-conservative torque.These reductions recover settings covered by the earlier work.Source: §3 · p. 6
• The analysis avoids any CFL-type coupling in space-time.This is described as an improved analysis accompanying the complementary extension.Source: §3 · p. 6

Model setting

• The setting is a bounded Lipschitz ferromagnetic domain in two or three dimensions, with the two-dimensional case corresponding to a thin-film limit.The physical stray-field problem is in three dimensions; in two dimensions the abstract analysis applies to reduced lower-order operators satisfying (A1).Source: §3.1 · p. 6

Effective field and torque

• The effective field is the negative variational derivative of the energy; its lower-order anisotropy and stray-field contributions are represented by π(m).The torque Q is not assumed to arise from the energy and is treated as a general non-conservative forcing term.Source: §3.1 · pp. 6–7

Assumptions

• The analysis covers anisotropy contributions whose effective field is linear in m, such as uniaxial or planar anisotropy, but not cubic anisotropy under (A1).The torque assumption can cover first-order terms depending on ∇m and zero-order terms depending only on m, after suitable extension if needed.Source: §3.1 · p. 7

LLG model

• The Gilbert-form LLG equation includes the effective field and a non-conservative torque, with initial unit-length magnetization and a chiral boundary condition when DMI is present.The damping parameter α is positive, and the magnetization’s unit modulus follows from pointwise orthogonality of the equation’s terms to m.Source: §3.1 · p. 7

3.2. Discretization.

• The time interval is uniformly partitioned with step size τ, and the sequence of discrete functions is equipped with first-order and BDF2 discrete time derivatives.The discretization also defines nodal weights and uses standard nodal norm equivalence.Source: §3.2 · p. 8
• The spatial mesh family is quasi-uniform and regular, with continuous piecewise-affine vector fields on triangular or tetrahedral elements in dimensions two or three.The mesh size is h, with element size and diameter controlled by h and the quasi-uniformity parameter.Source: §3.2 · p. 8
• Admissible discrete magnetizations satisfy the modulus constraint at mesh vertices, while the discrete tangent space enforces vertexwise orthogonality to the magnetization.These constructions mimic the continuous unit-modulus constraint and the property m·∂_t m = 0.Source: §3.2 · pp. 8–9

3.3. Numerical integrator.

• The discrete torque is reconstructed piecewise constantly in time, and the applied field is sampled at the next time level and interpolated at mesh nodes.The field reconstruction converges to the continuous field as h and τ tend to zero.Source: §3.3 · p. 9
• Assumptions (B1)–(B4) support stability and convergence, while the additional strong sequential consistency assumption (B5) is used only for the energy inequality.The discrete lower-order operator is assumed linear, uniformly bounded, and consistent with its continuous counterpart.Source: §3.3 · p. 9
• The scheme uses BDF2 differences for the tangent-space variational formulation, with exchange and DMI terms integrated by parts using the chiral boundary condition.Explicitly treated magnetization-dependent contributions use a second-order extrapolation; the same extrapolation is used for DMI and non-conservative torque contributions.Source: §3.3 · p. 10
• The resulting IMEX BDF2-type scheme preserves a linear discrete problem because the possibly nonlinear torque operator is evaluated at already known magnetizations.A non-projected first-order tangent-plane step supplies the additional initial value required by BDF2.Source: §3.3 · p. 10
• The bilinear form for the BDF2 steps is continuous and coercive for every τ > 0, so the Lax–Milgram lemma guarantees a unique solution.The first-step bilinear form is coercive under its stated condition on τ.Source: §3.3 · p. 10

3.4. Physical examples of non-conservative torque Q.

• The cutoff R is the metric projection onto the closed unit ball, is nonexpansive, and equals its input whenever that input has norm at most one.These properties provide a bounded extension for the discrete torque inputs.Source: §3.4 · p. 11
• The Zhang–Li torque combines a directional derivative along the current density with a term proportional to that derivative, using the cutoff R on the discrete magnetization.The cutoff leaves the continuous torque unchanged on unit-length magnetizations, but strong consistency is not generally obtained from the available compactness because the torque depends on ∇m.Source: §3.4 · p. 11
• The Slonczewski-type torque uses a polarization field and bounded Lipschitz functions of the magnetization-polarization inner product.Its damping-like and field-like contributions are contained in −ϕ × Q_SL(ϕ).Source: §3.4 · p. 11

3.5. Convergence theorem.

• The theorem’s weak solution has H1 space-time regularity, bounded-in-time H1 regularity, unit modulus almost everywhere, and the prescribed initial trace.These are among the defining properties of weak solutions used in the theorem.Source: §3.5 · p. 11
• For any (h, τ) → (0, 0) without a coupling condition, the discrete approximations admit subsequences converging to a weak solution of the full LLG model.The weak-solution framework is adapted to include DMI and general non-conservative torque contributions.Source: §3.5 · pp. 11–13
• The limit weak solution satisfies the energy inequality when the additional strong sequential consistency assumption (B5) holds.That assumption is invoked to pass to the non-conservative work term in the energy inequality.Source: §3.5 · p. 13
• If the weak solution is unique in the considered regularity class, the convergence statements hold for the full sequences rather than only subsequences.The full-sequence conclusion applies as h and τ tend to zero.Source: §3.5 · p. 13

3.6. Discrete energy identity.

• DMI contributions require a separate argument using the symmetry of H_DM, while the torque contribution is treated under assumption (B4).The lower-order field, DMI, and torque terms add contributions to the energy analysis.Source: §3.6 · p. 14
• The analysis extends the BDF2 energy argument using the associated G-norm.The G-norm is introduced in the lemma extending [AFP26a, Lemma 3.5].Source: §3.6 · p. 14
• Algorithm A provides a discrete energy identity alongside uniform bounds for the discrete approximations.The estimates extend the corresponding results for the reduced-field setting.Source: §3.6 · pp. 14, 18
• With an additional regularity bound independent of h and τ, the L1 modulus-constraint defect is bounded by τ^2 + h^2.The quadratic temporal decay is consistent with the behavior observed numerically.Source: §3.6 · p. 20
• The interpolants are uniformly bounded by estimates that do not depend on the explicit form of the effective field or non-conservative torque.The proof follows from the discrete bounds and interpolant definitions.Source: §3.6 · p. 20

3.7. Extraction of weakly convergent subsequences.

• A subsequence converges to a limit m in H1(ΩT) ∩ L∞(0, T; H1(Ω)), with |m| = 1 almost everywhere.This is the existence and unit-modulus conclusion of Proposition 11.Source: §3.7 · p. 21
• The interpolants converge strongly to m in Hs(ΩT) for 0 < s < 1 and in L2(0, T; Hs(Ω)) for 0 < s < 1.The compact embedding H1(ΩT) ⋐ Hs(ΩT) is used to obtain strong convergence.Source: §3.7 · p. 21
• The extracted weak and weak-star limits are identified with m through strong L2(ΩT) convergence.The continuous embedding of L2(0, T; H1(Ω)) into L2(ΩT) supports this identification.Source: §3.7 · pp. 21–22
• The time derivatives converge weakly in L2(ΩT) to the time derivative of m.The proof identifies the weak limit of the discrete time derivatives with ∂tm.Source: §3.7 · p. 22

3.8. Identification of the limit with a weak solution of LLG.

• The limit has initial trace m(0) = m0, verifying the initial-data condition for a weak solution.The proof uses strong convergence of the initial discrete data in H1(Ω).Source: §3.8 · p. 23
• The discrete variational formulation converges to its continuous counterpart, first for smooth tests and then for all H1(ΩT) tests by density.This establishes the variational identity required for the weak solution.Source: §3.8 · p. 25
• Together with the limit's regularity and unit-modulus property, these steps establish that m is a weak solution of LLG.The section concludes the proof of Theorem 3 by establishing the defining properties of a weak solution.Source: §3.8 · p. 25
• Under assumption (B5), the limit also satisfies the energy inequality.Assumption (B5) is used for the energy inequality, not for identifying the weak limit.Source: §3.8 · pp. 25, 30
• If the weak solution is unique in the considered regularity class, the convergences hold for the full sequences rather than only subsequences.The conclusion applies as (h, τ) → (0, 0).Source: §3.8 · p. 30

4.1. Experimental convergence rates and stability of an isolated skyrmion.

• The experiments adapt [HPP+19, Section 4.2] by replacing its disk with a square, using bulk instead of interfacial DMI, and adopting a relaxed thin-film approximation.The reduced model is evaluated for isolated-skyrmion stability.Source: §4 · p. 30
• Skyrmion-like spatial errors decrease irregularly rather than at expected orders because the discontinuous initial state violates H1 and hypothesis (B1), while mesh-dependent transition shapes contribute.Successive-mesh comparisons also remain irregular and do not establish improved accuracy or convergence.Source: §4.1 · pp. 35–36, Figure 4 · p. 36, Figure 5 · p. 36
• At Δt6 = 3.125 · 10−3 ps, BDF2’s unit-length deviation is about two orders of magnitude smaller than pfTPS’s; their deviation convergence orders are second and first, respectively.The comparison uses the same first time step for both schemes.Source: §4.1 · pp. 36–37
• pfTPS deviation tracks dissipated energy, whereas BDF2 accumulates squared second differences and is more accurate when the time step is small relative to magnetization changes.For the smallest step, BDF2’s start-up term accounts for more than 90% of final deviation; for the largest, it accounts for 60% and 26% for the two initial states.Source: §4.1 · p. 37
• Under bulk DMI, the skyrmion-like reversed core disappears at D = 1–2 mJ m−2, persists and grows at D = 3–4 mJ m−2, and more complex patterns appear for D ≥5 mJ m−2.The in-plane magnetization has a Bloch-type structure near the domain wall; stronger DMI promotes twisting and more complex domains.Source: §4.1 · p. 38, Figure 8 · p. 39
• With interfacial DMI, isolated reversed regions occur for D = 3–5 mJ m−2, deform at D = 6 mJ m−2, and give way to complex patterns at larger D.The walls are Néel-type; unlike, approximately uniform states persist through D ≤5 mJ m−2, and the comparison is qualitative because of model differences and the square boundary.Source: §4.1 · pp. 39–40, Figure 9 · p. 40

4.2. A benchmark problem: relaxation of a bulk-DMI skyrmion.

• The relaxed bulk-DMI profile remains close to the expected Bloch-type state, with nearly zero radial magnetization at the skyrmion radius.During energy-dissipating relaxation, the skyrmion adjusts its profile and radius while preserving its Bloch-type character.Source: §4.2 · p. 41
• The computed bulk-DMI skyrmion radius is 21.72 nm, close to the 22.03 nm analytical benchmark (1.4% discrepancy) and numerical values of 21.87 nm and 22.1 nm.The nearly vanishing radial component at the skyrmion radius also confirms the expected Bloch-type chirality.Source: Figure 10 · p. 41, §4.2 · p. 41

4.3. Current-driven steering of a Bloch skyrmion.

• Under repeated current-driven direction changes, the skyrmion remains localized and is steered without collapse or expulsion over the simulated interval.The tracked trajectory displays the resulting steering in response to the applied current.Source: §4.3 · p. 43, Figure 12 · p. 43

• The auxiliary results collected here support the paper’s proofs, while their detailed proofs are available in [AFP26a].Source: Appendix A · p. 43
• Lemma A.1 provides norm-equivalence results for the discrete approximations.Source: Appendix A · p. 44
• The appendix records useful identities satisfied by the discrete approximations generated by the BDF2 update.Source: Appendix A · p. 44
• Lemma A.3 states a discrete Gronwall inequality for sequences.Source: Appendix A · p. 45