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.
1.1 Energetics of a Ferromagnet
In the continuum theory of micromagnetism, whose origin dates back to the seminal work of Landau–Lifshitz  on small ferromagnetic particles, the amount of magnetic moment (per unit volume) of a rigid ferromagnetic body occupying a bounded region is represented by a classical vector field, the magnetization . Its module describes the so-called saturation magnetization. In single-crystal ferromagnets [2, 6], 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 , where is a vector field with values in the unit sphere of , and the observable magnetization states minimize the micromagnetic energy functional [16, 31]
defined for every .
The exchange energy penalizes spatial variations of the direction of the magnetization, with representing a material-dependent constant that summarizes the stiffness of short-range (symmetric) exchange interactions. The second term 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 affects the chirality of the ferromagnetic system [43, 42]. The third term is the magnetostatic self-energy, i.e., the energy due to the stray field induced by . From the mathematical point of view, can be characterized as the projection of on the closed subspace of gradient vector fields (see, e.g., [38, 21]). Here, denotes the vacuum permeability. Additionally, the micromagnetic energy includes two additional energy contributions: the magnetocrystalline anisotropy energy and the Zeeman energy . The energy density models the existence of easy directions of the magnetization due to the crystallographic structure of the ferromagnet, while models the tendency of a specimen to have the magnetization aligned with the external applied field , 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 , 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 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 . The easiest way is to assume a strong coupling condition [6, 5, 20]: although 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 and . 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 , , and 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
where denotes the standard basis of . It is therefore a special case of the energy density
with being 3-by-3 antisymmetric matrices, i.e., . Similarly, the symmetric exchange energy density can be generalized to the density
with being 3-by-3 invertible symmetric matrices, i.e., . Hence, for , it holds that
Note that 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
This phenomenological equation describes the magnetization dynamics as a dissipative precession driven by the effective field , and modulated by the gyromagnetic ratio of the electron and the Gilbert damping parameter . 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 . 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  proves unconditional convergence of the finite element approximation towards a weak solution of LLG in the sense of  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  focuses on the design and the analysis of effective approaches to treat the nonlocal field 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  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.
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.
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 and by . We denote by the scalar product in and by the duality pairing between and its dual. By , 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 be a bounded Lipschitz domain. The energy of is given by
where , while the bilinear form is defined, for all , by
Here, is a linear, bounded, and self-adjoint operator, while, for , the 3-by-3 matrices and have coefficients in , with being also symmetric and uniformly positive definite, i.e., it holds that and
The existence of minimizers of (2.1) in follows from the direct method of calculus of variations. Moreover, any minimizer satisfies the Euler–Lagrange equations
Here, is the (negative) Gâteaux derivative of the energy, i.e.,
Turning to the dynamical case, a non-equilibrium configuration evolves according to (1.3), which, after a suitable rescaling, reads
with being the Gilbert damping parameter. The dynamics is dissipative in the sense that any sufficiently smooth solution of (2.4) satisfies the energy law
Let . A vector field is called a global weak solution of (2.4) if and, for all , with , the following properties are satisfied:
in the sense of traces;
for all , it holds that(2.6)
it holds that(2.7)
We note that (2.3) implicitly includes natural boundary conditions on 𝒎, which are homogeneous Neumann boundary conditions if and for . For a more explicit presentation, we refer to . 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).
respectively, where .
(ii) The present setting covers and generalizes the model problems considered in previous mathematical works on the numerical integration of LLG.
With the choices and for , where is a characteristic length associated with DMI, , and , we obtain that
which is the setting analyzed in  for the simulation of chiral magnetic skyrmions by the means of a family of tangent plane integrators. Here, denotes the matrix such that for all .
3 Numerical Algorithms and Main Results
Let . For the spatial discretization, assuming Ω to be a polyhedral domain, we consider a 𝜅-quasi-uniform family of regular tetrahedral triangulations of Ω parametrized by the mesh size , i.e., for all . We denote by the set of vertices of . For any , we denote by the space of first-order polynomials on 𝐾. We consider the space of -piecewise affine and globally continuous finite elements
The classical basis for this finite-dimensional subspace of is the set of nodal hat functions , which satisfy for all . The nodal interpolant is defined by for all .
Let . For each time step, approximate solutions of (2.4) are sought in the set of admissible approximate magnetizations
which consists of all elements of satisfying the unit-length constraint at the nodes of the triangulation.
Besides the standard scalar product , given a mesh and the associated nodal interpolant , we consider the mass-lumped product defined by
Using the definition of the nodal interpolant, we see that
On , is a scalar product and the induced norm is equivalent to the standard norm of . In particular, it holds that
see [13, Lemma 3.9]. Finally, we define the mapping by
i.e., is the Riesz representative of in the Hilbert space .
For the time discretization, we consider a partition of the positive real axis with constant time-step size , i.e., for all . Given a sequence , we define
as well as the piecewise linear time reconstruction
which satisfies for any .
3.2 Ideal Midpoint Scheme
In the following algorithm, we adapt the scheme initially proposed in  to the present setting. The fundamental ingredients are the midpoint rule for the time discretization, the finite element space endowed with the mass-lumped scalar product 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 (Ideal Midpoint Scheme)
for all , compute such that(3.5)
for all .
sequence of approximations . ∎
With the sequence of approximations delivered by Algorithm 3.1, we define the piecewise linear time reconstruction 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.
The following statements hold.
Suppose that . Then, for all , (3.5) admits a solution . 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 , it holds that(3.6)
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 , 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.
Suppose that as . Then there exist thresholds and such that, for all and , the variational problem (3.5) admits a unique solution for all .
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 , turns out to be actually linear. Indeed, it holds that
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 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 .
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 by
We consider the local parts of the energy and the effective field given by
To treat the nonlocal contribution , we adopt the implicit-explicit (IMEX) approach introduced in . Let be an operator approximating 𝝅, assumed to be linear and uniformly bounded in in the sense that for some independent of ℎ. Moreover, we say that is consistent with 𝝅 if, for all and all with in as , it holds that
We define and
Then, in (3.8), we replace with its approximation to obtain
In particular, the nonlocal contribution, treated explicitly, becomes independent of the unknown . 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 , we write rather than 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 denote some prescribed tolerance. Set . For , given , compute such that, for all , it holds that
where denotes the vector-valued nodal interpolant. If 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 .
In the following proposition, we collect the properties of the proposed fixed-point iteration. The proof is postponed to Section 6.
For all , the variational problem (3.12) admits a unique solution . Moreover, it holds that .
If as , there exist a contraction constant and thresholds such that, for all and , it holds that(3.14)
The constants depend only on the mesh parameter 𝜅 and the problem data.
Algorithm 3.6 (Practical Midpoint Scheme, Constraint-Preserving Fixed-Point Iteration)
The following statements hold.
Suppose that . If as , there exist thresholds and such that, for all and , (3.15) admits solutions and with for all . In particular, the scheme preserves the unit-length constraint at the nodes of the triangulation for all time steps. The thresholds depend only on the mesh parameter 𝜅 and the problem data.
Under the assumptions of part (i), for all , , and , the scheme satisfies the discrete energy identity(3.16)
Let . Let be the smallest integer such that . Under the assumptions of part (i), if and is bounded in as , there exist thresholds , , and such that, for all , , , and , we have the stability estimate(3.17)
The constant and the thresholds depend only on the mesh parameter 𝜅, the final time 𝑇, and the problem data.
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 denote some tolerance. Set . For , given , compute such that, for all , it holds that
and define until
If is the smallest integer for which the stopping criterion (3.19) is satisfied, the approximate magnetization at the new time step is defined as .
For all , let . In view of the stopping criterion (3.19), it holds that . With this definition, the proposed linearization of Algorithm 3.1 based on the Newton method is covered by the following algorithm.
Algorithm 3.8 (Practical Midpoint Scheme, Newton Iteration)
The results on Algorithm 3.8 are stated in Lemma 3.9 ( -bound), Theorem 3.10 (stability), and Theorem 3.11 (well-posedness) below. Compared to Algorithm 3.6, our analysis is more involved: precisely, for , the proof of Theorem 3.11 requires 𝑖-independent bounds on and 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 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, for . However, assuming well-posedness of Algorithm 3.8, in the following lemma, we establish uniform -boundedness of the approximations obtained with Algorithm 3.8.
Suppose , , and let be the smallest integer such that . Let , and suppose that the Newton iteration (3.18)–(3.19) in Algorithm 3.8 terminates for all , i.e., the sequences are the output of Algorithm 3.8 and satisfy (3.20) with for all .
If as , then there exist a constant and thresholds , , and such that, for all , , and , it holds that uniformly for all . The thresholds depend only on the mesh parameter 𝜅 and the problem data, while the bound depends only on 𝜅, , and the final time , but not on the integer .
If , then there holds as .
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.
Let , and suppose that is a 𝜅-quasi-uniform family of triangulations. Suppose , , and let be the smallest integer such that . Let , and suppose that the Newton iteration (3.18)–(3.19) in Algorithm 3.8 terminates for all , i.e., the sequences are the output of Algorithm 3.8 and satisfy (3.20) with for all .
Under these assumptions, the scheme satisfies the discrete energy identity