Convergent tangent plane integrators for the simulation of chiral magnetic skyrmion dynamics

We consider the numerical approximation of the Landau-Lifshitz-Gilbert equation, which describes the dynamics of the magnetization in ferromagnetic materials. In addition to the classical micromagnetic contributions, the energy comprises the Dzyaloshinskii-Moriya interaction, which is the most important ingredient for the enucleation and the stabilization of chiral magnetic skyrmions. We propose and analyze three tangent plane integrators, for which we prove (unconditional) convergence of the finite element solutions towards a weak solution of the problem. The analysis is constructive and also establishes existence of weak solutions. Numerical experiments demonstrate the applicability of the methods for the simulation of practically relevant problem sizes.

experimentally observed [47,54] in several magnetic systems. The most important ingredient for the enucleation and the stabilization of magnetic skyrmions is the socalled Dzyaloshinskii-Moriya interaction (DMI) (see [26,46]). It is a short-range effect, sometimes also referred to as antisymmetric exchange, which exerts a torque on the magnetization inducing neighboring spins to be perpendicular to each other. It is thus in direct competition with the classical Heisenberg exchange interaction, which conversely favors uniform configurations. The DMI is modeled by an energy contribution, which is linear in the first spatial derivatives of the magnetization and is added to the micromagnetic energy for chiral ferromagnets. Magnetic skyrmions are currently subject of intense scientific research, which includes theoretical, computational, and experimental studies (see, e.g., [16,34,35,40,57]). As for the mathematical literature, the existence of isolated skyrmions emerging as energy minimizers of two-dimensional micromagnetic models and their dynamic stability have been investigated in [25,45], whereas chiral domain walls in ultrathin ferromagnetic films have been studied in [48]. The growing interest in skyrmions in the magnetic storage and magnetic logic community is connected with their potential as possible candidate to store the bits of future devices, with the information being encoded as the presence/absence of a skyrmion (see, e.g., [29,60]) for the proposal of skyrmion racetrack memories, which are believed to overcome the original domain wall-based device of [50] and pave new ways in magnetic data logic [36].
A well-accepted model for the magnetization dynamics is the Landau-Lifshitz-Gilbert equation (LLG) [33,42]. The numerical approximation of LLG poses several challenges: nonlinearities, a nonconvex pointwise constraint, an intrinsic energy law, which resembles the one of a gradient flow and combines conservative and dissipative effects, and the possible coupling with other partial differential equations (PDEs), e.g., the Maxwell equations.
The numerical integration of LLG has been the subject of several mathematical studies (see, e.g., [32,41,53]). A well-established approach is represented by the integrators usually referred to as tangent plane schemes. These methods are based on equivalent reformulations of the equation in the tangent space.
The integrator proposed in [5], which considers the case in which the energy only comprises the exchange contribution, requires only the solution of one linear system per time-step, is formally of first order in time, and is unconditionally convergent towards a weak solution of the problem, i.e., the numerical analysis of the scheme does not require to impose any restrictive CFL-type coupling condition on the timestep size and the spatial mesh size. The pointwise constraint is enforced by applying the nodal projection to the computed solution at each time-step. The scheme generalizes the explicit scheme proposed in [6] and analyzed in [14]. Implicit-explicit approaches of the algorithm of [5] for the full effective field were independently introduced and analyzed in [8,22]. Extensions of the scheme for the discretization of the coupling of LLG with other PDEs were studied in [10,11,43,44]. Inspired by [13], the projection-free version of the algorithm of [5], which avoids the use of the nodal projection, was introduced, analyzed, and applied to the decoupled integration of the coupling of LLG with a spin diffusion equation for the spin accumulation in [3]. The violation of the constraint at the nodes of the mesh occurring in this case is uniformly controlled by the time-step size. The projection-free tangent plane scheme of [3] was combined with a FEM-BEM coupling method for the discretization of the coupling of LLG with the magnetoquasistatic Maxwell equations in full space in [28]. There, assuming the existence of a unique sufficiently smooth solution, the authors proved optimal first-order convergence rates of the method. A tangent plane scheme characterized by an enhanced convergence order in time was proposed in [7]. The method is unconditionally convergent and formally of (almost) second order in time. A more efficient implicit-explicit version of this method has been proposed in [24]. Adapting ideas from [5,15], the recent work [39] proposes a similar predictor-corrector scheme based on a linear mass-lumped variational formulation of LLG.

Contributions and general outline of the present work
In this work, as a novel contribution, we introduce and analyze three tangent plane schemes for LLG in the presence of DMI. The integrators extend to this case the firstorder scheme of [5] (Algorithm 1), its projection-free variant from [3] (Algorithm 2), and the (almost) second-order scheme of [7] (Algorithm 3). For any algorithm, we prove that the sequence of finite element solutions, upon extraction of a subsequence, converges towards a weak solution of the problem. For the projection-free algorithm, we prove that the convergence is even unconditional, while the stability analysis requires a mild CFL-type condition on the discretization parameters and a geometric restriction on the underlying mesh for the other two approaches. The present extension of the LLG analysis is not straightforward, since the DMI term involves magnetization derivatives, is neither self-adjoint nor positive definite, and requires to impose different boundary conditions on LLG, which entail a careful treatment. A by-product of our constructive analysis is the proof of existence of weak solutions, which to our knowledge was missing in the literature. Finally, numerical experiments show that our approach can be used to study enucleation processes, stability, and dynamics of magnetic skyrmions.
The remainder of the work is organized as follows: For the convenience of the reader, we conclude this section by collecting the notation used throughout the paper. In Section 2, we propose an organic presentation of the physical background and the mathematical framework of the problem under consideration. In Section 3, we derive three tangent plane schemes and state the convergence result (Theorem 1). Section 4 is devoted to numerical experiments. Finally, in Section 5, we present the convergence analysis of the algorithms and, in particular, we establish the proof of Theorem 1.

Notation
We use the standard notation for Lebesgue, Sobolev, and Bochner spaces and norms (see, e.g., [27,Chapter 5] or [18,Chapter 2]). In the case of (spaces of) vector-valued or matrix-valued functions, we use bold letters, e.g., for any domain U , we denote both L 2 (U ; R 3 ) and L 2 (U ; R 3×3 ) by L 2 (U ). For the differential operators, we use the following notation: For a scalar function f , we denote by ∇f the gradient and by f the Laplace operator. For a vector-valued function f , we denote by ∇ · f the divergence, by ∇ × f the curl, by ∇f the Jacobian, and by f the vector-valued Laplace operator. Given another vector-valued function h, we also define (f · ∇)h by [(f · ∇)h] i = f · ∇h i for all 1 ≤ i ≤ 3. We denote the unit sphere by S 2 = {x ∈ R 3 : |x| = 1} and by {e i } 1≤i≤3 ∈ R 3 the standard basis of R 3 , i.e., (e i ) j = δ ij for all 1 ≤ i, j ≤ 3. Given a vector b ∈ R 3 and a matrix A ∈ R 3×3 (with columns a i ∈ R 3 for all 1 ≤ i ≤ 3), we denote by A × b ∈ R 3×3 the matrix whose columns are a i × b for all 1 ≤ i ≤ 3. By C > 0, we always denote a generic constant, which is independent of the discretization parameters, but not necessarily the same at each occurrence. We also use the notation to denote smaller than or equal to up to a multiplicative constant, i.e., we write A B if there exists a constant C > 0, which is clear from the context and always independent of the discretization parameters, such that A ≤ CB.

Let
⊂ R 3 be a bounded domain with boundary := ∂ . The dynamics of the normalized magnetization m = (m 1 , m 2 , m 3 ) ∈ S 2 is governed by LLG, which in the so-called Gilbert form reads (see [33,42]). Here, γ 0 ≈ 2.21e5 m/(As) is the rescaled gyromagnetic ratio of the electron, 0 < α ≤ 1 is the dimensionless Gilbert damping parameter, and H eff is the energy-based effective field (in A/m), i.e., it holds that where E (·) is the total energy, μ 0 = 4πe−7 N/A 2 is the vacuum permeability, and M s > 0 is the saturation magnetization (in A/m). In micromagnetics, the total energy is usually the sum of the following standard terms: -the Heisenberg exchange contribution where A > 0 denotes the exchange stiffness constant (in J/m); -the magnetocrystalline anisotropy contribution, which for the uniaxial case reads where K > 0 denotes the anisotropy constant (in J/m 3 ) and a ∈ S 2 is the easy axis, -the Zeeman contribution where H ext denotes an applied external field (in A/m); -the magnetostatic contribution where H s (m) = −∇u denotes the stray field (in A/m), with u being the magnetostatic potential (in A), which solves the full-space transmission problem Here, n : → S 2 denotes the outward-pointing unit normal vector to .
In this work, the total energy E (·) in (2) also comprises a term associated with the DMI [26,46]. This energy contribution is phenomenologically introduced for systems with a broken symmetry due to different interface crystal configurations as a linear combination of the so-called Lifshitz invariants, i.e., the components of the chirality tensor A = (a ij ) 1≤i,j ≤3 = ∇m × m (see [19,20]). The choice of the appropriate DMI form depends on the crystal structure of the material and on the geometry of the sample under consideration. This choice in turn determines the specific expression of the effective field (according to (2)) and the boundary conditions on , which are chosen in agreement with those satisfied by the solution of the Euler-Lagrange equations associated with the energy minimization problem In micromagnetics, two main DMI forms are usually considered: -For helimagnetic materials [46], the DMI is obtained by taking as Lifshitz invariants the trace of the matrix A, i.e., The so-called bulk DMI energy contribution then takes the form so that the resulting effective field term and boundary conditions on are given by -Another type of DMI is due to the interfaces between different materials which break inversion symmetry [23]. For a magnetic thin film aligned with the x 1 x 2plane, the Lifshitz invariants are given by The so-called interfacial DMI energy contribution thus takes the form (6) so that the resulting effective field term and boundary conditions on are given by In (5)-(6), the constant D ∈ R is the DMI constant (in J/m 2 ). 1 The sign of D determines the chirality of the system, which, in the case of a skyrmion state, defines the sense of rotation of the magnetization along the skyrmion diameter [40,57].

Problem formulation
To simplify the notation, in our analysis, we restrict ourselves to the case in which the energy solely comprises exchange and DMI. We refer to [3,7,8,22,24] for the design and the analysis of effective tangent plane integrators for the standard energy terms (exchange, anisotropy, Zeeman, and magnetostatic). Moreover, without loss of generality, we assume that D > 0 and consider the prototypical case of bulk DMI, because it is characterized by a short notation involving the curl operator. The same approach and the same results hold for all possible choices of the Lifshitz invariants.
After a suitable nondimensionalization, 2 the initial boundary value problem in which we are interested takes the form where the effective field h eff (m) is determined by the energy functional 1 Consider the energy w D in [19, equation (4)] for the crystallographic class C n (n = 3, 4, 6). The bulk DMI energy (5) corresponds to the last two terms of w D , i.e., D 1 = 0 and D 2 = D 3 = −D (crystallographic subclass D n ). The interfacial DMI energy (6) corresponds to the first term of w D , i.e., D 1 = D and D 2 = D 3 = 0 (crystallographic subclass C nv ). 2 We rescale the time according to the transformation t = γ 0 M s t. We define the rescaled effective field by h eff = H eff /M s and the rescaled energy by E = E /(μ 0 M 2 s ). However, to simplify the notation, we neglect all -superscripts.

according to the relation
The positive quantities ex = 2A/(μ 0 M 2 s ) and dm = 2D/(μ 0 M 2 s ) denote the exchange length and the DMI length (both measured in m), respectively.
Since ∇ × m L 2 ( ) ≤ √ 2 ∇m L 2 ( ) , using the weighted Young inequality for any a, b ∈ R and ε > 0, (10) it is easy to see that the energy (8) satisfies the condition For any T > 0, we define the space-time cylinder by T := × (0, T ). Moreover, we denote by ·, · the scalar product in L 2 ( ), by ·, · the duality pairing between H −1/2 ( ) and H 1/2 ( ), and by γ T : H (curl, ) → H −1/2 ( ) the tangential trace operator, which satisfies γ T [u] = u × n| for any smooth function u as well as the Green formula We conclude this section by extending the notion of a weak solution introduced in [9] to the present setting.

∂ t m(t), ϕ(t) dt
(iv) It holds that The variational formulation (13) comes from a weak formulation of (7a) in the space-time domain. The boundary conditions (7b) are enforced as natural boundary conditions. In particular, the term with ·, · arises from integrating by parts the exchange contribution and using (7b). The energy inequality (14) is a weak counterpart of the dissipative energy law satisfied by any sufficiently smooth solution of (7a)-(7c).

Remark 1
Taking the scalar product of (7a) with m, we deduce that m · ∂ t m = 0. In particular, since ∂ t |m| 2 = 2 m · ∂ t m = 0, it follows that a sufficiently smooth solution of (7a) satisfies the constraint |m| = 1, provided that it is satisfied by the initial condition.

Numerical algorithms and main result
In this section, we introduce three algorithms for the numerical approximation of the problem discussed in Section 2.2 and we state the main convergence result.

Preliminaries
For the time discretization, given an integer N > 0 and a final time T > 0, we consider a uniform partition of the time interval (0, T ) with time-step size k := T /N, i.e., t i := ik for all 0 ≤ i ≤ N. For the spatial discretization, we assume to be a polyhedral domain with Lipschitz boundary and consider a κ-quasi-uniform family {T h } h>0 of regular tetrahedral meshes of parametrized by the mesh size h > 0, i.e., there exists κ ≥ 1, independent of h, such that T h is κ-shape-regular and κ −1 h ≤ diam(K) 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 linear polynomials on K. We consider the space of piecewise linear and globally continuous functions from to R. The classi- We assume that all off-diagonal entries of the so-called stiffness matrix are nonpositive, i.e., This requirement is usually referred to as angle condition, since it is satisfied if the measure of all dihedral angles of all tetrahedra of the mesh is less than or equal to π/2. Any solution of LLG is characterized by the nonconvex pointwise constraint |m| = 1 and by the orthogonality property m · ∂ t m = 0. To mimic these properties at the discrete level, we require them to be satisfied only at the nodes of the mesh. To this end, we introduce the set of admissible discrete magnetizations which we call the discrete tangent space of ψ h .

Three tangent plane integrators
Using the well-known formula (7a) can be formally rewritten in the form Observing that this equation is linear with respect to the time derivative ∂ t m, we introduce the free variable v = ∂ t m. For any t ∈ (0, T ), v(t) belongs to the tangent space of S 2 at m(t). Taking this orthogonality and the expression (9) of the effective field into account, we obtain the following variational formulation: for all φ ∈ H 1 ( ) satisfying m(t) · φ = 0 a.e. in . To obtain (18), the boundary integral which arises from integrating by parts the exchange contribution and using the boundary conditions (7b) is rewritten as a volume integral by using (12). Note that, since the test function φ belongs to the tangent space of the sphere at m(t), in (18) the term corresponding to the last term (strongly nonlinear in m) on the righthand side of (17) vanishes. For any time-step 0 ≤ i ≤ N − 1, given the approximate current magnetization ) via a first-order time-stepping. To ensure that the discrete magnetization belongs to the set of admissible discrete magnetizations M h , the nodal projection is applied.
The resulting scheme, summarized in the following algorithm, extends the method proposed by [5] to the present situation. Algorithm 1 1st-order tangent plane scheme, TPS1.
In (19), the parameter 0 ≤ θ ≤ 1 modulates the "degree of implicitness" of the method in the treatment of the leading-order exchange contribution of the effective field.
In the following algorithm, we state a projection-free variant of Algorithm 1, where step (ii) is replaced by a simple linear first-order time-stepping. Note that, omitting the nodal projection, the pointwise constraint |m| = 1 is not explicitly enforced by the numerical scheme.
The idea of removing the nodal projection from the tangent plane scheme goes back to [3] for LLG and has been inspired by [13], where the same principle is applied to a certain class of geometrically constrained PDEs, e.g., the harmonic map heat flow.
In [7], the authors extend the tangent plane scheme of [5] to improve the formal convergence order in time of the method. If the tangential update v takes the form where P u : R 3 → span(u) ⊥ denotes the orthogonal projection onto span(u) ⊥ for any u ∈ S 2 , a sufficiently smooth solution m of LLG satisfies the Taylor expansion Differentiating (17) with respect to time and proceeding as in [7, Section 6], one obtains that, up to a residual term of order O(k 2 ), the tangential update (21) can be characterized as the solution of the following variational problem: To obtain an effective numerical method, we use the same predictor-corrector approach used for Algorithms 1-2: For However, in order to obtain a well-defined scheme, following [7], we perform two higher-order modifications of (23). Firstly, to ensure the well-posedness of the variational problem, we proceed as follows: Given M > 0, we define the cutoff function W M : R → R by By construction, it holds that (see, e.g., [24,Lemma 12]). In the variational formulation (23), we then replace Note that the function λ(m) is also the Lagrange multiplier associated with the constraint |m| = 1 in the constrained minimization problem (4). If M > 0 is sufficiently large, this modification introduces a consistency error of order O(k 2 ) in (23). In particular, to ensure this, we define M : R >0 → R >0 by M(k) := |k log k| −1 for all k > 0. Note that M satisfies the convergence properties Secondly, in the variational formulation (23), we replace 2 ex 2 k ∇v(t), ∇φ by where the stabilization function ρ : R >0 → R >0 is defined by ρ(k) := |k log k| for all k > 0. This artificial stabilization introduces a formal consistency error of order O(k 2−ε ) for any 0 < ε < 1.
In the following algorithm, we summarize the proposed extension of the tangent plane scheme of [7] to the present setting.
For the proof that the three proposed algorithms are well-posed (if the time-step size is sufficiently small in the case of TPS2), we refer to Proposition 1 below.

Remark 2
The natural starting point for a hypothetical projection-free version of TPS2 would be the expansion for which a nontangential update of the form would be required. In particular, it is not clear how to apply the tangent plane paradigm to this situation, where the update v(t) has a nonzero component parallel to m(t) in general, i.e., m(t) · v(t) = k 2 ∂ tt m(t) · m(t) = 0, which needs to be taken into account in order to achieve a second-order accuracy.

Convergence result
Using the sequences of approximations {m i h } 0≤i≤N and {v i h } 0≤i≤N −1 obtained from any algorithm, we define the piecewise linear time reconstruction m hk and the piecewise constant time reconstructions m ± hk and v − hk by The following theorem, which is the main result of the paper, states that the time reconstructions obtained by the three algorithms converge in an appropriate sense towards a weak solution of (7a)-(7c) as h, k → 0.

Theorem 1 Let the approximate initial condition satisfy the convergence property
Moreover, for each algorithm, consider the following specific assumptions: Then, for each algorithm, there exist a weak solution m of (7a)-(7c) and a subsequence of {m hk } which converges weakly in H 1 ( T ) towards m as h, k → 0.
Remark 3 (i) For the sake of brevity, we have considered the case in which the energy consists of only exchange and DMI. Adopting the implicit-explicit approaches of [3,8,22,24], the schemes and the convergence result of Theorem 1 can be extended to the case in which also the standard lower-order energy terms (and their discretizations) are included into the setting. (ii) One important aspect of the research on numerical integrators for LLG is related to the development of unconditionally convergent methods, for which the numerical analysis does not require to impose any CFL-type condition on h and k. Theorem 1 states that this goal is achieved by PF-TPS1. For TPS1 and TPS2, our analysis requires a mild CFL condition, which arises from the use of the nodal projection and the presence of the DMI. If the energy comprises only the standard micromagnetic contributions (exchange, uniaxial anisotropy, Zeeman, and magnetostatic), but no DMI, then the convergence towards a weak solution of LLG is unconditional also for TPS1 (for 1/2 < θ ≤ 1) and TPS2 (see [7,8,22,24]). (iii) Since the treatment of the DMI requires the CFL condition k/ h → 0 as h, k → 0 in our analysis, the result of Theorem 1 holds for TPS2 also without artificial stabilization, i.e., ρ ≡ 0 (see [7,24] for more details). In this case, TPS2 is of full second order in time.
(iv) In Theorem 1, we state the best result that we are able to prove in terms of stability, i.e., with the weakest CFL condition on the discretization parameters. The same convergence result can be also established at the price of more severe restrictions. In particular, the result of Theorem 1 holds (a) without angle condition (15) for TPS1 and TPS2, if

Numerical experiments
Before proceeding with the convergence analysis, we aim to show the effectivity of the proposed algorithms with three numerical experiments. The computations presented in this section have been performed with our micromagnetic software Commics [1,51]. Our Python code is based on the open-source finite element library Netgen/NGSolve [2]. The computation of the stray field, i.e., the numerical solution of the transmission problem (3a)-(3e), is based on the hybrid FEM-BEM method of [31], which requires the evaluation of the double-layer integral operator associated with the Laplace equation (see, e.g., [22,Section 4.4.1] or [52,Algorithm 12]). This part of the code exploits the open-source Galerkin boundary element library BEM++ [59]. For all three schemes, to discretize the classical lower-order contributions (anisotropy, Zeeman, and magnetostatic), we follow the explicit approaches of [3,22,24]. Magnetization configurations are visualized with ParaView [4].
For the spatial discretization, we consider three types of tetrahedral meshes (see Fig. 2). For meshes of type I, the domain is first uniformly decomposed into cubes. Then, each cube is split into six tetrahedra in such a way that any tetrahedron has three mutually perpendicular edges. Any mesh of this type satisfies the angle condition (15) (see [12,Lemma 3.5]). Meshes of type II are obtained from the previous Meshes of this type do not satisfy (15). For type III, we consider unstructured meshes obtained with Netgen, which are generated with the advancing front method (see [58] for details) and in general do not satisfy (15). In Fig. 3, we plot the time evolutions of the third component of the spatially averaged magnetization of the sample, i.e., m 3 (t) = | | −1 m 3 (x, t) dx, and the energy (8) obtained with the three algorithms for different mesh types and sizes, and a constant time-step size k = γ 0 e−7 ≈ 0.0221. In Fig. 3a and b, we compare the results obtained with TPS1 (θ = 1), PF-TPS1 (θ = 1), and TPS2 for the mesh types I-III. For each mesh type, we consider a mesh size of h ≈ 3.46 nm. Note that the meshes of types II-III violate (15). In Fig. 3c, we compare the results obtained with TPS1 (θ = 1/2) and TPS2. We consider a structured mesh of type I and compare the results obtained for different mesh sizes.
Although the convergence result of Theorem 1 does not cover meshes of types II-III for TPS1 and TPS2, the numerical results show that, in terms of stability, the methods behave identically, independently of the mesh type used. To better understand this aspect, we also monitored a posteriori the validity of the inequality which is the inequality effectively used in the stability analysis of TPS1 and TPS2 (see Propositions 2 and 4 below), respectively. It turned out that the inequality is always satisfied, even for meshes violating (15). The omission of the nodal projection in PF-TPS1 manifests itself as a phase error in the evolution of m 3 accumulating over time (see Fig. 3a), and as a lower energy level of the final magnetization configuration (see Fig. 3b). However, the overall qualitative outcome of the experiment is preserved.
The results also show that TPS1 with θ = 1 is more dissipative than TPS2. However, the choice of θ = 1/2, which would be favorable from an energetic point of view (no artificial damping), is not feasible, because it affects the stability of the   scheme (see Fig. 3b). The instability is more severe for smaller mesh sizes, giving numerical evidence of the CFL condition required for stability in this case (see Remark 3(ii)). Finally, in Fig. 4, we study the violation of the unit-length constraint which occurs for PF-TPS1. We consider a structured mesh of type I with mesh size h ≈ 4.33 nm and plot the error I h m + hk (T ) 2 − 1 L 1 ( ) for different time-step sizes k. Note that this error is identically zero for TPS1 and TPS2, because of the nodal projection. We observe a linear dependence of the error on k which is in total agreement with the theory, see estimate (50) in the proof of Proposition 6 below.

L 1 -error
For numerical experiments testing the experimental convergence rates in time of the schemes (in the absence of DMI), we refer to [56, Section 6.2.1] (for TPS1 and PF-TPS1) and [24, Section 7.1] (for TPS1 and TPS2). There, the observed rates with respect to a reference solution match the formal consistency error of the schemes, i.e., first-order convergence for TPS1 and PF-TPS1, second-order convergence for TPS2. A similar numerical study for the present model problem (which includes DMI) confirms the first-order convergence for TPS1 and PF-TPS1, but does not reveal a full second-order convergence for TPS2. We believe that this is due to a lack of regularity of the solution in time.

Stability of isolated skyrmions in nanodisks
We reproduce a numerical experiment from [57]. We investigate the relaxed states of a thin nanodisk of diameter 80 nm (aligned with x 1 x 2 -plane) and thickness 0.4 nm (x 3 -direction) centered at (0, 0, 0) for different values of the DMI constant and initial conditions. The effective field in (1) consists of exchange interaction, perpendicular uniaxial anisotropy, interfacial DMI, and stray field, i.e., The values of the involved material parameters mimic those of cobalt: M s = 5.8 ·  For all simulations, we choose T = 2 ns for D = 0, . . . , 6 mJ/m 2 and T = 5 ns for D = 7, 8 mJ/m 2 , which experimentally turn out to be sufficiently large times to relax the system. The computational domain is discretized by a regular partition consisting of 32,575 tetrahedra (mesh size of 1 nm). For the time discretization, we consider a uniform partition of the time interval (0, T ) with a time-step size of 0.1 ps.
In the case of the uniform out-of-plane initial condition, the stable state remains a quasi-uniform ferromagnetic state for the values D = 0, . . . , 6 mJ/m 2 and turns into a multidomain state for the values D = 7, 8 mJ/m 2 (see Fig. 5). For D = 0, . . . , 6 mJ/m 2 , the slight decrease of the total energy for increasing values of D corresponds to an inward tilt of the magnetization on the boundary of the disk. In the case of the skyrmion-like initial condition, the stable state is a quasi-uniform ferromagnetic state for the values D = 0, 1, 2 mJ/m 2 , a skyrmion for the values D = 3, . . . , 6 mJ/m 2 , and a multidomain state for the values D = 7, 8 mJ/m 2 (see Fig. 6). The skyrmion size, i.e., the diameter of the circle {m 3 = 0} in the x 1 x 2 -plane, increases from the minimum value of circa 14 nm for D = 3 mJ/m 2 to the maximum value of circa 48 nm for D = 6 mJ/m 2 . As observed in [57], the fact that for D = 3, . . . , 6 mJ/m 2 , which are realistic values for the DMI constant, both the ferromagnetic state and the skyrmion state can be stabilized is very relevant for applications. Indeed, this bistability can be exploited to code the information in future recording devices (the presence and the absence of a skyrmion can be used to encode one bit) (see, e.g., [29,60]).
In Fig. 7, we plot the total energy of the relaxed state for different values of the DMI constant. The energy values obtained with TPS1 (the results refer to the case θ = 1) and TPS2 are in perfect quantitative agreement with each other and with those reported in [57, Fig. 1]. The use of PF-TPS1 preserves the qualitative outcome of the experiment, but the quantitative agreement of the energy values with those of [57, Fig. 1], as a result of the violation of the pointwise constraint |m| = 1, is inevitably lost.

Field-induced dynamics of skyrmions in nanodisks
We numerically investigate the stability and the induced dynamics of isolated magnetic skyrmions in helimagnetic materials in response to an applied field pulse. The sample under consideration is a magnetic nanodisk of diameter 140 nm (x 1 x 2 -plane) and thickness 10 nm (x 3 -direction). The effective field in (1) consists of exchange interaction, bulk DMI, applied external field, and stray field, i.e., We use the material parameters of iron-germanium (FeGe), i.e., A = 8.78 · 10 −12 J/m, D = 1.58 · 10 −3 J/m 2 , and M s = 3.84 · 10 5 A/m (see, e.g., [16]). The initial condition for our experiment is obtained by setting H ext ≡ (0, 0, 0) and relaxing a uniform out-of-plane ferromagnetic state m 0 ≡ (0, 0, 1) for 3 ns. For the relaxation process, we choose the large value α = 1 for the Gilbert damping constant, since we are not interested in the precise magnetization dynamics. The resulting relaxed state is the skyrmion depicted in Fig. 8a. Starting from this configuration, we perturb the system from its equilibrium by applying an in-plane field pulse H ext (t) = (H (t), 0, 0) of maximum intensity H max > 0 for 150 ps (see Fig. 9). Then, we turn off the applied external field, i.e., H ext ≡ (0, 0, 0), and let the system relax to equilibrium. In order to capture all possible excitation modes, during the application of the field and the subsequent relaxation process, we set the value of the Gilbert damping constant to α = 0.002, which is considerably smaller than the experimental value of α = 0.28 measured for FeGe (see [16]). To probe the limit of the stability of the skyrmion, we test different values for the maximum intensity of the field H max , namely μ 0 H max = 1, 2, 5, 10, 20, 50, 100, 200 mT. For the spatial discretization, we consider a regular partition of the nanodisk consisting of 36,501 In Fig. 10, we plot the first 10.15 ns of the time evolution of the second component of the spatially averaged magnetization of the sample, i.e., m 2 (t) = | | −1 m 2 (x, t) dx. We see that, for the values μ 0 H max = 1, 2, 5, 10, 20, 50 mT, the induced dynamics is a periodic damped precession of the skyrmion around the center of the sample, which comes back to the initial stable configuration by the relaxation process. As expected, both the deflection from the stable symmetric initial state and the amplitude of the oscillations increase for larger values of H max . For the value μ 0 H max = 100 mT, the skyrmion is critically deformed by the applied field pulse, but the initial stable configuration is recovered by the relaxation process. Note that a different oscillating mode comes into play in this case. For μ 0 H max = 200 mT, the skyrmion is destroyed. After approximately 3.5 ns of chaotic dynamics, the magnetization configuration turns into a horseshoe state which then starts to rotate around the center of the sample (see Fig. 8b).
As observed for the experiment of Section 4.2, also in this case, the results obtained with TPS1 and TPS2 are in full quantitative agreement with each other. The use of PF-TPS1 preserves the qualitative outcome of the experiment, but the computed quantities, e.g., the amplitudes of the oscillations depicted in Fig. 10, are slightly perturbed.
The presented experiment is a preliminary study to investigate the stability and the dynamics of a skyrmion in the presence of a field pulse. This is done to explore the possibility to use time-resolved scanning Kerr microscopy [37] which is based on the interplay between laser and field pulses to directly map the dynamics of magnetic skyrmions [38].

Convergence analysis
In this section, we show that all proposed algorithms are well-posed and we present the proof of Theorem 1. To establish the convergence result, we use the standard energy method for proving existence of solutions of linear second-order parabolic problems (see, e.g., [27,Section 7.1.2]). The main difference is that, following [3,5,7], the construction of approximate solutions is not obtained by applying the Galerkin method based on a basis of appropriately normalized eigenfunctions of the Laplace operator, but rather by using the finite element solutions delivered by the numerical schemes.

Preliminaries
We introduce some further notation and collect some auxiliary results. We consider the nodal interpolant I h : It is well known that, for κ-shape-regular meshes and any integer 0 ≤ m ≤ 2, the nodal interpolant satisfies the approximation property where the constant C > 0 depends only on κ. We denote the vector-valued realization of the nodal interpolant by I h : C 0 ( ) → S 1 (T h ) 3 . The following classical inverse estimate requires the quasi-uniformity of the underlying family of meshes: For any 1 ≤ p ≤ ∞, it holds that where C inv > 0 depends only on κ and p. Standard scaling arguments show that, for any 1 ≤ p < ∞, we have the discrete norm equivalence where C norm > 0 depends only on κ and p.
To ensure that the approximate magnetization belongs to M h , TPS1 and TPS2 employ the nodal projection. We refer to [12,Lemma 3.2] for the proof that the nodal projection φ h → I h φ h / φ h does not increase the exchange energy of a discrete function if the underlying mesh fulfills a weak acuteness condition. Specifically, the angle condition (15) Owing to the application of the nodal projection, for all 0 ≤ i ≤ N − 1, the iterates of TPS1 and TPS2 satisfy m i+1 h L ∞ ( ) = 1 and the geometric estimates for any z h ∈ N h (see [6,14]). With (31), these nodewise inequalities are turned into where C geo > 0 depends only on κ and p. In the case of PF-TPS1, where the nodal projection is omitted, the geometric estimates (33a) and (33b) become trivial, but the equality m i+1 h L ∞ ( ) = 1 does not hold anymore. However, for any 1 ≤ j ≤ N, the linear time-stepping yields the recursive relation Together with (31), this leads to the estimate

Well-posedness
In the following proposition, we prove that the three proposed algorithms are all well-posed. (19) and (20). There exists a threshold time-step size k 0 > 0, which depends only on α, ex , and dm , such that, if k ≤ k 0 , there exists a unique solution v i h ∈ K h (m i h ) of (26). The time-steppings of all algorithms are well defined.
Proof For any 0 ≤ i ≤ N − 1, let a i hk (·, ·) be the bilinear form appearing on the left-hand side of (19) and (20). For any φ h ∈ S 1 (T h ) 3 , it holds that Hence, the bilinear form is elliptic, even on the full space S 1 (T h ) 3 (19) and (20) thus follow from the Lax-Milgram theorem.

. Existence and uniqueness of the solution
Similarly, for any 0 ≤ i ≤ N − 1, let b i hk (·, ·) denote the bilinear form on the lefthand side of (26). For any ε > 0 and φ h ∈ S 1 (T h ) 3 , using (10) and (24), we deduce that We choose ε = 2 ex / dm . With the properties (25) of the cutoff function M(k), it follows that both coefficients in front of the norms are positive if the time-step size k is sufficiently small. Then, the bilinear form is elliptic and (26) admits a unique solution v i h ∈ K h (m i h ). The linear time-stepping in step (ii) of PF-TPS1 is clearly well defined. In the case of TPS1 and TPS2, which include the nodal projection, since The time-steppings of TPS1 and TPS2 are therefore also well defined.

Discrete energy law and stability
In this section, we establish the discrete energy laws and study the stability for the discrete iterates delivered by the algorithms. We first observe that, given C 0 > 0, assumption (28) provides some h 0 > 0 such that In the following proposition, we prove the result for TPS1.
Proposition 2 (Discrete energy law and stability of TPS1) Let 1 ≤ j ≤ N. There exists a constant C > 0, which depends only on κ and dm , such that the iterates of TPS1 satisfy the discrete energy law (36) Moreover, there exists a constant C > 0, which depends only on α, κ, and dm , such that, if h ≤ h 0 and k ≤ C h, the iterates of TPS1 satisfy the stability estimate where the constant C > 0 depends only on α, C 0 , κ, ex , dm , and | |.
Since the angle condition (15) is satisfied, we obtain that Hence, it follows that We obtain the energy inequality With some simple algebraic manipulations, we rewrite the last four terms of the right-hand side (those which involve the curl operator) as Using the inverse estimate (30) and the geometric estimates (33a) and (33b), we infer that Similarly, it holds that Summation over 0 ≤ i ≤ j − 1 leads to (36).
To show (37), we first note that m j h L ∞ ( ) = 1 yields that We multiply (40) by 2 dm / 2 ex and add the resulting equation to (36). Using the characterization (11) of the energy, we obtain that Let C = α/(2C) and k ≤ C h. Since θ ≥ 1/2, all terms on the left-hand side are nonnegative. We obtain (37), where the constant C > 0 (which we do not compute explicitly) depends only on α, C 0 , κ, ex , dm , and | |.
In the following proposition, we prove the corresponding result for PF-TPS1.
Proposition 3 (Discrete energy law and stability of PF-TPS1) Let 1 ≤ j ≤ N and 1/2 < θ ≤ 1. The iterates of PF-TPS1 satisfy the discrete energy law Moreover, there exists a threshold time-step size k 0 > 0, which depends only on α, κ, ex , dm , and θ , such that, if h ≤ h 0 and k ≤ k 0 , the iterates of PF-TPS1 satisfy the stability estimate where the constant C > 0 depends only on α, C 0 , κ, ex , dm , and θ .
We follow the argument of the proof of Proposition 2: Due to the linear time-stepping of PF-TPS1, all the computations in (38a)-(38c) hold with equality sign and without resorting to the angle condition (15). Moreover, all but the last term on the right-hand side of (39) vanish. As a result, we obtain the energy identity Summing this identity over 0 ≤ i ≤ j − 1, we prove (41). To estimate the right-hand side, we apply the weighted Young inequality (10), which, for any ε > 0, yields that Choosing ε = 2 ex (θ − 1/2)/ dm , the characterization (11) of the energy shows that We multiply (35) ]}, then all terms on the left-hand side are nonnegative. This leads to (42), where the (explicitly computable) constant C > 0 depends only on α, C 0 , κ, ex , dm , and θ .
Finally, in the following proposition, we prove the stability result for TPS2. Proposition 4 (Discrete energy law and stability of TPS2) Let 1 ≤ j ≤ N. Suppose that the time-step size k is sufficiently small, so that (26) is well-posed by Proposition 1. Then, there exists a constant C > 0, which depends only on κ and dm , such that the iterates of TPS2 satisfy the discrete energy law Moreover, there exists a threshold time-step size k 0 > 0, which depends only on α, ex , and dm , and a constant C > 0, which depends only on α, κ, and dm , such that, if h ≤ h 0 , k ≤ k 0 , and k ≤ C h, the iterates of TPS2 satisfy the stability estimate where the constant C > 0 depends only on α, C 0 , κ, ex , dm , and | |.
Proof Let 0 ≤ i ≤ j − 1. We follow step by step the argument of the proof of Proposition 2 to obtain the inequality We reformulate the terms of the right-hand side which involve the curl operator, i.e., and proceed with their direct estimation: Using (30), (33a) and (33b), we obtain that

Extraction of weakly convergent subsequences
Exploiting the established stability estimates of the three algorithms, we are now able to prove that the time reconstructions defined by (27) are uniformly bounded.
Proof The estimate (46) follows directly from (37), (42), and (45). For TPS1 and TPS2, also the geometric estimate (33a) is used to conclude that ∂ t m hk L 2 ( T ) ≤ C. The convergence (47) for PF-TPS1 follows from Proposition 3. Indeed, it holds that For TPS1 (resp. TPS2), we first resort to an inverse estimate to obtain that The result then follows from Proposition 2 (resp. Proposition 4) and the fact that k/ h → 0 as h, k → 0 by assumption.
With this result, we can now extract weakly convergent subsequences.
Proof For the sake of clarity, we divide the proof into three steps.
Let v ∈ L 2 ( T ) such that v − hk v in L 2 ( T ). In the case of TPS1 and TPS2, which includes the nodal projection, it holds that k, which shows that v = ∂ t m a.e. in T . For PF-TPS1, the result directly follows from the equality m i+1 Step 3: m satisfies |m| = 1 a.e. in T .
In the case of TPS1 and TPS2, since I h m i h 2 = 1 and ∇m i h is piecewise constant for all 0 ≤ i ≤ N − 1, it holds that (46) h, which yields the convergence m − hk 2 → 1 in L 2 ( T ). Since m − hk → m pointwise a.e. in T , we deduce that |m| = 1 a.e. in T .
In the case of PF-TPS1, we start with a triangle inequality, which shows that The first two terms on the right-hand side converge to 0. Indeed, on the one hand, it holds that and m + hk → m in L 2 ( T ). On the other hand, using the approximation properties (29) of the nodal interpolant and the fact that ∇m i+1 h is piecewise constant, one shows that To conclude, it remains to show that For any t ∈ (0, T ), let 0 ≤ i ≤ N − 1 such that t ∈ [t i , t i+1 ). It holds that For the first term on the right-hand side, it holds that ≤ Ck.
Using the approximation properties of I h , we estimate the second term by Finally, since m 0 = 1 a.e. in by assumption, the third term satisfies that Thanks to (28), this yields the convergence m 0 h 2 → 1 in L 1 ( ). Altogether, this proves (49) and thus concludes the proof.

Identification of the limit with a weak solution of LLG
We start with establishing an auxiliary convergence result for the time reconstructions obtained by PF-TPS1. Since m hk → m in L 2 (0, T ; H s ( )) by Proposition 6, the result follows from (46)- (47).
We have collected all ingredients to finalize the proof of Theorem 1.
Proof of Theorem 1 By Proposition 6, for any algorithm, we deduce the desired convergence towards a function m ∈ L ∞ (0, T ; H 1 ( )) ∩ H 1 (0, T ; L 2 ( )) satisfying |m| = 1 a.e. in T . Since m hk m in H 1 ( T ), we also have the weak convergence of the traces, i.e., m hk (0) m(0) in H 1/2 ( ). By assumption (28), we deduce that m(0) = m 0 in the sense of traces. It remains to show that m fulfills the variational formulation (13) and the energy inequality (14). For the sake of clarity, we consider the three algorithms separately.
Step 1: Proof of the result for TPS1.

) × ϕ(t)] ∈ K h (m i h ). Integrating in
The proof follows the lines of the one for TPS1 discussed in Step 1. In the proof of the variational formulation (13), the only difference is the convergence of the second term on the left-hand side of (51), which is more subtle here, since omitting the nodal projection the uniform boundedness of m − hk in L ∞ ( T ) is lost. To show the desired convergence, we start with recalling the so-called Lagrange identity · d)(b · c) for all a, b, c Thanks to Lemma 1, we deduce that m − hk 2 → 1 in L 2 ( T ) as h, k → 0. Together with the weak convergence v − hk · ϕ ∂ t m · ϕ in L 2 ( T ), it follows that Finally, passing the discrete energy law (41) to the limit as h, k → 0, thanks to (28), (47), the available convergence results (48a)-(48g), and standard lower semicontinuity arguments, we obtain (14).
Step 3: Proof of the result for TPS2.
The verification of the variational formulation (13) follows by the same method used in Step 1 for TPS1. Given an arbitrary ϕ ∈ C ∞ ( T ), for any 0 ≤ i ≤ N − 1 and t ∈ (t i , t i+1 ), we test (26) where, in analogy with (27), we define the piecewise time reconstruction λ − hk by λ − hk (t) := λ i h for all 0 ≤ i ≤ N − 1 and t ∈ [t i , t i+1 ). With the available convergence result (48a)-(48g) and the convergence properties of W M(k) (·) and ρ(·), each of the three terms on the left-hand side converges towards the corresponding term of (13) as h, k → 0 (see [7] for details). We discuss the convergence of the two terms on the right-hand side. Since ∇ × m − Similarly, we show that the second term on the right-hand side converges towards On the other hand, we have that This proves (13) for any smooth test function ϕ. By density, we obtain the desired result. Finally, the energy inequality (14) is obtained by passing to the limit as h, k → 0 the discrete energy law (44) and using standard lower semicontinuity arguments.