The Mass-Lumped Midpoint Scheme for Computational Micromagnetics: Newton Linearization and Application to Magnetic Skyrmion Dynamics

Giovanni Di Fratta; Carl-Martin Pfeiler; Dirk Praetorius; Michele Ruggeri

doi:10.1515/cmam-2022-0060

Open Access Published by De Gruyter June 9, 2022

The Mass-Lumped Midpoint Scheme for Computational Micromagnetics: Newton Linearization and Application to Magnetic Skyrmion Dynamics

Giovanni Di Fratta , Carl-Martin Pfeiler , Dirk Praetorius and Michele Ruggeri

From the journal Computational Methods in Applied Mathematics

https://doi.org/10.1515/cmam-2022-0060

Abstract

We discuss a mass-lumped midpoint scheme for the numerical approximation of the Landau–Lifshitz–Gilbert equation, which models the dynamics of the magnetization in ferromagnetic materials. In addition to the classical micromagnetic field contributions, our setting covers the non-standard Dzyaloshinskii–Moriya interaction, which is the essential ingredient for the enucleation and stabilization of magnetic skyrmions. Our analysis also includes the inexact solution of the arising nonlinear systems, for which we discuss both a constraint-preserving fixed-point solver from the literature and a novel approach based on the Newton method. We numerically compare the two linearization techniques and show that the Newton solver leads to a considerably lower number of nonlinear iterations. Moreover, in a numerical study on magnetic skyrmions, we demonstrate that, for magnetization dynamics that are very sensitive to energy perturbations, the midpoint scheme, due to its conservation properties, is superior to the dissipative tangent plane schemes from the literature.

Keywords: Landau–Lifshitz–Gilbert Equation; Dzyaloshinskii–Moriya Interaction; Magnetic Skyrmions; Newton Linearization; Computational Micromagnetics; Finite Elements

MSC 2010: 35K55; 65M12; 65M22; 65M60; 65Z05

1 Introduction

1.1 Energetics of a Ferromagnet

In the continuum theory of micromagnetism, whose origin dates back to the seminal work of Landau–Lifshitz [34] on small ferromagnetic particles, the amount of magnetic moment (per unit volume) of a rigid ferromagnetic body occupying a bounded region Ω ⊂ R 3 is represented by a classical vector field, the magnetization M : Ω → R 3 . Its module M s := | M | describes the so-called saturation magnetization. In single-crystal ferromagnets [2, 6], M s depends only on the temperatureand is assumed to be constant when the specimen is well below the so-called Curie temperature of the material. In this case, the magnetization can be represented in the form M := M s ⁢ m , where m : Ω → S 2 is a vector field with values in the unit sphere of R 3 , and the observable magnetization states minimize the micromagnetic energy functional [16, 31]

(1.1) E ⁢ ( m ) := E Ω ⁢ ( m ) + K Ω ⁢ ( m ) + W Ω ⁢ ( m ) + A Ω ⁢ ( m ) + Z Ω ⁢ ( m ) := ∫ Ω A ⁢ | ∇ ⁡ m | 2 + D ⁢ ( ∇ × m ) ⋅ m - μ 0 2 ⁢ M s ⁢ H s ⁢ ( m ) ⋅ m + φ an ⁢ ( m ) - μ 0 ⁢ M s ⁢ H ext ⋅ m ⁢ d ⁢ x ,

defined for every m ∈ H 1 ⁢ ( Ω ; S 2 ) .

The exchange energy E Ω ⁢ ( m ) penalizes spatial variations of the direction of the magnetization, with A > 0 representing a material-dependent constant that summarizes the stiffness of short-range (symmetric) exchange interactions. The second term K Ω ⁢ ( m ) represents the bulk Dzyaloshinskii–Moriya interaction (DMI) [23, 35] and accounts for antisymmetric exchange interactions caused by possible lack of inversion symmetry in the crystal structure of the ferromagnet. The sign of the constant D ∈ R affects the chirality of the ferromagnetic system [43, 42]. The third term W Ω ⁢ ( m ) is the magnetostatic self-energy, i.e., the energy due to the stray field H s ⁢ ( m ) induced by M s ⁢ m . From the mathematical point of view, H s ⁢ ( m ) can be characterized as the projection of ( - M s ⁢ χ Ω ⁢ m ) ∈ L 2 ⁢ ( R 3 ; R 3 ) on the closed subspace of gradient vector fields ∇ ⁡ H 1 ⁢ ( R 3 ; R ) := { ∇ ⁡ u : u ∈ H 1 ⁢ ( R 3 ; R ) } (see, e.g., [38, 21]). Here, μ 0 denotes the vacuum permeability. Additionally, the micromagnetic energy includes two additional energy contributions: the magnetocrystalline anisotropy energy A Ω ⁢ ( m ) and the Zeeman energy Z Ω ⁢ ( m ) . The energy density φ an : S 2 → R ≥ 0 models the existence of easy directions of the magnetization due to the crystallographic structure of the ferromagnet, while Z Ω ⁢ ( m ) models the tendency of a specimen to have the magnetization aligned with the external applied field H ext ∈ L 2 ⁢ ( Ω ; R 3 ) , assumed to be unaffected by variations of 𝒎. The competition among the energy contributions in (1.1) explains most of the striking spin textures observable in ferromagnetic materials [31], in particular, the emergence of magnetic skyrmions [25, 26].

1.2 A More General Energy Functional

When a ferromagnetic system consists of several magnetic materials, the material-dependent quantities 𝐴, 𝐷, and M s are no longer constant in the region Ω occupied by the ferromagnet, and one has to model spin interactions among different magnetic materials at their touching interface [2]. The easiest way is to assume a strong coupling condition [6, 5, 20]: although M s can be discontinuous across an interface, the direction of the magnetization never jumps through it. Under this constitutive assumption, the analysis of the composite can be carried out under the classical conditions M s ∈ L ∞ ⁢ ( Ω ; R > 0 ) and m ∈ H 1 ⁢ ( Ω ; S 2 ) . In this setting, the observable states of a rigid ferromagnetic body can be characterized as the local minimizers of the micromagnetic energy functional still defined by (1.1), but with the quantities A = A ⁢ ( x ) , D = D ⁢ ( x ) , and M s = M s ⁢ ( x ) to be understood as functions defined on Ω.

In this paper, we are interested in a more general energy functional which, other than incorporating the previous one as a special case, also accounts for the presence of anisotropies in the lattice structures of the constituents. To introduce the model, we first observe that the bulk DMI energy density can be equivalently rewritten as

D ⁢ ( ∇ × m ) ⋅ m = D ⁢ ∑ d = 1 3 ( e d × ∂ d ⁡ m ) ⋅ m ,

where { e d } d = 1 , 2 , 3 denotes the standard basis of R 3 . It is therefore a special case of the energy density

g asym ⁢ ( x , s , ξ ) = ∑ d = 1 3 K d ⁢ ( x ) ⁢ ξ d ⋅ s for all ⁢ x ∈ Ω , s ∈ R 3 , and ⁢ ξ = ( ξ 1 , ξ 2 , ξ 3 ) ∈ R 3 × 3 ,

with { K d } d = 1 , 2 , 3 being 3-by-3 antisymmetric matrices, i.e., K d = - K d T . Similarly, the symmetric exchange energy density can be generalized to the density

g sym ⁢ ( x , ξ ) = 1 2 ⁢ ∑ d = 1 3 A d ⁢ ( x ) ⁢ ξ d ⋅ ξ d for all ⁢ x ∈ Ω ⁢ and ⁢ ξ = ( ξ 1 , ξ 2 , ξ 3 ) ∈ R 3 × 3 ,

with { A d } d = 1 , 2 , 3 being 3-by-3 invertible symmetric matrices, i.e., A d = A d T . Hence, for g := g sym + g asym , it holds that

(1.2) g ⁢ ( x , s , ξ ) = 1 2 ⁢ ∑ d = 1 3 ( A d ⁢ ( x ) ⁢ ξ d ⋅ ξ d - 2 ⁢ K d ⁢ ( x ) ⁢ s ⋅ ξ d ) = 1 2 ⁢ ∑ d = 1 3 A d ⁢ ( x ) ⁢ ( ξ d - A d - 1 ⁢ ( x ) ⁢ K d ⁢ ( x ) ⁢ s ) ⋅ ( ξ d - A d - 1 ⁢ ( x ) ⁢ K d ⁢ ( x ) ⁢ s ) + 1 2 ⁢ ∑ d = 1 3 K d ⁢ ( x ) ⁢ A d - 1 ⁢ ( x ) ⁢ K d ⁢ ( x ) ⁢ s ⋅ s .

Note that K d ⁢ ( x ) ⁢ A d - 1 ⁢ ( x ) ⁢ K d ⁢ ( x ) is a symmetric matrix. This discussion suggests the opportunity to investigate an energy functional covering the above generalized form; see (2.1) below. It is worth to notice that the structure of this energy functional does not only allow for the description of a mixture of ferromagnetic materials, but also covers typical homogeneous models arising as Γ-limit of composite ferromagnetic materials with highly oscillating heterogeneities [5, 20].

1.3 Landau–Lifshitz–Gilbert Equation and Its Numerical Integration

When the magnetization 𝒎 does not minimize the micromagnetic energy functional, the ferromagnetic system is in a non-equilibrium state. A well-accepted model for its time evolution is the Landau–Lifshitz–Gilbert equation (LLG) [34, 28], which in the so-called Gilbert form reads

(1.3) ∂ t ⁡ m = - γ 0 ⁢ m × H eff ⁢ ( m ) + α ⁢ m × ∂ t ⁡ m .

This phenomenological equation describes the magnetization dynamics as a dissipative precession driven by the effective field H eff ⁢ ( m ) := - μ 0 - 1 ⁢ M s - 1 ⁢ δ ⁢ E ⁢ ( m ) δ ⁢ m , and modulated by the gyromagnetic ratio of the electron γ 0 > 0 and the Gilbert damping parameter α > 0 . The numerical approximation of LLG is not a trivial task. Nonlinearities, the numerical realization of the unit-length constraint, the possible coupling with other (nonlinear) partial differential equations, and the need of unconditionally stable numerical schemes make the problem very challenging. For this reason, in the last twenty years, the problem has been the subject of several mathematical studies; see, e.g., [40, 7, 33, 15, 27, 4, 14, 18, 19, 9, 8, 17, 1, 24, 32, 30, 22, 3].

In this work, we consider the mass-lumped midpoint scheme proposed in [15]. The method is based on a mass-lumped first-order finite element method for the spatial discretization and the second-order midpoint rule for the time discretization, and involves the solution of one nonlinear system per time step. Besides introducing the method, the work [15] proves unconditional convergence of the finite element approximation towards a weak solution of LLG in the sense of [10] and proposes a fixed-point iteration to linearize the nonlinear problem arising from the scheme. The scheme has also been the subject of further research. On the one hand, the works [12, 19] incorporate the inexact solution of the nonlinear system into the convergence result. On the other hand, the work [39] focuses on the design and the analysis of effective approaches to treat the nonlocal field contributions.

1.4 Contributions

In this work, as a novel contribution, we extend the midpoint scheme and its analysis to more general energy contributions; see the discussion in Section 1.2. In particular, the present analysis covers DMI, which is not covered by the analysis in [15, 12, 19, 39]. We note that DMI is the essential ingredient for the enucleation and the stabilization of magnetic skyrmions. At this point, it is worth pointing out that DMI contributions represent a challenging testing ground for numerical schemes for LLG. Indeed, besides requiring accurate adaptations in the numerical analysis, they determine magnetization configurations – magnetic skyrmions – that turn out to be very sensitive to small perturbations of the micromagnetic energy. In addition, we also discuss the linearization of the nonlinear scheme: we extend the fixed-point iteration proposed in [12] to the present setting and propose an approach based on the Newton method, for which we provide a first full analysis (well-posedness, stability, convergence). Finally, in a collection of numerical experiments, first, we accurately test the energy conservation properties of the mass-lumped midpoint scheme and extensively compare it with the tangent plane schemes from [4, 8, 30, 22]. Second, we emphasize the faster convergence of the Newton linearization over the fixed-point solver in a numerical study, which additionally indicates that the restrictive conditions required by our analysis of the Newton linearization are likely to be suboptimal and might be weakened to match those for the fixed-point solver by an improved analysis.

1.5 Outline

The remainder of the work is organized as follows. We conclude this section by collecting the notation used throughout the paper. In Section 2, we describe the mathematical problem under consideration. In Section 3, we present the proposed algorithms and state their stability and convergence results. Section 4 is devoted to numerical experiments. Finally, in Sections 5–6, we collect the proofs of the results stated in Section 3.

1.6 Notation

Throughout the paper, we use the standard notation for Lebesgue, Sobolev, and Bochner spaces and norms. To highlight (spaces of) vector-valued or matrix-valued functions, we use bold letters, e.g., we denote both L 2 ⁢ ( Ω ; R 3 ) and L 2 ⁢ ( Ω ; R 3 × 3 ) by L 2 ⁢ ( Ω ) . We denote by ⟨ ⋅ , ⋅ ⟩ Ω the scalar product in L 2 ⁢ ( Ω ) and by ⟨ ⋅ , ⋅ ⟩ the duality pairing between H 1 ⁢ ( Ω ) and its dual. By C > 0 , we always denote a generic constant, which is independent of the discretization parameters, but not necessarily the same at each occurrence.

2 Problem Formulation

Let Ω ⊂ R 3 be a bounded Lipschitz domain. The energy of m ∈ H 1 ⁢ ( Ω ; S 2 ) is given by

(2.1) E ⁢ ( m ) = 1 2 ⁢ a ⁢ ( m , m ) - ⟨ f , m ⟩ Ω ,

where f ∈ L 2 ⁢ ( Ω ) , while the bilinear form a : H 1 ⁢ ( Ω ) × H 1 ⁢ ( Ω ) → R is defined, for all ψ , ϕ ∈ H 1 ⁢ ( Ω ) , by

(2.2) a ⁢ ( ψ , ϕ ) = ∑ d = 1 3 ⟨ A d ⁢ ( ∂ d ⁡ ψ - J d ⁢ ψ ) , ∂ d ⁡ ϕ - J d ⁢ ϕ ⟩ Ω - ⟨ π ⁢ ( ψ ) , ϕ ⟩ Ω .

Here, π : L 2 ⁢ ( Ω ) → L 2 ⁢ ( Ω ) is a linear, bounded, and self-adjoint operator, while, for d = 1 , 2 , 3 , the 3-by-3 matrices A d and J d have coefficients in L ∞ ⁢ ( Ω ) , with A d being also symmetric and uniformly positive definite, i.e., it holds that A d T = A d and

A d ⁢ ( x ) ⁢ u ⋅ u ≥ A 0 ⁢ | u | 2 for almost all ⁢ x ∈ Ω ⁢ and all ⁢ u ∈ R 3 ,

where A 0 > 0 is a fixed constant. The energy (2.1) covers the extensions of the classical micromagnetic functional discussed in Section 1.2; cf. the expression in (1.2).

The existence of minimizers of (2.1) in H 1 ⁢ ( Ω ; S 2 ) follows from the direct method of calculus of variations. Moreover, any minimizer m ∈ H 1 ⁢ ( Ω ; S 2 ) satisfies the Euler–Lagrange equations

⟨ h eff ⁢ ( m ) , ϕ ⟩ = 0 for all ⁢ ϕ ∈ H 1 ⁢ ( Ω ) ⁢ such that ⁢ m ⋅ ϕ = 0 ⁢ a.e. in ⁢ Ω .

Here, h eff ⁢ ( m ) := - δ ⁢ E ⁢ ( m ) δ ⁢ m is the (negative) Gâteaux derivative of the energy, i.e.,

(2.3) - ⟨ h eff ⁢ ( m ) , ϕ ⟩ = ⟨ δ ⁢ E ⁢ ( m ) δ ⁢ m , ϕ ⟩ = lim δ → 0 ⁡ E ⁢ ( m + δ ⁢ ϕ ) - E ⁢ ( m ) δ ⁢ = ( ⁢ 2.1 ⁢ ) ⁢ a ⁢ ( m , ϕ ) - ⟨ f , ϕ ⟩ Ω .

Turning to the dynamical case, a non-equilibrium configuration m ⁢ ( t ) ∈ H 1 ⁢ ( Ω ; S 2 ) evolves according to (1.3), which, after a suitable rescaling, reads

(2.4) ∂ t ⁡ m = - m × h eff ⁢ ( m ) + α ⁢ m × ∂ t ⁡ m ,

with α > 0 being the Gilbert damping parameter. The dynamics is dissipative in the sense that any sufficiently smooth solution of (2.4) satisfies the energy law

(2.5) d d ⁢ t ⁢ E ⁢ ( m ⁢ ( t ) ) = - α ⁢ ∥ ∂ t ⁡ m ⁢ ( t ) ∥ L 2 ⁢ ( Ω ) 2 ≤ 0 for all ⁢ t > 0 .

We conclude this section by recalling the notion of a weak solution of (2.4); see [10].

Definition 2.1

Let m 0 ∈ H 1 ⁢ ( Ω ; S 2 ) . A vector field m : Ω × R > 0 → S 2 is called a global weak solution of (2.4) if m ∈ L ∞ ⁢ ( R > 0 ; H 1 ⁢ ( Ω ; S 2 ) ) and, for all T > 0 , with Ω T := Ω × ( 0 , T ) , the following properties are satisfied:

m ∈ H 1 ⁢ ( Ω T ) ;
m ⁢ ( 0 ) = m 0 in the sense of traces;
for all φ ∈ H 1 ⁢ ( Ω T ) , it holds that
(2.6) ∫ 0 T ⟨ ∂ t ⁡ m ⁢ ( t ) , φ ⁢ ( t ) ⟩ Ω ⁢ d t = - ∫ 0 T ⟨ h eff ⁢ ( m ⁢ ( t ) ) , φ ⁢ ( t ) × m ⁢ ( t ) ⟩ ⁢ d t + α ⁢ ∫ 0 T ⟨ m ⁢ ( t ) × ∂ t ⁡ m ⁢ ( t ) , φ ⁢ ( t ) ⟩ Ω ⁢ d t ;
it holds that
(2.7) E ⁢ ( m ⁢ ( T ) ) + α ⁢ ∫ 0 T ∥ ∂ t ⁡ m ⁢ ( t ) ∥ L 2 ⁢ ( Ω ) 2 ⁢ d t ≤ E ⁢ ( m 0 ) .

We note that (2.3) implicitly includes natural boundary conditions on 𝒎, which are homogeneous Neumann boundary conditions ∂ n = 0 if A d = ℓ ex 2 ⁢ Id and J d = 0 for d = 1 , 2 , 3 . For a more explicit presentation, we refer to [30]. The variational formulation (2.6) comes from a weak formulation of (2.4) in the space-time cylinder. The energy inequality (2.7) is a weak counterpart of the dissipative energy law (2.5).

Remark 2.2

(i) For ease of presentation, we restrict ourselves to the case of a time-independent field f ∈ L 2 ⁢ ( Ω ) . For time-dependent fields, the strong form (2.5) and the weak form (2.7) of the energy law of LLG read

d d ⁢ t ⁢ E ~ ⁢ ( m ⁢ ( t ) ) = - α ⁢ ∥ ∂ t ⁡ m ⁢ ( t ) ∥ L 2 ⁢ ( Ω ) 2 + ⟨ f ⁢ ( t ) , ∂ t ⁡ m ⁢ ( t ) ⟩ Ω

and

E ~ ⁢ ( m ⁢ ( T ) ) + α ⁢ ∫ 0 T ∥ ∂ t ⁡ m ⁢ ( t ) ∥ L 2 ⁢ ( Ω ) 2 ⁢ d t - ∫ 0 T ⟨ f ⁢ ( t ) , ∂ t ⁡ m ⁢ ( t ) ⟩ Ω ⁢ d t ≤ E ~ ⁢ ( m 0 ) ,

respectively, where E ~ ⁢ ( m ) = E ⁢ ( m ) + ⟨ f , m ⟩ Ω = a ⁢ ( m , m ) / 2 .

(ii) The present setting covers and generalizes the model problems considered in previous mathematical works on the numerical integration of LLG.

With the choices A d = ℓ ex 2 ⁢ Id and J d = 0 for d = 1 , 2 , 3 , where ℓ ex > 0 is the so-called exchange length and Id is the 3-by-3 identity matrix, π ≡ 0 , and f ≡ 0 , we obtain that
⟨ h eff ⁢ ( ψ ) , φ ⟩ = - ℓ ex 2 ⁢ ⟨ ∇ ⁡ ψ , ∇ ⁡ ϕ ⟩ Ω .
This is the so-called small particle limit of LLG, for which finite element schemes have been proposed, e.g., in the seminal papers [15, 4].
With the choices A d = ℓ ex 2 ⁢ Id and J d = 0 for d = 1 , 2 , 3 , we obtain that
⟨ h eff ⁢ ( ψ ) , ϕ ⟩ = - ℓ ex 2 ⁢ ⟨ ∇ ⁡ ψ , ∇ ⁡ ϕ ⟩ Ω + ⟨ π ⁢ ( ψ ) , ϕ ⟩ Ω + ⟨ f , ϕ ⟩ Ω .
This setting covers the classical energy contributions considered in micromagnetics (exchange, uniaxial anisotropy, magnetostatic, Zeeman) and numerical integrators for this case have been analyzed, e.g., in [9, 17, 8, 39, 22].
With the choices A d = ℓ ex 2 ⁢ Id and J d = ℓ dm ⁢ [ e d ] × / ( 2 ⁢ ℓ ex 2 ) for d = 1 , 2 , 3 , where ℓ dm ∈ R is a characteristic length associated with DMI, π ⁢ ( m ) = ℓ dm ⁢ m / ( 2 ⁢ ℓ ex 2 ) , and f ≡ 0 , we obtain that
⟨ h eff ⁢ ( ψ ) , ϕ ⟩ = - ℓ ex 2 ⁢ ⟨ ∇ ⁡ ψ , ∇ ⁡ ϕ ⟩ Ω - ℓ dm 2 ⁢ ⟨ ∇ × ψ , ϕ ⟩ Ω - ℓ dm 2 ⁢ ⟨ ψ , ∇ × ϕ ⟩ Ω ,
which is the setting analyzed in [30] for the simulation of chiral magnetic skyrmions by the means of a family of tangent plane integrators. Here, [ e d ] × ∈ R 3 × 3 denotes the matrix such that [ e d ] × ⁢ u = e d × u for all u ∈ R 3 .

3 Numerical Algorithms and Main Results

3.1 Preliminaries

Let κ ≥ 1 . For the spatial discretization, assuming Ω to be a polyhedral domain, we consider a 𝜅-quasi-uniform family { T h } h > 0 of regular tetrahedral triangulations of Ω parametrized by the mesh size h = max K ∈ T h ⁡ diam ⁡ ( K ) , i.e., κ - 1 h ≤ vol ( K ) 1 / 3 ≤ h for all K ∈ T h . We denote by N h the set of vertices of T h . For any K ∈ T h , we denote by P 1 ⁢ ( K ) the space of first-order polynomials on 𝐾. We consider the space of T h -piecewise affine and globally continuous finite elements

S 1 ⁢ ( T h ) := { v h ∈ C 0 ⁢ ( Ω ¯ ) : v h | K ∈ P 1 ⁢ ( K ) ⁢ for all ⁢ K ∈ T h } .

The classical basis for this finite-dimensional subspace of H 1 ⁢ ( Ω ) is the set of nodal hat functions { φ z } z ∈ N h , which satisfy φ z ⁢ ( z ′ ) = δ z , z ′ for all z , z ′ ∈ N h . The nodal interpolant I h : C 0 ⁢ ( Ω ¯ ) → S 1 ⁢ ( T h ) is defined by I h ⁢ [ u ] = ∑ z ∈ N h u ⁢ ( z ) ⁢ φ z for all u ∈ C 0 ⁢ ( Ω ¯ ) .

Let V h := S 1 ⁢ ( T h ) 3 . For each time step, approximate solutions of (2.4) are sought in the set of admissible approximate magnetizations

M h := { ϕ h ∈ V h : | ϕ h ⁢ ( z ) | = 1 ⁢ for all ⁢ z ∈ N h } ,

which consists of all elements of V h satisfying the unit-length constraint at the nodes of the triangulation.

Besides the standard scalar product ⟨ ⋅ , ⋅ ⟩ Ω , given a mesh T h and the associated nodal interpolant I h ⁢ [ ⋅ ] , we consider the mass-lumped product ⟨ ⋅ , ⋅ ⟩ h defined by

⟨ ψ , ϕ ⟩ h = ∫ Ω I h ⁢ [ ψ ⋅ ϕ ] ⁢ d x for all ⁢ ψ , ϕ ∈ C 0 ⁢ ( Ω ¯ ) .

Using the definition of the nodal interpolant, we see that

⟨ ψ , ϕ ⟩ h = ∑ z ∈ N h β z ⁢ ψ ⁢ ( z ) ⋅ ϕ ⁢ ( z ) , where ⁢ β z := ∫ Ω φ z ⁢ d x > 0 .

On V h , ⟨ ⋅ , ⋅ ⟩ h is a scalar product and the induced norm ∥ ⋅ ∥ h is equivalent to the standard norm of L 2 ⁢ ( Ω ) . In particular, it holds that

(3.1) ∥ ϕ h ∥ L 2 ⁢ ( Ω ) ≤ ∥ ϕ h ∥ h ≤ 5 ⁢ ∥ ϕ h ∥ L 2 ⁢ ( Ω ) for all ⁢ ϕ h ∈ V h ;

see [13, Lemma 3.9]. Finally, we define the mapping P h : H 1 ⁢ ( Ω ) ⋆ → V h by

(3.2) ⟨ P h ⁢ u , ϕ h ⟩ h = ⟨ u , ϕ h ⟩ for all ⁢ u ∈ H 1 ⁢ ( Ω ) ⋆ ⁢ and ⁢ ϕ h ∈ V h ,

i.e., P h ⁢ u ∈ V h is the Riesz representative of ⟨ u , ⋅ ⟩ ∈ V h ⋆ in the Hilbert space ( V h , ⟨ ⋅ , ⋅ ⟩ h ) .

For the time discretization, we consider a partition of the positive real axis R > 0 with constant time-step size k > 0 , i.e., t i := i ⁢ k for all i ∈ N 0 . Given a sequence { ϕ h i } i ∈ N 0 ⊂ V h , we define

(3.3) ϕ h i + 1 / 2 := ϕ h i + 1 + ϕ h i 2 and d t ⁢ ϕ h i + 1 := ϕ h i + 1 - ϕ h i k for all ⁢ i ∈ N 0 ,

as well as the piecewise linear time reconstruction

(3.4) ϕ h ⁢ k ⁢ ( t ) := t - t i k ⁢ ϕ h i + 1 + t i + 1 - t k ⁢ ϕ h i for all ⁢ i ∈ N 0 ⁢ and ⁢ t ∈ [ t i , t i + 1 ] ,

which satisfies ϕ h ⁢ k ∈ H 1 ⁢ ( Ω × ( 0 , T ) ) for any T > 0 .

3.2 Ideal Midpoint Scheme

In the following algorithm, we adapt the scheme initially proposed in [15] to the present setting. The fundamental ingredients are the midpoint rule for the time discretization, the finite element space V h endowed with the mass-lumped scalar product ⟨ ⋅ , ⋅ ⟩ h for the spatial integration, and the mapping (3.2) for the discrete realization of the effective field. We refer to the method as ideal midpoint scheme in the sense that, as we will see in the next section, practical implementations require suitable modifications.

Algorithm 3.1

Algorithm 3.1 (Ideal Midpoint Scheme)

Input:

m h 0 ∈ M h .

Loop:

for all i ∈ N 0 , compute m h i + 1 ∈ V h such that

(3.5) ⟨ d t ⁢ m h i + 1 , ϕ h ⟩ h = - ⟨ m h i + 1 / 2 × P h ⁢ h eff ⁢ ( m h i + 1 / 2 ) , ϕ h ⟩ h + α ⁢ ⟨ m h i + 1 / 2 × d t ⁢ m h i + 1 , ϕ h ⟩ h

for all ϕ h ∈ V h .

Output:

sequence of approximations { m h i } i ∈ N 0 . ∎

With the sequence of approximations { m h i } i ∈ N 0 delivered by Algorithm 3.1, we define the piecewise linear time reconstruction m h ⁢ k via (3.4). In the following theorem, we establish the stability and convergence of the approximations obtained with Algorithm 3.1. The proof is postponed to Section 5.

Theorem 3.2

The following statements hold.

Suppose that m h 0 ∈ M h . Then, for all i ∈ N 0 , (3.5) admits a solution m h i + 1 ∈ M h . In particular, the scheme preserves the unit-length constraint at any time step at the nodes of the triangulation.
The scheme is unconditionally stable in the sense that, for all J ∈ N , it holds that
(3.6) E ⁢ ( m h J ) + α ⁢ k ⁢ ∑ i = 0 J - 1 ∥ d t ⁢ m h i ∥ h 2 = E ⁢ ( m h 0 ) .
Suppose that m h 0 → m 0 in H 1 ⁢ ( Ω ) as h → 0 . Then there exist a global weak solution m : Ω × R > 0 → S 2 of (2.4) in the sense of Definition 2.1 and a subsequence of { m h ⁢ k } (not relabeled) which unconditionally converges towards 𝒎. Specifically, as h , k → 0 , m h ⁢ k ⁢ ⇀ ⋆ ⁢ m in L ∞ ⁢ ( R > 0 ; H 1 ⁢ ( Ω ; S 2 ) ) and m h ⁢ k | Ω T ⇀ m | Ω T in H 1 ⁢ ( Ω T ) for all T > 0 .

Remark 3.3

Note that, differently from the corresponding estimates for tangent plane schemes [4, 8, 30, 22], the stability result for Algorithm 3.1 (Theorem 3.2 (ii)) does not require any geometric assumption on the mesh. Moreover, (3.6) holds with equality and without any artificial dissipative term on the left-hand side.

Theorem 3.2 (i) establishes unconditional existence of a solution of (3.5), but does not provide information about its uniqueness. If k = o ⁢ ( h 2 ) , one can show that a suitable fixed-point iteration is a contraction provided that the discretization parameters are sufficiently small. With the Banach fixed-point theorem, this implies that each time step of Algorithm 3.1 is well-posed.

Proposition 3.4

Suppose that k = o ⁢ ( h 2 ) as h , k → 0 . Then there exist thresholds h 0 > 0 and k 0 > 0 such that, for all h < h 0 and k < k 0 , the variational problem (3.5) admits a unique solution m h i + 1 ∈ M h for all i ∈ N 0 .

The proof of Proposition 3.4 can be obtained simplifying the argument of the proofs of Proposition 3.5 and Theorem 3.7 below; therefore, we omit it.

3.3 Practical Midpoint Schemes

Each time step of Algorithm 3.1 requires the solution of a nonlinear system and the computation of nonlocal field contributions.

Nonlinearity is a consequence of the first term on the right-hand side of (3.5). The second term on the right-hand side, at first glance also nonlinear in m h i + 1 , turns out to be actually linear. Indeed, it holds that

m h i + 1 / 2 × d t ⁢ m h i + 1 ⁢ = ( ⁢ 3.3 ⁢ ) ⁢ m h i + 1 + m h i 2 × m h i + 1 - m h i k = - 1 k ⁢ m h i + 1 × m h i .

However, using an arbitrary off-the-shelf nonlinear solver for (3.5), the conservation and stability properties of Algorithm 3.1 established in Theorem 3.2 (i)–(ii) are in general lost. Moreover, 𝝅 can be nonlocal and non-exactly computable (e.g., for the stray field) so that the field contribution π ⁢ ( m h i + 1 / 2 ) must be numerically approximated. Hence, a direct implementation of Algorithm 3.1 should be based on an inner iteration performing the solution of the nonlinear system (3.5) and the approximate computation of π ⁢ ( m h i + 1 / 2 ) .

In the remainder of this section, we discuss and analyze an effective treatment of the nonlocal contribution, which we combine with two approaches for the linearization of (3.5), from which we will obtain two practical midpoint schemes.

To start with, we define the bilinear form a loc : H 1 ⁢ ( Ω ) × H 1 ⁢ ( Ω ) → R by

(3.7) a loc ⁢ ( ψ , ϕ ) = ∑ d = 1 3 ⟨ A d ⁢ ( ∂ d ⁡ ψ - J d ⁢ ψ ) , ∂ d ⁡ ϕ - J d ⁢ ϕ ⟩ Ω for all ⁢ ψ , ϕ ∈ H 1 ⁢ ( Ω ) .

We consider the local parts of the energy and the effective field given by

E loc ⁢ ( m ) := 1 2 ⁢ a loc ⁢ ( m , m ) - ⟨ f , m ⟩ Ω ⁢ = ( ⁢ 2.1 ⁢ ) ⁢ E ⁢ ( m ) + 1 2 ⁢ ⟨ π ⁢ ( m ) , m ⟩ Ω

and

- ⟨ h eff loc ⁢ ( m ) , ϕ ⟩ := ⟨ δ ⁢ E loc ⁢ ( m ) δ ⁢ m , ϕ ⟩ = ( ⁢ 3.7 ⁢ ) ⁢ a loc ⁢ ( m , ϕ ) - ⟨ f , ϕ ⟩ Ω = ( ⁢ 2.2 ⁢ ) ⁢ a ⁢ ( m , ϕ ) + ⟨ π ⁢ ( m ) , ϕ ⟩ Ω - ⟨ f , ϕ ⟩ Ω = ( ⁢ 2.3 ⁢ ) - ⟨ h eff ⁢ ( m ) , ϕ ⟩ + ⟨ π ⁢ ( m ) , ϕ ⟩ Ω ,

respectively. Then, for i ∈ N 0 , we rewrite (3.5) in terms of the new unknown η h i := m h i + 1 / 2 ∈ V h . Since d t ⁢ m h i + 1 = 2 ⁢ ( η h i - m h i ) / k , it is easy to see that (3.5) is equivalent to the following problem: first, compute η h i ∈ V h such that, for all ϕ h ∈ V h , it holds that

(3.8) ⟨ η h i , ϕ h ⟩ h + k 2 ⁢ ⟨ η h i × P h ⁢ h eff loc ⁢ ( η h i ) , ϕ h ⟩ h + k 2 ⁢ ⟨ η h i × P h ⁢ π ⁢ ( η h i ) , ϕ h ⟩ h + α ⁢ ⟨ η h i × m h i , ϕ h ⟩ h = ⟨ m h i , ϕ h ⟩ h .

Then define

(3.9) m h i + 1 := 2 ⁢ η h i - m h i .

To treat the nonlocal contribution π ⁢ ( η h i ) , we adopt the implicit-explicit (IMEX) approach introduced in [39]. Let π h : V h → V h be an operator approximating 𝝅, assumed to be linear and uniformly bounded in L 2 ⁢ ( Ω ) in the sense that ∥ π h ∥ L ⁢ ( L 2 ⁢ ( Ω ) ; L 2 ⁢ ( Ω ) ) ≤ C π for some C π > 0 independent of ℎ. Moreover, we say that π h is consistent with 𝝅 if, for all ϕ ∈ L 2 ⁢ ( Ω ) and all ( ϕ h ) h > 0 ⊂ V h with ϕ h → ϕ in L 2 ⁢ ( Ω ) as h → 0 , it holds that

(3.10) π h ⁢ ( ϕ h ) → π ⁢ ( ϕ ) in ⁢ L 2 ⁢ ( Ω ) as ⁢ h → 0 .

We define m h - 1 := m h 0 and

Π h ⁢ ( m h i , m h i - 1 ) := 3 2 ⁢ π h ⁢ ( m h i ) - 1 2 ⁢ π h ⁢ ( m h i - 1 ) .

Then, in (3.8), we replace π ⁢ ( η h i ) with its approximation Π h ⁢ ( m h i , m h i - 1 ) to obtain

(3.11) ⟨ η h i , ϕ h ⟩ h + k 2 ⁢ ⟨ η h i × P h ⁢ h eff loc ⁢ ( η h i ) , ϕ h ⟩ h + k 2 ⁢ ⟨ η h i × P h ⁢ Π h ⁢ ( m h i , m h i - 1 ) , ϕ h ⟩ h + α ⁢ ⟨ η h i × m h i , ϕ h ⟩ h = ⟨ m h i , ϕ h ⟩ h .

In particular, the nonlocal contribution, treated explicitly, becomes independent of the unknown η h i . We now discuss two strategies to linearize (3.11) in order to arrive at two practical midpoint schemes. To emphasize the inexact solution of (3.11) up to some accuracy ε > 0 , we write m h ⁢ ε i rather than m h i for the iterates of the practical (linearized) midpoint schemes.

3.3.1 Constraint-Preserving Fixed-Point Iteration

We solve (3.11) with the following fixed-point iteration. Let ε > 0 denote some prescribed tolerance. Set η h i , 0 := m h ⁢ ε i . For ℓ ∈ N 0 , given η h i , ℓ ∈ V h , compute η h i , ℓ + 1 ∈ V h such that, for all ϕ h ∈ V h , it holds that

(3.12) ⟨ η h i , ℓ + 1 , ϕ h ⟩ h + k 2 ⁢ ⟨ η h i , ℓ + 1 × P h ⁢ h eff loc ⁢ ( η h i , ℓ ) , ϕ h ⟩ h + k 2 ⁢ ⟨ η h i , ℓ + 1 × P h ⁢ Π h ⁢ ( m h ⁢ ε i , m h ⁢ ε i - 1 ) , ϕ h ⟩ h + α ⁢ ⟨ η h i , ℓ + 1 × m h ⁢ ε i , ϕ h ⟩ h = ⟨ m h ⁢ ε i , ϕ h ⟩ h ,

until

(3.13) ∥ I h ⁢ [ η h i , ℓ + 1 × P h ⁢ ( h eff loc ⁢ ( η h i , ℓ + 1 ) - h eff loc ⁢ ( η h i , ℓ ) ) ] ∥ h ≤ ε ,

where I h ⁢ [ ⋅ ] denotes the vector-valued nodal interpolant. If ℓ * ∈ N 0 is the smallest integer for which the stopping criterion (3.13) is satisfied, in view of (3.9), the approximate magnetization at the new time step is defined as m h ⁢ ε i + 1 := 2 ⁢ η h i , ℓ * + 1 - m h ⁢ ε i .

In the following proposition, we collect the properties of the proposed fixed-point iteration. The proof is postponed to Section 6.

Proposition 3.5

Let i ∈ N 0 .

For all ℓ ∈ N 0 , the variational problem (3.12) admits a unique solution η h i , ℓ + 1 ∈ V h . Moreover, it holds that ∥ η h i , ℓ + 1 ∥ L ∞ ⁢ ( Ω ) ≤ 1 .
If k = o ⁢ ( h 2 ) as h , k → 0 , there exist a contraction constant 0 < q < 1 and thresholds h 0 , k 0 > 0 such that, for all h < h 0 and k < k 0 , it holds that
(3.14) ∥ η h i , ℓ + 2 - η h i , ℓ + 1 ∥ h ≤ q ⁢ ∥ η h i , ℓ + 1 - η h i , ℓ ∥ h for all ⁢ ℓ ∈ N 0 .
The constants q , h 0 , k 0 depend only on the mesh parameter 𝜅 and the problem data.
Under the assumptions of part (ii), the stopping criterion (3.13) is met in a finite number of iterations. If ℓ * ∈ N 0 denotes the smallest integer for which (3.13) is satisfied, the new approximation
m h ⁢ ε i + 1 := 2 ⁢ η h i , ℓ * + 1 - m h ⁢ ε i ∈ V h
belongs to M h .

For all i ∈ N 0 , let r h ⁢ ε i := P h ⁢ ( h eff loc ⁢ ( η h i , ℓ * + 1 ) - h eff loc ⁢ ( η h i , ℓ * ) ) ∈ V h . Because of the stopping criterion (3.13), it holds that ∥ I h ⁢ [ m h ⁢ ε i + 1 / 2 × r h ⁢ ε i ] ∥ h ≤ ε . With this definition, the proposed linearization of Algorithm 3.1 is covered by the following algorithm.

Algorithm 3.6

Algorithm 3.6 (Practical Midpoint Scheme, Constraint-Preserving Fixed-Point Iteration)

Input:

m h ⁢ ε 0 := m h 0 ∈ M h .

Loop:

for all i ∈ N 0 , use the constraint-preserving fixed-point iteration (3.12)–(3.13) to compute m h ⁢ ε i + 1 ∈ M h and r h ⁢ ε i ∈ V h with ∥ I h ⁢ [ m h ⁢ ε i + 1 / 2 × r h ⁢ ε i ] ∥ h ≤ ε such that

(3.15) ⟨ d t ⁢ m h ⁢ ε i + 1 , ϕ h ⟩ h = - ⟨ m h ⁢ ε i + 1 / 2 × P h ⁢ h eff ⁢ ( m h ⁢ ε i + 1 / 2 ) , ϕ h ⟩ h + α ⁢ ⟨ m h ⁢ ε i + 1 / 2 × d t ⁢ m h ⁢ ε i + 1 , ϕ h ⟩ h + ⟨ m h ⁢ ε i + 1 / 2 × [ r h ⁢ ε i + P h ⁢ ( π ⁢ ( m h ⁢ ε i + 1 / 2 ) - Π h ⁢ ( m h ⁢ ε i , m h ⁢ ε i - 1 ) ) ] , ϕ h ⟩ h

for all ϕ h ∈ V h .

Output:

sequence of approximations { m h ⁢ ε i } i ∈ N 0 .∎

In the following theorem, we establish the stability and convergence of the approximations obtained with Algorithm 3.6. The proof is postponed to Section 6.

Theorem 3.7

The following statements hold.

Suppose that m h 0 ∈ M h . If k = o ⁢ ( h 2 ) as h , k → 0 , there exist thresholds h 0 > 0 and k 0 > 0 such that, for all h < h 0 and k < k 0 , (3.15) admits solutions m h ⁢ ε i + 1 ∈ M h and r h ⁢ ε i ∈ V h with ∥ I h ⁢ [ m h ⁢ ε i + 1 / 2 × r h ⁢ ε i ] ∥ h ≤ ε for all i ∈ N 0 . In particular, the scheme preserves the unit-length constraint at the nodes of the triangulation for all time steps. The thresholds h 0 , k 0 depend only on the mesh parameter 𝜅 and the problem data.
Under the assumptions of part (i), for all h < h 0 , k < k 0 , and J ∈ N , the scheme satisfies the discrete energy identity
(3.16) E ( m h ⁢ ε J ) + α k ∑ i = 0 J - 1 ∥ d t m h ⁢ ε i + 1 ∥ h 2 = E ( m h 0 ) - k ∑ i = 0 J - 1 ⟨ m h ⁢ ε i + 1 / 2 × [ r h ⁢ ε i + P h ⁢ ( π ⁢ ( m h ⁢ ε i + 1 / 2 ) - Π h ⁢ ( m h ⁢ ε i , m h ⁢ ε i - 1 ) ) ] , P h h eff ( m h ⁢ ε i + 1 / 2 ) - α d t m h ⁢ ε i + 1 ⟩ h .
Let T > 0 . Let J ∈ N be the smallest integer such that T ≤ k ⁢ J . Under the assumptions of part (i), if ε = O ⁢ ( h ) and { m h 0 } h > 0 is bounded in H 1 ⁢ ( Ω ) as h , ε → 0 , there exist thresholds 0 < h 0 * ≤ h 0 , 0 < k 0 * ≤ k 0 , and ε 0 * > 0 such that, for all h < h 0 * , k < k 0 * , ε < ε 0 * , and 1 ≤ j ≤ J , we have the stability estimate
(3.17) ∥ m h ⁢ ε j ∥ H 1 ⁢ ( Ω ) 2 + k ⁢ ∑ i = 0 j - 1 ∥ d t ⁢ m h ⁢ ε i + 1 ∥ L 2 ⁢ ( Ω ) 2 ≤ C .
The constant C > 0 and the thresholds h 0 * , k 0 * , ε 0 * depend only on the mesh parameter 𝜅, the final time 𝑇, and the problem data.
In addition to the assumptions of part (iii), assume m h 0 → m 0 in H 1 ⁢ ( Ω ) as h → 0 , and suppose that π h is consistent (3.10) with 𝝅. Then there exist m ∈ H 1 ⁢ ( Ω T ) ∩ L ∞ ⁢ ( 0 , T ; H 1 ⁢ ( Ω ; S 2 ) ) and a subsequence of { m h ⁢ ε ⁢ k } (not relabeled) which converges towards 𝒎 as h , k , ε → 0 . Specifically,
m h ⁢ ε ⁢ k | Ω T ⁢ ⇀ ⋆ ⁢ m in ⁢ L ∞ ⁢ ( 0 , T ; H 1 ⁢ ( Ω ; S 2 ) ) and m h ⁢ ε ⁢ k | Ω T ⇀ m in ⁢ H 1 ⁢ ( Ω T ) as ⁢ h , k , ε → 0 .
The limit function 𝒎 satisfies conditions (i)–(iv) of Definition 2.1.

3.3.2 Newton Iteration

Based on the Newton scheme, in [11, Section 1.4.1], the authors employ a linearization of the nonlinear system (3.5) in the ideal midpoint scheme with simplified effective field, i.e., without nonlocal contributions and without DMI. Their 2D numerical experiments give hope for a less restrictive CFL condition than for the fixed-point iteration from Section 3.3.1.

For three-dimensional micromagnetics and considering the full effective field (2.3), in Section 7.3, we apply Newton’s method to the nonlinear system of equations (3.11) resulting in the following iteration. Let ε > 0 denote some tolerance. Set η h i , 0 := m h ⁢ ε i . For ℓ ∈ N 0 , given η h i , ℓ ∈ V h , compute u h i , ℓ ∈ V h such that, for all ϕ h ∈ V h , it holds that

(3.18) ⟨ u h i , ℓ , ϕ h ⟩ h + k 2 ⁢ ⟨ u h i , ℓ × P h ⁢ ( h eff loc ⁢ ( η h i , ℓ ) + Π h ⁢ ( m h ⁢ ε i , m h ⁢ ε i - 1 ) ) , ϕ h ⟩ h + α ⁢ ⟨ u h i , ℓ × m h ⁢ ε i , ϕ h ⟩ h + k 2 ⁢ ⟨ η h i , ℓ × P h ⁢ ( h eff loc ⁢ ( u h i , ℓ ) - f ) , ϕ h ⟩ h = ⟨ m h ⁢ ε i - η h i , ℓ , ϕ h ⟩ h - k 2 ⁢ ⟨ η h i , ℓ × P h ⁢ ( h eff loc ⁢ ( η h i , ℓ ) + Π h ⁢ ( m h ⁢ ε i , m h ⁢ ε i - 1 ) ) , ϕ h ⟩ h - α ⁢ ⟨ η h i , ℓ × m h ⁢ ε i , ϕ h ⟩ h ,

and define η h i , ℓ + 1 := η h i , ℓ + u h i , ℓ until

(3.19) ∥ I h ⁢ [ u h i , ℓ × P h ⁢ ( h eff loc ⁢ ( u h i , ℓ ) - f ) ] ∥ h ≤ ε .

If ℓ * ∈ N 0 is the smallest integer for which the stopping criterion (3.19) is satisfied, the approximate magnetization at the new time step is defined as m h ⁢ ε i + 1 := 2 ⁢ η h i , ℓ * + 1 - m h ⁢ ε i .

For all i ∈ N 0 , let r h ⁢ ε i = I h ⁢ [ u h i , ℓ * × P h ⁢ ( h eff loc ⁢ ( u h i , ℓ * ) - f ) ] ∈ V h . In view of the stopping criterion (3.19), it holds that ∥ r h ⁢ ε i ∥ h ≤ ε . With this definition, the proposed linearization of Algorithm 3.1 based on the Newton method is covered by the following algorithm.

Algorithm 3.8

Algorithm 3.8 (Practical Midpoint Scheme, Newton Iteration)

Input:

m h ⁢ ε 0 := m h 0 ∈ V h .

Loop:

for all i ∈ N 0 , use Newton’s method (3.18)–(3.19) with initial guess η h i , 0 := m h ⁢ ε i to compute m h ⁢ ε i + 1 ∈ V h and r h ⁢ ε i ∈ V h with ∥ r h ⁢ ε i ∥ h ≤ ε such that

(3.20) ⟨ d t ⁢ m h ⁢ ε i + 1 , ϕ h ⟩ h = - ⟨ m h ⁢ ε i + 1 / 2 × P h ⁢ h eff ⁢ ( m h ⁢ ε i + 1 / 2 ) , ϕ h ⟩ h + α ⁢ ⟨ m h ⁢ ε i + 1 / 2 × d t ⁢ m h ⁢ ε i + 1 , ϕ h ⟩ h + ⟨ r h ⁢ ε i , ϕ h ⟩ h + ⟨ m h ⁢ ε i + 1 / 2 × P h ⁢ ( π ⁢ ( m h ⁢ ε i + 1 / 2 ) - Π h ⁢ ( m h ⁢ ε i , m h ⁢ ε i - 1 ) ) , ϕ h ⟩ h

for all ϕ h ∈ V h .

Output:

sequence of approximations { m h ⁢ ε i } i ∈ N 0 .∎

The results on Algorithm 3.8 are stated in Lemma 3.9 ( L ∞ -bound), Theorem 3.10 (stability), and Theorem 3.11 (well-posedness) below. Compared to Algorithm 3.6, our analysis is more involved: precisely, for i < J , the proof of Theorem 3.11 requires 𝑖-independent bounds on ∥ m h ⁢ ε i ∥ L ∞ ⁢ ( Ω ) and E ⁢ ( m h ⁢ ε i ) in order to guarantee well-posedness of Algorithm 3.8, while Lemma 3.9 and Theorem 3.10 require termination of (3.18)–(3.19) so that m h ⁢ ε i + 1 is well-defined.

In contrast to the fixed-point iteration (3.12)–(3.13), the Newton iteration (3.18)–(3.19) does not inherently preserve discrete unit length, i.e., in general, | m h ⁢ ε i + 1 ⁢ ( z ) | ≠ | m h ⁢ ε i ⁢ ( z ) | for z ∈ N h . However, assuming well-posedness of Algorithm 3.8, in the following lemma, we establish uniform L ∞ ⁢ ( Ω ) -boundedness of the approximations obtained with Algorithm 3.8.

Lemma 3.9

Suppose h , k , ε > 0 , m h 0 ∈ M h , and let J ∈ N be the smallest integer such that T ≤ k ⁢ J . Let 0 ≤ i < J , and suppose that the Newton iteration (3.18)–(3.19) in Algorithm 3.8 terminates for all 0 ≤ n ≤ i , i.e., the sequences { m h ⁢ ε n } n = 0 i + 1 , { r h ⁢ ε n } n = 0 i ⊂ V h are the output of Algorithm 3.8 and satisfy (3.20) with ∥ r h ⁢ ε n ∥ h ≤ ε for all 0 ≤ n ≤ i .

If ε = O ⁢ ( h 3 / 2 ) as h , k , ε → 0 , then there exist a constant C ∞ > 0 and thresholds h 0 > 0 , k 0 > 0 , and ε 0 > 0 such that, for all h < h 0 , k < k 0 , and ε < ε 0 , it holds that ∥ m h ⁢ ε n + 1 ∥ L ∞ ⁢ ( Ω ) ≤ C ∞ uniformly for all 0 ≤ n ≤ i . The thresholds h 0 , k 0 , ε 0 depend only on the mesh parameter 𝜅 and the problem data, while the bound C ∞ > 0 depends only on 𝜅, ε ⁢ h - 3 / 2 ≲ 1 , and the final time T > 0 , but not on the integer i < J .
If ε = o ⁢ ( h 3 / 2 ) , then there holds ∥ I h ⁢ ( | m h ⁢ ε ⁢ k | 2 ) - 1 ∥ L ∞ ⁢ ( [ 0 , t i + 1 ] × Ω ) → 0 as h , k , ε → 0 .

Assuming well-posedness of Algorithm 3.8, in the following theorem, we establish the stability and convergence of the approximations obtained with Algorithm 3.8. The proof is postponed to Section 7.1.

Theorem 3.10

Let T > 0 , and suppose that { T h } h > 0 is a 𝜅-quasi-uniform family of triangulations. Suppose h , k , ε > 0 , m h 0 ∈ M h , and let J ∈ N be the smallest integer such that T ≤ k ⁢ J . Let 0 ≤ i < J , and suppose that the Newton iteration (3.18)–(3.19) in Algorithm 3.8 terminates for all 0 ≤ n ≤ i , i.e., the sequences { m h ⁢ ε n } n = 0 i + 1 , { r h ⁢ ε n } n = 0 i ⊂ V h are the output of Algorithm 3.8 and satisfy (3.20) with ∥ r h ⁢ ε n ∥ h ≤ ε for all 0 ≤ n ≤ i .

Under these assumptions, the scheme satisfies the discrete energy identity
E ( m h ⁢ ε i + 1 ) + α k ∑ n = 0 i ∥ d t m h ⁢ ε n + 1 ∥ h 2 = E ( m h 0 ) - k ∑ n = 0 i ⟨ r h ⁢ ε n + m h ⁢ ε n + 1 / 2 × P h ⁢ ( π ⁢ ( m h ⁢ ε n + 1 / 2 ) - Π h ⁢ ( m h ⁢ ε n , m h ⁢ ε n - 1 ) ) , P h h eff ( m h ⁢ ε n + 1 / 2 ) - α d t m h ⁢ ε n + 1