Spin-polarized transport in ferromagnetic multilayers: An unconditionally convergent FEM integrator

We propose and analyze a decoupled time-marching scheme for the coupling of the Landau-Lifshitz-Gilbert equation with a quasilinear diffusion equation for the spin accumulation. This model describes the interplay of magnetization and electron spin accumulation in magnetic and non-magnetic multilayer structures. Despite the strong nonlinearity of the overall PDE system, the proposed integrator requires only the solution of two linear systems per time-step. Unconditional convergence of the integrator towards weak solutions is proved.


Introduction
The interaction between electric current and magnetization in magnetic nanostructure devices and the control of this interaction have been realized through the prediction of the spin-transfer torque by Slonczewski and Berger [14,31]. The transfer of spin angular momentum between the spin-polarized electrical current and the local magnetization has been observed in various magnetic devices, such as metallic spin-valves systems, magnetic tunnel junctions, and magnetic domain walls in permalloy nanowires [29,33]. Based on these experiments, a number of technological applications have been proposed, e.g., STT-MRAMs, racetrack memories, and magnetic vortex oscillators [26,27].
The fundamental physics underlying these phenomena is understood as due to a spin torque that arises from the transfer of the spin angular momentum between conduction free electrons and magnetization. In the original works of Berger and Slonczewski [14,31], a homogeneous spin accumulation is assumed due to a current which flows through a first magnetic layer perpendicular to the interface into a second magnetic layer. The spin torque effect leads to an interaction between the spin-polarized current and the magnetization in the second layer. For magnetic multilayers it has been shown that a proper description of the magnetoresistance is essential to take into account the interplay between successive interfaces [22,34,32]. In order to calculate the spin torque transfer, the spin transport properties have to be calculated far beyond the interface.
The original model of Berger and Slonczewski has been extended by taking into account the diffusion process of the spin accumulation by Shpiro et al. for one-dimensional systems [30] and by García-Cervera and Wang [21,20] for three-dimensional systems. There, the overall system of PDEs (SDLLG) is a quasilinear diffusion equation for the evolution of the spin accumulation coupled to the Landau-Lifshitz-Gilbert equation (LLG) for the magnetization dynamics. Existence of global weak solutions to LLG goes back to [5], while, in the same spirit, existence of global weak solutions to SDLLG is proved in [21].
The reliable numerical integration of LLG (and, in particular, SDLLG) faces several challenges due to the nonuniqueness of weak solutions, the explicit nonlinearity, and an inherent nonconvex modulus constraint. Numerical approximation schemes for weak solutions of LLG are first proposed in [3,12]. First unconditional convergence results can be found in [13,2], which consider the small-particle limit of LLG with exchange only. On the one hand, the integrator of [13] relies on the midpoint rule and reduced integration, and thus has to solve one nonlinear system of equations per time-step. On the other hand, the tangent plane integrator of [2], which extends the prior works [3,12], relies on a reformulation of LLG which is solved for the discrete time derivative. Each time-step consists of the solution of one linear system of equations plus nodal projection. It has been generalized to linear-implicit time integration and full effective field in [4,15].
Numerical integration of the coupling of LLG to other time-dependent PDEs has been analyzed in [6,7] for the full Maxwell equations (MLLG), in [25,24] for the eddy current formulation, and in [8] for LLG with magnetostriction. While [6] analyzes an extension of the midpoint scheme of [13], the works [7,25,24,8] extend the tangent plane scheme from [2], and emphasis is on the decoupling of the time-marching scheme in [7,24,8].
In the models and works mentioned, e.g., MLLG, the coupling of LLG and Maxwell equations is weak in the sense that the magnetization of LLG only contributes to the right-hand side of the Maxwell system, while the magnetic field from the Maxwell equations gives a contribution to the effective field of LLG. In SDLLG the principal part of the differential operator of the spin diffusion equation depends nonlinearly on the magnetization. A first numerical integrator for SDLLG is proposed and empirically validated in [20]. While this scheme appears to be unconditionally stable, the work does not prove convergence of the discrete solution towards a weak solution of SDLLG.
In our work, we extend the tangent plane integrator to SDLLG and prove unconditional convergence. Altogether, the contributions of the current work can be summarized as follows: • The proposed integrator is proven to converge (at least for a subsequence) towards a weak solution of SDLLG. This convergence is unconditional, i.e., there is no CFL-type coupling of the time and space discretizations. Despite the nonlinearity of SDLLG, each timestep requires only the solution of two successive linear systems, one for (the discrete time derivative of) the magnetization and one for the spin accumulation.
• Unlike prior work on the tangent plane integrator, we adopt an idea from [11] and show that the nodal projection step of the tangent plane scheme is not necessary. In particular and unlike the cited works, our analysis can therefore avoid a technical angle condition on the triangulations used. This result also transfers to the models and analysis of [2,4,15,7,8,24] and simplifies their (extended) tangent plane integrators.

Outline
The paper is organized as follows: In Section 2, we introduce and accurately describe the mathematical model, see (9) for the nondimensional formulation of SDLLG. In Section 3, we formulate a decoupled time-marching scheme (Algorithm 6) for the numerical integration of SDLLG and prove its well-posedness (Proposition 9). Section 4 contains the main result of our work (Theorem 12), which states unconditional convergence of the scheme towards weak solutions of SDLLG. Following [21], weak solutions of SDLLG have finite energy. In Section 5, we prove that any weak limit obtained by the proposed numerical integrator shows the same energy behavior as formal strong solutions of SDLLG (Theorem 24). Numerical examples as well as the empirical validation of the proposed algorithm are postponed to a forthcoming paper [1].

Notation
We use the standard notation [16] for Lebesgue and Sobolev spaces and norms. For any domain D, we denote the L 2 scalar product by (r, s) D = D rs for all r, s ∈ L 2 (D). In the case of (spaces of) vector-valued functions, we use bold letters. For a sequence {x n } n≥1 in a Banach space X and x ∈ X, we write x n → x (resp. x n x) in X if the sequence converges strongly (resp. weakly) to x in X. Similarly, we write x n sub −→ x (resp. x n sub − x) in X if there exists a subsequence of {x n } n≥1 which converges strongly (resp. weakly) to x in X. Throughout the paper, C denotes a generic positive constant, independent of the discretization parameters, not necessarily the same at each occurrence. Alternatively, we write A B to abbreviate A ≤ C B. Given a, b ∈ R 3 , we denote by a ⊗ b ∈ R 3×3 the tensor product defined by (a ⊗ b) jk = a j b k for all 1 ≤ j, k ≤ 3. By |·|, we denote both the Frobenius norm of a matrix and the Euclidean norm of a vector. Since the meaning is clear from the argument, this does not lead to any ambiguity.

Model problem
In this section, we present the mathematical model, for which we introduce a nondimensional formulation, as well as the notion of a weak solution. We use physical units in the International System of Units (SI).

Physical background
We consider a magnetic multilayer. Let ω ⊂ Ω be polyhedral Lipschitz domains in R 3 , where Ω corresponds to the volume occupied by the multilayer, and ω corresponds to the ferromagnetic part. A possible experimental setup is shown in Figure 2 In micromagnetics, the quantity of interest is the magnetization M : ω T → R 3 , measured in ampere per meter (A/m). If the temperature is constant and far below from the Curie temperature of the ferromagnetic material, M is a vector field of constant modulus |M| = M s , with M s being the saturation magnetization (in A/m). In the absence of spin currents, the dynamics of M is Figure 1: Schematic of a magnetic nanopillar structure (trilayer) consisting of two ferromagnetic films, ω 1 and ω 2 , separated by a nonmagnetic interlayer ω 0 . The current is assumed to flow perpendicularly from ω 1 to a bottom electrode connected to ω 2 . In this case, ω = ω 1 ∪ ω 2 and Ω = ω 1 ∪ ω 0 ∪ ω 2 .
described by the Landau-Lifshitz-Gilbert equation (LLG), which, in the so-called Gilbert form, Here, γ = 1.76 × 10 11 rad/(s T) (radian per second per tesla) and µ 0 = 4π × 10 −7 N/A 2 (newton per square ampere) are the gyromagnetic ratio and the permeability of vacuum, respectively, while α > 0 is the nondimensional empiric Gilbert damping parameter. The effective field H eff : Ω T → R 3 , measured in A/m, depends on M and is proportional to the negative functional derivative of the total magnetic Gibbs free energy with respect to M, i.e., In (2) the energy functional reads and consists of four terms, which correspond to the exchange energy, the anisotropy energy, Zeeman's energy, and the magnetostatic energy, respectively. In (3), A > 0 is the so-called exchange stiffness constant, measured in joule per meter (J/m), and K > 0 is the anisotropic constant (in J/m 3 ), while φ : S 2 → R is a (nondimensional) smooth function, which takes into account the anisotropy of the ferromagnetic material. Moreover, H e is a given external field (in A/m), while u : R 3 → R refers to the magnetostatic potential, which is the unique solution of the full-space transmission problem Combining (2) and (3), we obtain the following expression for the effective field where H s = −∇u denotes the stray field (in A/m). 4 The dynamics of the spin accumulation S : Ω T → R 3 , measured in A/m, is described by the diffusion equation where D 0 : Ω → R is the diffusion coefficient (in m 2 /s), λ s f is the characteristic length of the spin-flip relaxation, and λ J is related to the mean free path of an electron (both measured in m). The spin current J S : Ω T → R 3×3 , measured in A/s, is defined by where µ B = 9.2741 × 10 −24 A m 2 is the Bohr magneton, e = −1.602 × 10 −19 A s is the charge of the electron, and J e : Ω T → R 3 is the applied current density field (in A/m 2 ), while the constants 0 < β, β < 1 are the nondimensional spin polarization parameters of the magnetic layers. In (6) we denote by ∇S · M ∈ R 3 the matrix-vector product between the transpose of the Jacobian ∇S and M, i.e., In (5)- (6), it is implicitly assumed that M = 0 in the nonmagnetic but conducting material Ω \ ω.
To describe the dynamics of the magnetization, we take into account the interaction between the spin accumulation and the magnetization. Thus, we consider an augmented version of (1), where the constant J in N/A 2 is the strength of the interaction between the spin accumulation and the magnetization. Finally, to complete the setting, (5)- (7)

Nondimensional form of the problem
We introduce a nondimensional form of the system (5)- (7). We perform the substitution t = γµ 0 M s t, with t being the so-called (nondimensional) reduced time, and set T = γµ 0 M s T . We rescale the spatial variable by x = x/L, with L being a characteristic length of the problem (measured in m), e.g., the intrinsic length scale L = 2A/µ 0 M 2 s . However, to simplify our notation, we write t, T , x, ω, and Ω, instead of t , T , x , ω/L, and Ω/L, respectively. We introduce the nondimensional vector unknowns m = M/M s , so that the modulus constraint becomes |m| = 1, and s = S/M s . Furthermore, we set With these substitutions, the nondimensional augmented form of LLG becomes where the effective field is given by To simplify our notation and without loss of generality, we assume that L = λ s f = λ J .
Remark 4. The boundary term in (11b) is missing in [21]. This error has recently been noticed and corrected, so that the overall result of [21] remains valid [18]. The present analysis provides an alternate proof for the existence of solutions and hence validity of the results of [21,18].
The following lemma highlights the parabolic nature of equation (9b).
As a consequence, since D 0 ≥ D * and 0 < ββ < 1, we get This establishes (12) and concludes the proof.

Numerical algorithm
For the time discretization, we consider a uniform partition Given a sequence of functions ϕ j 0≤ j≤N , such that any ϕ j is associated with the time-step For the spatial discretization, let T Ω h h>0 be a shape-regular and (globally) quasi-uniform family of regular tetrahedral triangulations of Ω, parameterized by the meshsize h We assume that ω is resolved, i.e., Let us denote by S 1 (T Ω h ) 3 the standard finite element space of globally continuous and piecewise affine functions from Ω to R 3 . Correspondingly, we also consider 3 , we denote the nodal interpolation operators onto these spaces. Since ω is resolved, these operators coincide on ω, i.e., I Ω h (ϕ)| ω = I ω h (ϕ| ω ) for all ϕ ∈ C(Ω). In particular, there is no ambiguity, if we denote both operators by I h . The set of nodes of the triangulation T ω h is denoted by N ω h . We recall that, under the constraint |m| = 1, the strong form of (9a) can equivalently be stated as This formulation is used to construct the upcoming numerical scheme. Since (14) is linear in ∂ t m, the main idea is to introduce an additional free variable v = ∂ t m. To discretize v, we introduce the discrete tangent space defined by 3 . Moreover, we consider the set These sets reflect two main properties of m and v, namely the orthogonality m · v = 0 and the unit-length constraint |m| = 1.
We consider the nodal projection map h and φ h ∈ U h . A simple argument based on the elementwise use of barycentric coordinates shows that Π h φ h L ∞ (ω) = 1 for all φ h ∈ U h . Moreover, we have the estimate where the constant c Π > 0 depends only on the shape-regularity of the triangulation, cf., e.g., [11,Lemma 2.2]. With an additional angle condition on T ω h , it is well known that (15) holds even with c Π = 1, cf. [10].
Let m 0 h ∈ M h and s 0 h ∈ S 1 (T Ω h ) 3 be suitable approximations of the initial conditions. Moreover, we consider a numerical realization π h : L 2 (ω) → L 2 (ω) of π, which is assumed to fulfill a certain set of properties, see (H2)-(H3) below. This allows us to include the approximation errors, e.g., those which arise from the numerical computation of the stray field, into the overall convergence analysis. For ease of presentation, we assume that f and j are continuous in time, i.e., f ∈ C(0, T ; L 2 (ω)) and j ∈ C(0, T ; H 1 (Ω)), so that the expressions f j = f(t j ) and j j = j(t j ) 8 are meaningful for all 0 ≤ j ≤ N. It is even possible to replace f j and j j by some numerical approximation f j h and j j h as long as some weak convergence properties are fulfilled, cf. [15]. Analogously to what we have done in Section 2 for the continuous problem, for 0 ≤ i ≤ N − 1 we define the bilinear form a i+1 h : For the numerical integration of the SDLLG system (9), we propose the following algorithm.
The overall system (9) is a nonlinearly coupled system of a linear diffusion equation for s with the nonlinear LLG equation for m. However, our scheme only requires the solution of two linear systems per time-step, since the treatment of the micromagnetic part and the spin diffusion part is completely decoupled for the time-integration. This greatly simplifies an actual numerical implementation as well as the possible preconditioning of iterative solvers.
Remark 7. Unlike this work, earlier results on the tangent plane integrator [2,4,15,6,8,24] define m i+1 h := Π h (m i h +kv i h ) in (16b). Unconditional convergence in the sense of Theorem 12 can then be achieved with an additional angle condition on the triangulation T ω h , which ensures (15) with c Π = 1. This assumption is avoided in the present work.
The following result follows from standard scaling arguments.
for all w h ∈ S 1 (T h ).
The constant C > 0 depends only on r, but is independent of the meshsize h. 9 The following proposition states that the above algorithm is well defined, cf. [11, Proposition 3.1 and Proposition 4.1] for corresponding results in the frame of harmonic maps and the harmonic map heat flow.
Proposition 9. Algorithm 6 is well defined in the following sense: For each time as well as where the constant C * > 0 depends only on the shape-regularity of T ω h h>0 , but is independent of h and k.
Proof. Let 0 ≤ i ≤ N − 1. For step (i) of the algorithm, it is straightforward to show that problem (16a) is characterized by a positive definite bilinear form. Unique solvability thus follows from linearity and finite space dimension.
Step (ii) is clearly well defined. For all z ∈ N h , the nodewise orthogonality from K m i h proves This proves (17). The norm equivalence from Lemma 8 in the case r = 2 yields This establishes (18). For step (iii), we use the same argument as for step (i). Due to (17), the nodewise projections in (16c) are well defined. Let b i h : As 0 < ββ < 1 and D * > 0, b i h (·, ·) is positive definite and problem (16c) is thus well posed.

Convergence analysis
In this section, we consider the convergence properties of Algorithm 6 and show that it is indeed unconditionally convergent towards a weak solution of SDLLG in the sense of Definition 2. We emphasize that the proof is constructive in the sense that it even shows existence of weak solutions. We start by collecting some general assumptions: (H2) The general field contribution π h is bounded, i.e., with a constant C π > 0 which depends only on |ω|.
Within this setting, and with an appropriate modification of assumption (H3), the hybrid FEM-BEM method from [19] for the computation of the stray field can also be included into our analysis. Then, the proof of Proposition 19 below becomes more technical, but the assertion remains true. We refer to the argument of [15] which can be adapted accordingly.
From now on, we consider the time-approximations m hk , m ± hk , s hk , s ± hk defined by (13). The next theorem is the main result of this work. (b) In addition to the above, let assumption (H3) be satisfied. Then, it holds where (m, s) is a weak solution of SDLLG.
Remark 13. In particular, Theorem 12 yields existence of weak solutions, and each accumulation point of (m hk , s hk ) is a weak solution of SDLLG in the sense of Definition 2.
The proof of Theorem 12 will roughly be done in three steps, namely (i) boundedness of the discrete quantities and energies, (ii) existence of weakly convergent subsequences via compactness, (iii) identification of the limits with weak solutions of SDLLG.
For the sake of readability, we split our argument into several lemmata.
To start with, we recall the following result, which states a well-known and simple algebraic trick which often simplifies the computation and the estimation of sums.
Lemma 14 (Abel's summation by parts). Let X be a vector space endowed with a symmetric bilinear form (·, ·). Given an integer j ≥ 1, let {v i } 0≤i≤ j ⊂ X. Then, it holds The first ingredient for step (i) is the following proposition.
The constant C > 0 depends only on the data, but is in particular independent of the discretization parameters h and k.
Proof. Let 0 ≤ i ≤ j − 1. For (16c), we choose ζ h = s i+1 h as test function. After multiplication by k, we obtain cf. the proof of Lemma 5. Summing up over i = 0, . . . , j − 1, and exploiting Abel's summation by parts from Lemma 14 for the term j−1 h Ω , we get Exploiting 0 < 1 − ββ < 1 on the left-hand side, the Cauchy-Schwarz inequality and the Young inequality on the right-hand side, we obtain, for any choice of ε > 0, Here the constant C > 0 is the stability constant of the trace operator. It follows that If we choose ε < 2D * (1 − ββ )/Cβ, then all the coefficients on the left-hand side are positive. From (H1) and the regularity of j, we know that the right-hand side is uniformly bounded with respect to h and k. This yields the estimate (19).
Proof. The result follows from the boundedness of the discrete functions s i+1 h 0≤i≤N−1 from Proposition 15.
is quasi-uniform, it is well known that P h is stable in H 1 (Ω), i.e., We also refer to [9,23] for recent results on H 1 -stability on locally refined meshes. With this, we obtain uniform boundedness of ∂ t s hk . 13 Proposition 17. The sequence {∂ t s hk } is uniformly bounded in L 2 (0, T ; H −1 (Ω)), i.e., where the constant C > 0 depends only on the data, but is in particular independent of the discretization parameters h and k. Proof . From (16c) and the H 1 -stability (20) of P h , we get Dividing by w H 1 (Ω) and taking the supremum over w ∈ H 1 (Ω) \ {0}, we obtain Squaring, integrating over (t i , t i+1 ), and summing up over 0 ≤ i ≤ N − 1, we get The boundedness from Proposition 15 thus yields (21).
We derive the corresponding estimates for the discrete quantities Exploiting the vector identity 2 (a + b) · a = |a| 2 + |a + b| 2 − |b| 2 for all a, b ∈ R 3 with the choice a = k∇v i h and b = ∇m i h , and taking into account (16b), we obtain (22).
The constant C > 0 depends only on the data and k 0 , but is otherwise independent of the discretization parameters h and k.
Proof. Let 1 ≤ j ≤ N. From Lemma 18,multiplying (22) by k/C exch , summing up over 0 ≤ i ≤ j − 1 and exploiting the telescopic sum, we obtain The Cauchy-Schwarz inequality and the Young inequality, together with assumption (H2), yield for any ε > 0 From Proposition 9, we deduce where the constant C > 0 depends only on |ω|, T and C * . We thus obtain Note that θ > 1/2. If we choose ε < 2α/(2 + c), for k < k 0 := ε (2α − (2 + c)ε) / C 2 π C all the coefficients on the left-hand side are positive. From the regularity of f, assumption (H1), and the boundedness from Proposition 15, we know that the right-hand side is uniformly bounded. This yields the estimate (23).
Corollary 20. Under the assumptions of Proposition 19, and if k < k 0 , the sequences {m hk }, m ± hk , Π h m + hk and v − hk are uniformly bounded. In particular, it holds where the constant C > 0 depends only on the data and k 0 , but is independent of the discretization parameters h and k.
Proof. The result follows from the boundedness of the discrete functions v i h , m i+1 h 0≤i≤N−1 from Proposition 9 and Proposition 19, and from (15).
We can now proceed with step (ii) of the proof and conclude the existence of weakly convergent subsequences.
, we can identify the limits of the subsequences of {m hk } and m ± hk . As ∂ t m hk = v − hk , it clearly holds that v = ∂ t m a.e. in ω T . We now prove that the limiting function m satisfies the unit-length constraint. First, we observe that for (h, k) → 0. For 0 ≤ i ≤ N − 1 and K ∈ T ω h , a standard interpolation estimate for the piecewise linear function m i+1 From Proposition 19, we obtain For all 1 ≤ j ≤ N, Proposition 9 and the discrete norm equivalence of Lemma 8 with r = 2 yield Then, from Proposition 19, we deduce Combining (26)-(27), the triangle inequality thus yields that m + hk 2 sub −→ 1 in L 1 (ω T ) for (h, k) → 0, whence |m| = 1 a.e. in ω T follows from (25).
For x ∈ R 3 with |x| ≥ 1, it holds that Due to (17), this yields for whence by virtue of Proposition 19 k.
Finally, from Proposition 17, we deduce the existence of a weakly convergent subsequence of {∂ t s hk }, and it is easy to see that its limit is precisely ∂ t s, cf. [16, Section 7.1.2, Theorem 3]. This establishes (24e)-(24f) and thus concludes the proof.
Up to an additive constant, the nondimensional counterpart of (29) reads The following theorem proves an energy estimate which generalizes (30) to weak solutions.