Electrons in a material pose local magnetic orders, but typically do not exhibit a spontaneous macroscopic magnetic ordering unless a collective motion of these local magnetic orders is present. This results in a net magnetization even in the absence of an external magnetic field. Such a material is called a ferromagnet. It has binary stable configurations, which makes it an ideal material for data recording and storage. Recent advances in experiment and theory Žutić et al. (2004) have demonstrated effective and precise control of ferromagnetic configurations by means of external fields.
A very common phenomenological model for magnetization dynamics is the Landau-Lifshitz-Gilbert (LLG) equation Landau and Lifshitz (1935); Gilbert (1955). This model has been successfully used to interpret various experimental observations. The LLG equation is technically quasilinear, nonlocal and has a non-convex constraint, which poses interesting challenges in designing efficient numerical methodologies. In addition, the magnetization reversal process requires numerical methods to resolve different length and temporal scales in the presence of domain walls and vortices, due to their important roles in the switching process (Shi et al., 1999, 2000). Therefore, numerical methods for the LLG equation with high accuracy and efficiency are highly demanding.
There has been a continuous progress of developing numerical algorithms in the past few decades; see for example (Kruzík and Prohl, 2006; Cimrák, 2008) and references therein. The spatial derivative is typically approximated by the finite element method (FEM) (Abert et al., 2014; Bruckner et al., 2014; Baňas et al., 2014; Aurada et al., 2015; Le et al., 2015; Baňas et al., 2015; Praetorius et al., 2018; Alouges and Jaisson, 2006; Alouges, 2008; Alouges et al., 2014a, b) and the finite difference method (E and Wang, 2000; Fuwa et al., 2012; Jeong and Kim, 2010; Kim and Lipnikov, 2017; Romeo et al., 2008). As for the temporal discretization, explicit schemes (Alouges and Jaisson, 2006; Jiang et al., 2001), fully implicit schemes (Prohl, 2001; Bartels and Prohl, 2006; Fuwa et al., 2012), and semi-implicit schemes (E and Wang, 2000; Wang et al., 2001; García-Cervera and E, 2003; Lewis and Nigam, 2003; Cimrák, 2005; Gao, 2014; Boscarino et al., 2016) have been extensively studied. Explicit schemes suffer from severe stability constraints. Fully implicit schemes can overcome this severe stability constraints. However, a (nonsymmetric) nonlinear system of equations needs to be solved at each time step, which is time-consuming. A nonlinear multigrid method is used to handle the nonlinearity at each time step in (Jeong and Kim, 2014), and the fixed point iteration technique is used to deal with the nonlinearity in (Cimrák and Slodička, 2004). In (Fuwa et al., 2012), the existence and uniqueness of a solution to the nonlinear system is proved under the condition that the temporal stepsize be proportional to the square of the spatial gridsize. This, however, is contrary to the unconditional stability of implicit schemes.
Semi-implicit schemes achieve a desired balance between stability and efficiency. One of the most popular methods is the Gauss-Seidel projection method (GSPM) developed by Wang, García-Cervera, and E (E and Wang, 2000; Wang et al., 2001; García-Cervera and E, 2003). This method is based on a combination of a Gauss-Seidel implementation of a fractional step implicit solver for the gyromagnetic term, and the projection method for the heat flow of harmonic maps to overcome the difficulties associated with the stiffness and nonlinearity. Only several linear systems of equations need to be solved at each step, whose complexity is comparable to solving the scalar heat equation implicitly. It is tested that GSPM is unconditionally stable with first-order accuracy in time. In order to get second-order accuracy in time, two nonsymmetric linear systems of equations need to be solved at each step. Note that a projection step is needed to preserve the pointwise length constraint.
In this work, we propose two second-order semi-implicit projection methods for LLG equation based on the second-order backward differentiation formula and the second-order interpolation formula using the information at previous two temporal steps. The unconditional unique solvability of both methods is proved, with their second-order accuracy verified through numerical examples in both 1D and 3D. The efficiency of both methods is compared to that of another two popular schemes in the literature. In addition, we test the robustness of both methods using the first benchmark problem for a ferromagnetic thin film material developed by the micromagnetic modeling activity group from National Institute of Standards and Technology (NIST). It is worth mentioning that a modification of the proposed method has been proved to be second-order accurate in time under mild conditions Chen et al. (2019).
The rest of the paper is organized as follows. In Section 2, we first introduce the micromagnetics model based on the LLG equation. The second-order semi-implicit projection methods are described in Section 2.2 with their unique solvability given in Section 2.3. The calculation of the demagnetization field (stray filed) is given in Section 2.4. Numerical results in Section 3 are used to test the accuracy and the efficiency of the proposed methods in both 1D and 3D. Moreover, the first benchmark for a ferromagnetic thin film material developed by the micromagnetics modeling activity group from NIST is used to check the applicability of the proposed methods in Section 4. Conclusions are drawn in Section 5.
2 Second-order semi-implicit methods
2.1 Landau-Lifshitz-Gilbert equation
on . Here the magnetization
is a three-dimensional vector field with a pointwise constraintand is the unit outward normal vector. is a bounded domain occupied by the ferromagnetic material. The first term on the right hand side in eq. 1 is the gyromagnetic term and the second term is the damping term with being the dimensionless damping coefficient.
The effective field consists of the exchange field, the anisotropy field, the external field and the demagnetization or stray field . For a uniaxial material,
Here, the dimensionless parameters are and with the diameter of the ferromagnetic body, the anisotropy constant, the exchange constant, the permeability of vacuum, and the saturation magnetization. , and denotes the standard Laplacian operator. is the applied (external) magnetic field and the detailed description of will be given in §2.4. Typical values of the physical parameters for Permalloy are included in Table 1.
|Physical Parameters for Permalloy|
2.2 Second-order semi-implicit projection methods
Denote the temporal step-size by , the spatial mesh size by , the standard second-order centered difference for Laplacian operator by , and , with the final time. For convenience, we use the finite difference method to approximate the spatial derivatives in eqs. 6 and 5. For the temporal discretization, we employ the second-order backward differentiation formulas (BDFs) to approximate the temporal derivative
Such a discretization results the following fully implicit scheme for eq. 5:
As expected, at each time step, a nonlinear system of equations needs to be solved in eq. 8. Moreover, the nonsymmetric structure of the system introduces additional difficulties.
To overcome this severe difficulty, we approximate the nonlinear prefactors in front of the discrete Laplacian term using the information from previous time steps (one-sided interpolation) with its accuracy the same as the corresponding BDF scheme. For (8), we have
and . However, such a scheme cannot preserve the magnitude of magnetization, we therefore add a projection step and obtain the following scheme for eq. 5
where is the intermediate magnetization. For eq. 6, using the same idea, we have
Using the same idea, we can construct high-order semi-implicit projection schemes for eq. 5 and eq. 6 using high-order BDFs and high-order one-sided interpolations. However, if the order is greater than , then such a scheme cannot be A-stable. Numerical tests for LLG equation also indicates the conditional stability of higher order semi-implicit projection schemes. First-order semi-implicit projection schemes are A-stable, but they do not have obvious advantages over the Gauss-Seidel projection method Wang et al. (2001).
To kick start Scheme A and Scheme B in implementation, we use the first-order semi-implicit projection scheme with the first-order BDF and the first-order one-sided interpolation for the first temporal step, and thus the whole method is still second-order accurate.
2.3 Unconditionally unique solvability
For simplicity of illustration, we assume that the spatial mesh size . An extension to the general case is straightforward.
We firstly introduce the discrete inner product.
Definition 1 (Discrete inner product).
For grid functions and over the uniform numerical grid, we define
where is the index set and is the index which closely depends on the dimension .
For the grid function , we define the average of summation as
For the grid function with , we define the discrete -norm as
For ease of notation, we drop the temporal indices and rewrite Scheme A and Scheme B in a more compact form
where , , and are given.
Theorem 2.1 (Solvability for Scheme A).
Given and , the numerical scheme eq. 12 admits a unique solution.
and take the average on both sides. Therefore, can be represented as follows:
To proceed, we need the following lemma.
Let be a real, reflexive Banach space and let (the dual space of ) be bounded, continuous, coercive (i.e., , as ) and monotone. Then for any there exists a solution of the equation .
Furthermore, if the operator is strictly monotone, then the solution is unique.
Proof of Theorem 2.1.
For any , with , we denote
and derive the following monotonicity estimate:
Note that the following equality and inequality have been applied in the second step:
The third step is based on the fact that both and are constants, and , so that . Moreover, for any , with , we get
and the equality only holds when .
Therefore, an application of Lemma 2.2 implies a unique solution of Scheme A.
Theorem 2.3 (Solvability for Scheme B).
Given and , the numerical scheme eq. 13 admits a unique solution.
Proof of Theorem 2.3.
We first rewrite eq. 13 in a compact form , where
is the identity matrix,is the antisymmetric matrix corresponding to the discrete operator , and is the symmetric positive definite matrix corresponding to which admits a decomposition with being nonsingular.
Thus, we have
Due to the antisymmetry of matrix , we have
It follows from the spectral lemma for antisymmetric matrices (Prasolov, 1994)
that the eigenvalues ofare either or purely imaginary, thus the eigenvalues of are either or purely imaginary as well. This unique solvability comes as a consequence of the fact that all eigenvalues of have as real parts and . ∎
Note that the unique solvability of Scheme A and Scheme B does not impose any condition on and . This is in contrast with earlier results for the fully implicit schemes where is needed for the unique solution of the nonlinear system of equations at each time step; see (Fuwa et al., 2012) for example.
2.4 Computation of the stray field
The stray field with the scalar function in which satisfies
together with jump conditions
Here denotes the jump of at the material boundary as
and thus the stray field
where is the Newtonian potential.
It follows from eq. 18 combined with the divergence theorem that
The evaluation of the stray field can be carried out by performing an integration over the entire material for every point . In terms of the computational complexity, a direct evaluation of (18) requires with
the degree of freedoms. Moreover, we need to evaluate (18) at each time step. Therefore, the direct evaluation is computationally expensive and thus a fast solver is highly desirable. The complexity for solving stray field using FFT is Wang et al. (2001); García-Cervera and E (2003).
3 Accuracy and efficiency test
We use examples in both 1D and 3D to show the second-order accuracy of Scheme A and Scheme B.
In order to have the exact magnetization profile, we consider the simplified LLG equation with only the exchange term and set with the forcing term
where with the exact solution. The ferromagetic body in 1D and in 3D. The final time . Since the exchange term is the stiffest term in the original LLG equation, it is adequate to use (20) to test accuracy and efficiency of the proposed methods.
The exact solution in 1D is
which satisfies the homogeneous Neumann boundary condition.
The exact solution in 3D is
where , , and .
Since both schemes are semi-implicit, we compare their efficiency with another two semi-implicit methods: Gauss-Seidel projection method Wang et al. (2001) and the second-order implicit-explicit method (Boscarino et al., 2016). For completeness, we first state these two methods.
3.1 The Gauss-Seidel projection method
Heat flow without constraints :
Projection onto :
3.2 The second-order implicit-explicit method
For simplicity, we only formulate the second-order implicit-explicit (IMEX2) scheme (Boscarino et al., 2016) without Gilbert damping. Define
The second-order IMEX scheme (IMEX2) reads as
where . Note that IMEX2 solves two linear systems of equations at each step.
The computational cost for GSPM, IMEX2, and BDF2 at each temporal step is as follows. Seven symmetric linear systems of equations with constant coefficients and dimension need to be solved for GSPM, and two nonsymmetric linear systems of equations with variable coefficients and dimension need to be solved for IMEX2. BDF2 only needs to solve one nonsymmetric linear systems of equations with variable coefficients and dimension for both Scheme A and Scheme B. Here is the number of unknowns in each spatial dimension.
3.3 Accuracy test
Since Scheme A and Scheme B are comparable numerically, we only show results of Scheme A, termed as BDF2. For comparison, we also list results of the other two semi-implicit methods: GSPM and IMEX2.
In 1D, for (21), we fix and record the temporal error in terms of the temporal stepsize in Table 2 and Figure (a)a to get the temporal accuracy. Both BDF2 and IMEX2 are second-order accurate, while GSPM is first-order accurate. To get the spatial accuracy, we fix and record the spatial error in terms of in Table 3 and Figure (b)b. BDF2, IMEX2, and GSPM are all second-order accurate.
In 3D, for (22), we fix and record the temporal error in terms of in table 4 and fig. (a)a to get the temporal accuracy. Both BDF2 and IMEX2 are second-order accurate, while GSPM is first-order accurate. To get the spatial accuracy, we fix the temporal stepsize and record the spatial error in terms of in table 5 and fig. (b)b. BDF2, IMEX2, and GSPM are all second-order accurate.
3.4 Efficiency comparison
In 1D, for the same tolerance, costs of BDF2, GSPM, and IMEX2 in fig. (a)a satisfy: BDF2IMEX2GSPM. In 3D, for the same tolerance, costs of BDF2, GSPM, and IMEX2 in fig. (b)b satisfy: BDF2IMEX2GSPM. For both cases, BDF2 is slightly better than IMEX2 since two linear systems of equations need to be solved in IMEX2 while only one linear system needs to solved in BDF2. Both BDF2 and IMEX2 are better than GSPM.
When the spatial mesh is very fine, we observe that to achieve the same tolerance, costs of BDF2, GSPM, and IMEX2 satisfy: GSPMBDF2IMEX2. The reason is that fast solvers can solve symmetric linear systems with constant coefficients in GSPM, while nonsymmetric linear systems with variable coefficients are involved in both BDF2 and IMEX2. It becomes increasingly difficult to solve such systems using the Generalized Minimum Residual Method (GMRES), for example. This issue will be further explored in a subsequent work.
4 Benchmark problem from NIST
To examine our methods in the realistic case, we simulate the first standard problem established by the micromagnetic modelling activity group at National Institute of Standards and Technology (NIST) Micromagnetic Modeling Activity Group (2000). This problem asks for simulating the hysteresis loop of a thin-film element with material parameters that are not too different from Permalloy. The hysteresis loop is obtained in the following way: A positive external field of strength , in the unit of is applied. The magnetization is able to reach a steady state. Once this steady state is approached, the applied external field is reduced by a certain amount, and the material sample is again allowed to reach a steady state. The process continues until the hysteresis system attains a negative field of strength . The process then is repeated, increasing the field in small steps until it reaches the initial applied external field. As a consequence, we are able to plot the average magnetization at the steady state as a function of the external filed strength during the hysteresis loop. For BDF2, the unsymmetric linear systems of equations are solved by GMRES from library called high performance preconditioners Center for Applied Scientific Computing, Lawrence Livermore National Laboratory (1952) which was developed by Center for Applied Scientific Computing, Lawrence Livermore National Laboratory.
4.1 Magnetization profile
For comparison, we use the same setup of the available code mo96a of the first standard problem from NIST. Its setup is grid points and the canting angle of applied field from nominal axis. The calculation of the demagnetization field is done by FFT. The initial state is uniform. In the loop, successive steps between and are adopted for both x-loop and y-loop.
Due to the presence of meta-stable symmetric states, the applied fields should be rotated one degree counterclockwise off the nominal axis. The damping coefficient , the temporal stepsize and the cell size is . A stopping criterion is used to determine that a steady state is reached when the relative change in the total energy is less than . For mo96a, LABEL:NIST_long_magnetization and LABEL:NIST_short_magnetization plot the average remanent magnetization on the bottom surface of the sample in the xy- plane when when the applied fields are approximately parallel (canting angle ) to the y- (long) axis and the x- (short) axis, respectively. For BDF2, the corresponding results are shown in LABEL:BDF2_long_magnetization and LABEL:BDF2_short_magnetization, respectively. The in-plane magnetization components are represented by arrows in LABEL:NIST_BDF2_magnetization.
Furthermore, magnetization components are visualized by the grayscale value in fig. 5. For mo96a, fig. (a)a - fig. (d)d plot the x- component and the y- component when the applied field is along the y- axis, the x- component and the y- component when the applied field is along the x- axis, respectively. Results of BDF2 are shown in Figures (h)h, (g)g, (f)f and (e)e.
From LABEL:NIST_BDF2_magnetization and fig. 5, we observe that results of BDF2 are in qualitative agreements with those of mo96a.
4.2 Hysteresis loop
Hysteresis loops generated by the code mo96a are shown in fig. (a)a when the applied field is approximately parallel to the long axis and in fig. (b)b when the applied field is approximately parallel to the short axis, respectively. The average remanent magnetization in reduced units is for the y-loop and for the x-loop. The coercive fields are in fig. (a)a and in fig. (b)b.
For BDF2, hysteresis loops are presented in fig. (c)c when the applied field is approximately parallel to the long axis and in fig. (d)d when the applied field is approximately parallel to the short axis, respectively. The average remanent magnetization in reduced units is for the y-loop and for the x-loop. The coercive fields are in fig. (c)c and in fig. (d)d.
Based on these results, we conclude that results of BDF2 agree well with those of NIST, both qualitatively and quantitatively.
In this paper, we proposed two second-order semi-implicit projection methods to solve the Landau-Lifshitz-Gilbert equation, which possess unconditionally unique solvability at each time step. Examples in both 1D and 3D are used to verify the accuracy and the efficiency of both schemes. Additionally, the first benchmark problem from NIST is used to check the applicability of our methods in the realistic situation using the full Landau-Lifshitiz-Gilbert equation. Results show that our schemes can produce the correct hysteresis loop with quantities of interest agreeing with other methods in a quantitative manner.
One issue associated to the proposed methods is that nonsymmetric linear systems with variable coefficients are involved at each time step. It becomes increasingly difficult to solve such systems using the Generalized Minimum Residual Method. However, such linear systems have some unique features, as a consequence of the unique structure at the continuous level. This shall be used to develop more efficient linear solvers. Meanwhile, the technique presented here may be applicable to the model for current-driven domain wall dynamics Chen et al. (2015) and the Schrödinger-Landau-Lifshitz system Chen et al. (2016), which shall be explored later.
This work is supported in part by the grants NSFC 21602149 (J. Chen), AFORS grant FA9550-18-1-0095 (C. J. García-Cervera), NSF DMS-1418689 (C. Wang), the Postgraduate Research and Practice Innovation Program of Jiangsu Province via grant KYCX19_1947 (C. Xie), and NSFC 11801016 (Z. Zhou).
- Žutić et al. (2004) I. Žutić, J. Fabian, S. Das Sarma, Spintronics: Fundamentals and applications, Rev. Mod. Phys. 76 (2004) 323–410.
- Landau and Lifshitz (1935) L. Landau, E. Lifshitz, On the theory of the dispersion of magnetic permeability in ferromagnetic bodies, Phys. Z. Sowjet. 8 (1935) 153–169.
- Gilbert (1955) T. Gilbert, Phys. Rev. 100 (1955) 1243. [Abstract only; full report, Armor Research Foundation Project No. A059, Supplementary Report, May 1, 1956 (unpublished)].
- Shi et al. (1999) J. Shi, S. Tehrani, M. R. Scheinfein, Magnetization vortices and anomalous switching in patterned NiFeCo submicron arrays, Appl. Phys. Lett. 74 (1999) 2525–2527.
- Shi et al. (2000) J. Shi, T. Zhu, S. Tehrani, Y. Zheng, J.-G. Zhu, Geometry dependence of magnetization vortices in patterned submicron NiFe elements, Appl. Phys. Lett. 76 (2000) 2588–2590.
- Kruzík and Prohl (2006) M. Kruzík, A. Prohl, Recent developments in the modeling, analysis, and numerics of ferromagnetism, SIAM Rev. 48 (2006) 439–483.
- Cimrák (2008) I. Cimrák, A survey on the numerics and computations for the Landau-Lifshitz equation of micromagnetism, Arch. Comput. Method Eng. 15 (2008) 277–309.
- Abert et al. (2014) C. Abert, G. Hrkac, M. Page, D. Praetorius, M. Ruggeri, D. Suess, Spin-polarized transport in ferromagnetic multilayers: An unconditionally convergent FEM integrator, Comput. Math. Appl. 68 (2014) 639–654.
- Bruckner et al. (2014) F. Bruckner, D. Suess, M. Feischl, T. Führer, P. Goldenits, M. Page, D. Praetorius, M. Ruggeri, Multiscale modeling in micromagnetics: Existence of solutions and numerical integration, Math. Models Methods Appl. Sci. 24 (2014) 2627–2662.
- Baňas et al. (2014) L. Baňas, M. Page, D. Praetorius, J. Rochat, A decoupled and unconditionally convergent linear FEM integrator for the Landau–Lifshitz–Gilbert equation with magnetostriction, IMA J. Numer. Anal. 34 (2014) 1361–1385.
- Aurada et al. (2015) M. Aurada, J. Melenk, D. Praetorius, FEM–BEM coupling for the large-body limit in micromagnetics, J. Comput. Appl. Math. 281 (2015) 10–31.
- Le et al. (2015) K.-N. Le, M. Page, D. Praetorius, T. Tran, On a decoupled linear FEM integrator for eddy-current-LLG, Appl. Anal. 94 (2015) 1051–1067.
- Baňas et al. (2015) L. Baňas, M. Page, D. Praetorius, A convergent linear finite element scheme for the Maxwell-Landau-Lifshitz-Gilbert equation, Electron. Trans. Numer. Anal. 44 (2015) 250–270.
- Praetorius et al. (2018) D. Praetorius, M. Ruggeri, B. Stiftner, Convergence of an implicit–explicit midpoint scheme for computational micromagnetics, Comput. Math. Appl. 75 (2018) 1719–1738.
- Alouges and Jaisson (2006) F. Alouges, P. Jaisson, Convergence of a finite element discretization for the Landau-Lifshitz equations in micromagnetism, Math. Models Methods Appl. Sci. 16 (2006) 299–316.
- Alouges (2008) F. Alouges, A new finite element scheme for Landau-Lifchitz equations, Discrete Contin. Dyn. Syst. Ser. S 1 (2008) 187–196.
- Alouges et al. (2014a) F. Alouges, E. Kritsikis, J. Steiner, J.-C. Toussaint, A convergent and precise finite element scheme for Landau–Lifschitz–Gilbert equation, Numer. Math. 128 (2014a) 407–430.
- Alouges et al. (2014b) F. Alouges, A. de Bouard, A. Hocquet, A semi-discrete scheme for the stochastic Landau–Lifshitz equation, Stoch. Partial Differ. Equ. Anal. Comput. 2 (2014b) 281–315.
- E and Wang (2000) W. E, X. Wang, Numerical methods for the Landau-Lifshitz equation, SIAM J. Numer. Anal. 38 (2000) 1647–1665.
- Fuwa et al. (2012) A. Fuwa, T. Ishiwata, M. Tsutsumi, Finite difference scheme for the Landau-Lifshitz equation, Japan J. Indust. Appl. Math. 29 (2012) 83–110.
- Jeong and Kim (2010) D. Jeong, J. Kim, A Crank–Nicolson scheme for the Landau–Lifshitz equation without damping, J. Comput. Appl. Math. 234 (2010) 613–623.
- Kim and Lipnikov (2017) E. Kim, K. Lipnikov, The mimetic finite difference method for the Landau–Lifshitz equation, J. Comput. Phys. 328 (2017) 109–130.
- Romeo et al. (2008) A. Romeo, G. Finocchio, M. Carpentieri, L. Torres, G. Consolo, B. Azzerboni, A numerical solution of the magnetization reversal modeling in a permalloy thin film using fifth order Runge–Kutta method with adaptive step size control, Phys. B: Condens. Matter 403 (2008) 464 – 468.
- Jiang et al. (2001) J. S. Jiang, H. G. Kaper, G. K. Leaf, Hysteresis in layered spring magnets, Discrete Continuous Dyn. Syst. Ser. B. 1 (2001) 219–232.
- Prohl (2001) A. Prohl, Computational micromagnetism, volume 1, Springer, 2001.
- Bartels and Prohl (2006) S. Bartels, A. Prohl, Convergence of an implicit finite element method for the Landau-Lifshitz-Gilbert equation, SIAM J. Numer. Anal. 44 (2006) 1405–1419.
- Wang et al. (2001) X. Wang, C. J. García-Cervera, W. E, A Gauss-Seidel projection method for micromagnetics simulations, J. Comput. Phys. 171 (2001) 357–372.