# Convergence analysis of a second-order accurate, linear numerical scheme for
the Landau-Lifshitz equation with large damping parameters

Yongyong Cai  Laboratory of Mathematics and Complex Systems and School of
Mathematical Sciences
Beijing Normal University
Beijing
China<EMAIL_ADDRESS>,  Jingrun Chen  School of Mathematical
Sciences
Soochow University
Suzhou
China<EMAIL_ADDRESS>(Corresponding author) ,  Cheng Wang
Mathematics Department
University of Massachusetts
North Dartmouth
MA 02747
USA<EMAIL_ADDRESS>and  Changjian Xie  Department of Mathematics
The Hong Kong University of Science and Technology
Clear Water Bay
Kowloon
Hong Kong
China<EMAIL_ADDRESS>
###### Abstract.

A second order accurate, linear numerical method is analyzed for the Landau-
Lifshitz equation with large damping parameters. This equation describes the
dynamics of magnetization, with a non-convexity constraint of unit length of
the magnetization. The numerical method is based on the second-order backward
differentiation formula in time, combined with an implicit treatment of the
linear diffusion term and explicit extrapolation for the nonlinear terms.
Afterward, a projection step is applied to normalize the numerical solution at
a point-wise level. This numerical scheme has shown extensive advantages in
the practical computations for the physical model with large damping
parameters, which comes from the fact that only a linear system with constant
coefficients (independent of both time and the updated magnetization) needs to
be solved at each time step, and has greatly improved the numerical
efficiency. Meanwhile, a theoretical analysis for this linear numerical scheme
has not been available. In this paper, we provide a rigorous error estimate of
the numerical scheme, in the discrete
$\ell^{\infty}(0,T;\ell^{2})\cap\ell^{2}(0,T;H_{h}^{1})$ norm, under suitable
regularity assumptions and reasonable ratio between the time step-size and the
spatial mesh-size. In particular, the projection operation is nonlinear, and a
stability estimate for the projection step turns out to be highly challenging.
Such a stability estimate is derived in details, which will play an essential
role in the convergence analysis for the numerical scheme, if the damping
parameter is greater than 3.

###### Key words and phrases:

Landau-Lifshitz equation, large damping parameters, linear numerical scheme,
second order accuracy, convergence analysis, stability estimate for the
projection step

###### 2010 Mathematics Subject Classification:

35K61, 65M06, 65M12

##  1\. Introduction

The Landau-Lifshitz (LL) equation [Landau1935On], with quasilinearity and the
constraint of unit length of magnetization, describes the evolution of the
magnetization in ferromagnetic materials with applications of information
storage in the magnetic-based recording devices [Prohl2001Computational]. The
nonlinear conservative term of the LL equation preserves the unit length of
magnetization and drives the system. The remaining nonlinear part related to
the harmonic mapping of the LL equation is dissipated by a factor of damping
parameters. Such a parameter plays an important role for energy evolution,
which can be calculated [Miranda2021MechanismsBL].

There have been extensive numerical works for the LL equation
[GilbertKelly1955, LiZengPRLGiant2019, Miranda2021MechanismsBL,
Tanaka2014MicrowaveAssistedMR]. One of the most popular temporal
discretization is the semi-implicit method [jingrun2019analysis,
cimrak2005error, gao2014optimal], which turns out to remarkably relax
restrictions of temporal step-size. For example, the linearly implicit
backward Euler approach has been well studied in [Lubich2021,
alouges2006convergence, jingrun2019analysis, gao2014optimal]. A combination of
this numerical idea with a high-order non-conforming finite element
discretization in space has been proposed and analyzed in
[alouges2006convergence], in which a projection is applied to an approximate
tangent space to the normality constraint. Moreover, a convergence analysis in
both space and time has been established in [Lubich2021], by evaluating the
approximated error of time derivative term which is orthogonal to the
magnetization. The error estimates for linearly implicit schemes, based on
either backward Euler or Crank-Nicolson method, combined with finite
element/finite difference spatial discretization, have been obtained in
[An2021, An2016OptimalEE, gao2014optimal]. The backward differentiation
formula (BDF)-based linearly implicit methods have been analyzed in
[Lubich2021, jingrun2019analysis], and a second-order accuracy have been
rigorously proved under the same condition that temporal step-size
proportional to the spatial mesh-size in both space and time.

Meanwhile, it is noticed that, all these existing numerical methods lead to an
asymmetric linear system of equations with the coefficients dependent of the
updated magnetization. An efficient numerical solver for such an asymmetric
linear system is highly non-trivial, which usually results in a very expensive
computation cost in the three-dimensional simulation. Therefore, a numerical
scheme only involved with a constant coefficient linear system, so that the
coefficients are independent of the updated magnetization, is highly
desirable. In fact, such a linear numerical scheme has been proposed and
studied in a recent work [Xie2021BDF2Linear]. In more details, the second-
order BDF stencil is used in the temporal discretization, the perfect Laplace
term (in the harmonic mapping part) is treated implicitly, while the nonlinear
terms are approximated by explicit extrapolation formulas. After an
intermediate magnetization vector is obtained by this linear algorithm, a
projection of magnetization onto the unit sphere is applied, to satisfy the
non-convex constraint of unit length. Of course, this numerical approach leads
to a linear system with constant coefficients independent the updated
magnetization at each time step. Because of this subtle fact, the linear
numerical scheme has demonstrated great advantages in the simulation of
ferromagnetic materials for large damping parameters [Xie2021BDF2Linear] .
Furthermore, extensive simulation experiments have indicated that the proposed
linear numerical method preserves better stability property as the damping
parameter takes large values, in comparison with all the existing works
[Lubich2021, alouges2006convergence, An2016OptimalEE, Boscarino2016High,
jingrun2019analysis], etc.

On the other hand, a theoretical analysis of the proposed linear numerical
scheme has not been available, in spite of its extensive advantages in the
practical computations for large damping parameters. The key theoretical
difficulty is associated with the fact that, a fully explicit treatment of the
nonlinear gyromagnetic term (by an extrapolation formula) breaks its
(energetic) conservative feature at the numerical level, so that a direct
control of this nonlinear term becomes a very challenging issue. The only
hopeful approach is to control this term by the linear diffusion term in the
harmonic mapping part, while the fact that the nonlinear gyromagnetic term and
the linear diffusion term are updated by different temporal discretization
makes this estimate highly non-trivial. In this article, we provide the
convergence analysis and the optimal rate error estimate for the proposed
linear numerical scheme, in the discrete
$\ell^{\infty}(0,T;\ell^{2})\cap\ell^{2}(0,T;H_{h}^{1})$ norm, if the damping
parameter is greater than 3. To overcome the above-mentioned difficulties, we
build the stability estimate of the projection step, which will play a crucial
role in the rigorous error estimate for the original error function. In
particular, a standard $\ell^{2}$ stability estimate for the projection step
is not sufficient to recover the convergence analysis, and an $H_{h}^{1}$
stability estimate turns out to be necessary at the projection step, which
comes from the technique of controlling the nonlinear gyromagnetic term by the
linear diffusion term. In more details, the a-priori $W_{h}^{1,\infty}$
estimate for the numerical solution and a-priori $H_{h}^{1}$ estimate for the
intermediate numerical error function at the previous time step is needed in
the error analysis. Meanwhile, the a-priori $H_{h}^{1}$ estimate for the
numerical error for the magnetization vector can be controlled by a growth
factor $1+\delta$ acting on the $H_{h}^{1}$ estimate of the intermediate
magnetization error function, with $\delta$ being arbitrary positive number.
Such a $W_{h}^{1,\infty}$ bound for numerical solution and $H_{h}^{1}$
estimate for the intermediate numerical error function can be recovered at the
next time step, with the help of the inverse inequality and a mild temporal
constraint and large damping parameter (larger than three). In turn, an
estimate for numerical error function becomes a straightforward consequence of
an application of discrete Gronwall inequality, combined with the fine
estimate of a growth factor $1+\delta$ acting on the $H_{h}^{1}$ estimate of
the intermediate error function.

The rest of this paper is organized as follows. In Section 2, we review the
fully discrete numerical schemes and state the main theoretical result of
convergence. The detailed proof is provided in Section 3. Finally, some
concluding remarks are made in Section 4.

##  2\. The mathematical model and the numerical scheme

###  2.1. The Landau-Lifshitz equation

The LL equation is formulated as

(2.1)
$\begin{cases}\partial_{t}{\mbox{\boldmath$m$}}(\mbox{\boldmath$x$},t)=-{\mbox{\boldmath$m$}}\times\Delta{\mbox{\boldmath$m$}}+\alpha\Delta\mbox{\boldmath$m$}+\alpha|\nabla\mbox{\boldmath$m$}|^{2}\mbox{\boldmath$m$},&\mbox{\boldmath$x$}\in\Omega,\,t>0,\\\
\partial_{\bf
n}\mbox{\boldmath$m$}(\mbox{\boldmath$x$},t)|_{\partial\Omega}=0,&\mbox{\boldmath$x$}\in\partial\Omega,\,t\geq
0,\\\
{\mbox{\boldmath$m$}}(\mbox{\boldmath$x$},0)=\mbox{\boldmath$m$}_{0}(\mbox{\boldmath$x$}),&\mbox{\boldmath$x$}\in\Omega,\end{cases}$

where $x$ and $t$ are the variables of space and time, respectively,
$\Omega\subset\mathbb{R}^{d}$ ($d=1,2,3$, with $d$ being the spatial
dimension) is a bounded domain and ${\bf n}$ is the unit outward normal vector
along $\partial\Omega$,
${\mbox{\boldmath$m$}}:=\mbox{\boldmath$m$}(\mbox{\boldmath$x$},t)=(m_{1},m_{2},m_{3})^{T}\,:\,\Omega\subset\mathbb{R}^{d}\to
S^{2}$ represents the magnetization vector field with
$|{\mbox{\boldmath$m$}}|=1,\;\forall\mbox{\boldmath$x$}\in\Omega$, and
$\alpha>0$ is the damping parameter. The notations $\partial_{t}$, $\nabla$,
and $\Delta$ represent the temporal derivative, the gradient and the
Laplacian, respectively. The homogeneous Neumann boundary condition is
considered. The first term on the right hand side of eq. 2.1 is the
gyromagnetic term, and the remaining term related to $\alpha$ is the damping
term. In comparison with the ferromagnetic model [Landau1935On], eq. 2.1 only
includes the exchange term which poses the main difficulties in numerical
analysis, as done in the literature [bartels2006convergence,
weinan2001numerical, gao2014optimal]. To simplify the presentation, we set
$\Omega=[0,1]^{d}$ and consider the 3D case in this paper, while the results
hold for the 1D and 2D models.

###  2.2. Finite difference discretization and the numerical method

We set the temporal step-size as $k>0$, so that the time step instant becomes
$t^{n}=nk$ ($n\leq\left\lfloor\frac{T}{k}\right\rfloor$, with $T$ being the
final time, $\left\lfloor\cdot\right\rfloor$ being the floor operator). The
spatial mesh-size is given by $h_{x}=\frac{1}{N_{x}}$,
$h_{y}=\frac{1}{N_{y}}$, $h_{z}=\frac{1}{N_{z}}$, with $N_{x}$, $N_{y}$ and
$N_{z}$ being the number of grid points of uniform partition along $x$, $y$
and $z$ directions, respectively. We use the half grid points
$({x}_{i-\frac{1}{2}},{y}_{j-\frac{1}{2}},{z}_{\ell-\frac{1}{2}})$ (also
written as $(\hat{x}_{i},\hat{y}_{j},\hat{z}_{\ell})$), with
$x_{i-\frac{1}{2}}=(i-\frac{1}{2})h_{x}$,
$y_{j-\frac{1}{2}}=(j-\frac{1}{2})h_{y}$ and
$z_{\ell-\frac{1}{2}}=(\ell-\frac{1}{2})h_{z}$ ( $i=0,1,\cdots,N_{x}+1$;
$j=0,1,\cdots,N_{y}+1$; $\ell=0,1,\cdots,N_{z}+1$). The numerical domain
becomes
$\Omega_{h}=\\{(\hat{x}_{i},\hat{y}_{j},\hat{z}_{\ell})|i=0,1,\cdots,N_{x}+1;j=0,1,\cdots,N_{y}+1;\ell=0,1,\cdots,N_{z}+1\\}$
and the interior domain is
$\Omega_{h}^{0}=\\{(\hat{x}_{i},\hat{y}_{j},\hat{z}_{\ell})|i=1,\cdots,N_{x};j=1,\cdots,N_{y};\ell=1,\cdots,N_{z}\\}$,
and $\Omega_{h}/\Omega_{h}^{0}$ is the set of boundary (ghost) points. We
introduce the notation of the discrete vector grid function
$\mbox{\boldmath$m$}_{h}(\mbox{\boldmath$x$})\in\mathbb{R}^{3}$ defined for
$\mbox{\boldmath$x$}\in\Omega_{h}$ with
$\mbox{\boldmath$m$}_{i,j,\ell}=\mbox{\boldmath$m$}_{h}(\hat{x}_{i},,\hat{y}_{j},\hat{z}_{\ell})$
(similar notations for the scalar functions ), and the discrete homogeneous
Neumann boundary condition reads for $i_{x}=0,N_{x}$, $j_{y}=0,N_{y}$,
$\ell_{z}=0,N_{z}$ and $0\leq i\leq N_{x}+1,0\leq j\leq N_{y}+1,0\leq\ell\leq
N_{z}+1$,

(2.2)
$\mbox{\boldmath$m$}_{i_{x},j,\ell}=\mbox{\boldmath$m$}_{i_{x}+1,j,\ell},\quad\mbox{\boldmath$m$}_{i,j_{y},\ell}=\mbox{\boldmath$m$}_{i,j_{y}+1,\ell},\quad\mbox{\boldmath$m$}_{i,j,\ell_{z}}=\mbox{\boldmath$m$}_{i,j,\ell_{z}+1}.$

Let
$X=\\{f_{h}(\mbox{\boldmath$x$})\in\mathbb{R},\mbox{\boldmath$x$}\in\Omega_{h},\,f_{h}\,\mbox{satisfies
boundary condition}~{}\eqref{BC-1}\\}$ be the scalar function space and
$\mbox{\boldmath$X$}=\\{\mbox{\boldmath$m$}_{h}(\mbox{\boldmath$x$})\in\mathbb{R}^{3},\mbox{\boldmath$x$}\in\Omega_{h},\,\mbox{\boldmath$m$}_{h}\,\mbox{satisfies
boundary condition}~{}\eqref{BC-1}\\}$ be the vector-valued function space.
The corresponding continuous version is denoted by
$X_{e},\mbox{\boldmath$X$}_{e}$. The standard second-order centered difference
approximation for the Laplace operator results in

$\Delta_{h}\mbox{\boldmath$m$}_{i,j,\ell}=\delta_{x}^{2}\mbox{\boldmath$m$}_{i,j,\ell}+\delta_{y}^{2}\mbox{\boldmath$m$}_{i,j,\ell}+\delta_{z}^{2}\mbox{\boldmath$m$}_{i,j,\ell},\quad\delta_{x}^{2}\mbox{\boldmath$m$}_{i,j,\ell}=\frac{\mbox{\boldmath$m$}_{i+1,j,\ell}-2\mbox{\boldmath$m$}_{i,j,\ell}+\mbox{\boldmath$m$}_{i-1,j,\ell}}{h_{x}^{2}},$

where $\delta_{y}^{2}$, $\delta_{z}^{2}$ for the approximation of
$\partial_{yy}$, $\partial_{zz}$ could be similarly defined. The forward
finite difference operators $D_{x},D_{y}$ and $D_{z}$ are defined for
$f_{h}\in X$:

$D_{x}f_{i,j,\ell}=\frac{f_{i+1,j,\ell}-f_{i,j,\ell}}{h_{x}},\,D_{y}f_{i,j,\ell}=\frac{f_{i,j+1,\ell}-f_{i,j,\ell}}{h_{y}},\,D_{z}f_{i,j,\ell}=\frac{f_{i,j,\ell+1}-f_{i,j,\ell}}{h_{z}}.$

These finite difference operators could be applied to the scalar or vector
grid functions in the same way. The discrete gradient operator (forward)
$\nabla_{h}\mbox{\boldmath$m$}_{h}$ with
$\mbox{\boldmath$m$}_{h}=(u_{h},v_{h},w_{h})^{T}\in\mbox{\boldmath$X$}$ reads
as

$\displaystyle\nabla_{h}\mbox{\boldmath$m$}_{i,j,\ell}=\begin{bmatrix}D_{x}u_{i,j,\ell}&D_{x}v_{i,j,\ell}&D_{x}w_{i,j,\ell}\\\
D_{y}u_{i,j,\ell}&D_{y}v_{i,j,\ell}&D_{y}w_{i,j,\ell}\\\
D_{z}u_{i,j,\ell}&D_{z}v_{i,j,\ell}&D_{z}w_{i,j,\ell}\end{bmatrix}.$

A semi-implicit numerical scheme has been proposed in [Xie2018], and used in
the numerical simulation for small damping parameter models. In more details,
semi-implicit approximations are applied to the two nonlinear terms, namely
$\mbox{\boldmath$m$}\times\Delta\mbox{\boldmath$m$}$ and
$\mbox{\boldmath$m$}\times(\mbox{\boldmath$m$}\times\Delta\mbox{\boldmath$m$})$,
in which $\Delta\mbox{\boldmath$m$}$ is treated implicitly, while the
coefficient variables are explicitly updated via an extrapolation formula. The
theoretical convergence analysis has been established in a more recent work
[jingrun2019analysis]. However, this numerical scheme involves a large linear
system with time-dependent coefficients, related to the updated magnetization
at each time step, and the symmetry is not available in the linear system, due
to the nonlinear structure. To overcome this subtle difficulty, which leads to
significant computational costs (especially in the 3D case), we make use of
the alternate PDE formulation (2.1), and treat the linear diffusion term
$\alpha\Delta\mbox{\boldmath$m$}$ implicitly, while the two nonlinear terms,
namely $-\mbox{\boldmath$m$}\times\Delta\mbox{\boldmath$m$}$ and
$\alpha|\nabla\mbox{\boldmath$m$}|^{2}\mbox{\boldmath$m$}$, are discretized in
a fully explicit way. Subsequently, a point-wise projection is applied to the
intermediate field, so that the numerical solution of $m$ has a unit length at
the point-wise level. In more details, the following numerical scheme has been
studied, namely Algorithm 2.1: for given
$\mbox{\boldmath$m$}_{h},\tilde{\mbox{\boldmath$m$}}_{h}^{n},\mbox{\boldmath$m$}_{h}^{n+1},\tilde{\mbox{\boldmath$m$}}_{h}^{n+1}\in\mbox{\boldmath$X$}$,
find
$\tilde{\mbox{\boldmath$m$}}_{h}^{n+2},\mbox{\boldmath$m$}_{h}^{n+2}\in\mbox{\boldmath$X$}$
by solving

$\displaystyle\frac{\frac{3}{2}\tilde{\mbox{\boldmath$m$}}_{h}^{n+2}(\mbox{\boldmath$x$})-2\tilde{\mbox{\boldmath$m$}}_{h}^{n+1}(\mbox{\boldmath$x$})+\frac{1}{2}\tilde{\mbox{\boldmath$m$}}_{h}^{n}(\mbox{\boldmath$x$})}{{k}}$
$\displaystyle=-\hat{\mbox{\boldmath$m$}}_{h}^{n+2}\times\Delta_{h}\hat{\tilde{\mbox{\boldmath$m$}}}_{h}^{n+2}$
$\displaystyle+\alpha\Delta_{h}\tilde{\mbox{\boldmath$m$}}_{h}^{n+2}+\alpha|{\mathcal{A}}_{h}\nabla_{h}\hat{\mbox{\boldmath$m$}}_{h}^{n+2}|^{2}\hat{\mbox{\boldmath$m$}}_{h}^{n+2},\quad\mbox{\boldmath$x$}\in\Omega_{h}^{0},$

followed by the point-wise projection
$\mbox{\boldmath$m$}_{h}^{n+2}=\frac{\tilde{\mbox{\boldmath$m$}}_{h}^{n+2}}{|\tilde{\mbox{\boldmath$m$}}_{h}^{n+2}|}$,
where the extrapolation formula is defined by
$\hat{\mbox{\boldmath$m$}}_{h}^{n+2}=2\mbox{\boldmath$m$}_{h}^{n+1}-\mbox{\boldmath$m$}_{h}^{n},\,\,\,\hat{\tilde{\mbox{\boldmath$m$}}}_{h}^{n+2}=2\tilde{\mbox{\boldmath$m$}}_{h}^{n+1}-\tilde{\mbox{\boldmath$m$}}_{h}^{n}$,
and ${\mathcal{A}}_{h}\nabla_{h}$ (second approximation to the gradient
operator) is an average gradient operator defined for the gird function
$\mbox{\boldmath$m$}_{h}=(u_{h},v_{h},w_{h})^{T}\in\mbox{\boldmath$X$}$ as
$\mathcal{A}_{h}\nabla_{h}\mbox{\boldmath$m$}_{h}=\nabla_{h}\mathcal{A}_{h}\mbox{\boldmath$m$}_{h}$
and
$\mathcal{A}_{h}\mbox{\boldmath$m$}_{h}=(\mathcal{A}_{x}u_{h},\mathcal{A}_{y}v_{h},\mathcal{A}_{z}w_{h})$:

$\mathcal{A}_{x}u_{i,j,\ell}=\frac{u_{i,j,\ell}+u_{i-1,j,\ell}}{2},\,\mathcal{A}_{y}v_{i,j,\ell}=\frac{v_{i,j,\ell}+v_{i,j-1,\ell}}{2},\,\mathcal{A}_{z}w_{i,j,\ell}=\frac{w_{i,j,\ell}+w_{i,j,\ell-1}}{2}.$

In addition, a modified version is proposed by [Xie2021BDF2Linear], namely
Algorithm 2.2: for given
$\mbox{\boldmath$m$}_{h},\mbox{\boldmath$m$}_{h}^{n+1}\in\mbox{\boldmath$X$}$,
denote
$\hat{\mbox{\boldmath$m$}}_{h}^{n+2}=2\mbox{\boldmath$m$}_{h}^{n+1}-\mbox{\boldmath$m$}_{h}^{n}$,
and find
$\mbox{\boldmath$m$}_{h}^{n+2},\tilde{\mbox{\boldmath$m$}}_{h}^{n+2}\in\mbox{\boldmath$X$}$
by solving

$\displaystyle\frac{\frac{3}{2}\tilde{\mbox{\boldmath$m$}}_{h}^{n+2}(\mbox{\boldmath$x$})-2\mbox{\boldmath$m$}_{h}^{n+1}(\mbox{\boldmath$x$})+\frac{1}{2}\mbox{\boldmath$m$}_{h}^{n}(\mbox{\boldmath$x$})}{{k}}$
$\displaystyle=-\hat{\mbox{\boldmath$m$}}_{h}^{n+2}\times\Delta_{h}\hat{\mbox{\boldmath$m$}}_{h}^{n+2}+\alpha\Delta_{h}\tilde{\mbox{\boldmath$m$}}_{h}^{n+2}$
(2.3)
$\displaystyle\quad+\alpha|{\mathcal{A}}_{h}\nabla_{h}\hat{\mbox{\boldmath$m$}}_{h}^{n+2}|^{2}\hat{\mbox{\boldmath$m$}}_{h}^{n+2},\qquad\mbox{\boldmath$x$}\in\Omega_{h}^{0},$
(2.4) $\displaystyle\mbox{\boldmath$m$}_{h}^{n+2}$
$\displaystyle=\frac{\tilde{\mbox{\boldmath$m$}}_{h}^{n+2}}{|\tilde{\mbox{\boldmath$m$}}_{h}^{n+2}|}.$

######  Remark 2.1.

The initial data is given by
$\mbox{\boldmath$m$}_{h}^{0}=\mathcal{P}_{h}\mbox{\boldmath$m$}_{0}\in X$,
where ${\mathcal{P}}_{h}:[C(\Omega)]^{3}\to X$ is the point-wise interpolation
as

(2.5)
$\mathcal{P}_{h}\mbox{\boldmath$m$}_{0}(\mbox{\boldmath$x$})=\mbox{\boldmath$m$}_{0}(\mbox{\boldmath$x$}),\quad\mbox{\boldmath$x$}\in\Omega_{h}^{0}.$

In addition, the first-order semi-implicit projection scheme could be applied
to obtain $\mbox{\boldmath$m$}_{h}^{1}$, so that the two-step numerical method
could be jump started. Such a single-step first order algorithm will preserve
the overall second-order accuracy in time; see the detailed analysis in the
related works [guo16, guo2021] for Cahn-Hilliard equation, in which a single
step, first order semi-implicit algorithm creates a second order accurate
numerical solution in the first step.

######  Remark 2.2.

The primary difference between Algorithm 2.1 and Algorithm 2.2 is focused on
the temporal-derivative approximation, where Algorithm 2.1 uses the
intermediate approximate magnetization, while Algorithm 2.2 takes the previous
projected values. Extensive numerical experiments have demonstrated that,
Algorithm 2.2 provides a much better stability than Algorithm 2.1 in the
simulation of the realistic ferromagnetic material with large damping
parameters, as reported in [Xie2021BDF2Linear]. In this article, we present a
theoretical justification of the stability and convergence analysis for
Algorithm 2.2.

###  2.3. Main theoretical results

For simplicity of presentation, we make an assumption that
$N_{x}=N_{y}=N_{z}=N$ (with $d=3$) so that $h_{x}=h_{y}=h_{z}=h$. An extension
to the general case is straightforward. For the gird function
$\mbox{\boldmath$f$}_{h},\mbox{\boldmath$g$}_{h}\in\mbox{\boldmath$X$}$, the
discrete $\ell^{2}$ inner product $\langle\cdot\rangle$, discrete
$\|\cdot\|_{2}$ norm and discrete $\|\cdot\|_{\infty}$ norm are defined as

$\langle{\bm{f}}_{h},{\bm{g}}_{h}\rangle=h^{3}\sum_{\mbox{\boldmath$x$}\in\Omega_{h}^{0}}\mbox{\boldmath$f$}_{h}(\mbox{\boldmath$x$})\cdot\mbox{\boldmath$g$}_{h}(\mbox{\boldmath$x$}),\,\|{\bm{f}}_{h}\|_{2}=\sqrt{\langle{\bm{f}}_{h},{\bm{f}}_{h}\rangle},\,\|\mbox{\boldmath$f$}_{h}\|_{\infty}=\max_{\mbox{\boldmath$x$}\in\Omega_{h}^{0}}|\mbox{\boldmath$f$}_{h}(\mbox{\boldmath$x$})|.$

In addition, the discrete $H_{h}^{1}$-norm is given by
$\|\mbox{\boldmath$f$}_{h}\|_{H_{h}^{1}}^{2}:=\|\mbox{\boldmath$f$}_{h}\|_{2}^{2}+\|\nabla_{h}\mbox{\boldmath$f$}_{h}\|_{2}^{2}$,
and the discrete $\ell^{p}$ ($1<p<\infty$) norm is defined as
$\|\mbox{\boldmath$f$}_{h}\|_{p}^{p}=h^{3}\sum_{\mbox{\boldmath$x$}\in\Omega_{h}^{0}}|\mbox{\boldmath$f$}_{h}(\mbox{\boldmath$x$})|^{p}$.
Such norms induce the discrete spaces

$\displaystyle\ell^{p}$
$\displaystyle=\\{\mbox{\boldmath$f$}_{h}\in\mbox{\boldmath$X$}\;|\;\|\mbox{\boldmath$f$}_{h}\|_{p}<\infty\\},\quad
0<p\leq\infty$ $\displaystyle H_{h}^{1}$
$\displaystyle=\\{\mbox{\boldmath$f$}_{h}\in\mbox{\boldmath$X$}\;|\;\|\mbox{\boldmath$f$}_{h}\|_{H_{h}^{1}}<\infty\\},$
$\displaystyle W^{1,\infty}_{h}$
$\displaystyle=\\{\mbox{\boldmath$f$}_{h}\in\mbox{\boldmath$X$}\;|\;\|\mbox{\boldmath$f$}_{h}\|_{\infty}+\|\nabla_{h}\mbox{\boldmath$f$}_{h}\|_{\infty}<\infty\\},$
$\displaystyle W^{1,4}_{h}$
$\displaystyle=\\{\mbox{\boldmath$f$}_{h}\in\mbox{\boldmath$X$}\;|\;\|\mbox{\boldmath$f$}_{h}\|_{4}+\|\nabla_{h}\mbox{\boldmath$f$}_{h}\|_{4}<\infty\\},$
$\displaystyle\ell^{\infty}(0,T;\ell^{2})$
$\displaystyle=\\{\mbox{\boldmath$f$}_{h}^{n}\in\mbox{\boldmath$X$}\;|\;\max_{n}\|\mbox{\boldmath$f$}^{n}_{h}\|_{2}<\infty,\;n\in\left[0,\left\lfloor\frac{T}{k}\right\rfloor\right]\\},$
$\displaystyle\ell^{2}(0,T;H_{h}^{1})$
$\displaystyle=\left\\{\mbox{\boldmath$f$}_{h}^{n}\in\mbox{\boldmath$X$}\;\left|\right.\;\left(\sum_{n=0}^{N}\|\mbox{\boldmath$f$}^{n}_{h}\|_{2}^{2}\right)^{\frac{1}{2}}<\infty,\;n\in\left[0,\left\lfloor\frac{T}{k}\right\rfloor\right]\right\\}.$

Meanwhile, we define the continuous spaces for the function
$\mbox{\boldmath$f$}(\mbox{\boldmath$x$},t)=(f_{1},f_{2},f_{3})$ as below,

$\displaystyle C^{3}([0,T];[C^{0}(\Omega)]^{3})$
$\displaystyle:=\\{\mbox{\boldmath$f$}(\mbox{\boldmath$x$},t)\in\mbox{\boldmath$X$}_{e}\;|\;\frac{d^{3}}{dt^{3}}f_{i}\in
C^{0}([0,T]),f_{i}\in C^{0}(\Omega),\;i=1,2,3\\},$ $\displaystyle
C^{2}([0,T];[C^{2}(\bar{\Omega})]^{3})$
$\displaystyle:=\\{\mbox{\boldmath$f$}(\mbox{\boldmath$x$},t)\in\mbox{\boldmath$X$}_{e}\;|\;\frac{d^{2}}{dt^{2}}f_{i}\in
C^{0}([0,T]),\frac{d^{2}}{d\mbox{\boldmath$x$}^{2}}f_{i}\in
C^{0}(\bar{\Omega}),\;i=1,2,3\\},$ $\displaystyle
L^{\infty}([0,T];[C^{4}(\bar{\Omega})]^{3})$
$\displaystyle:=\\{\mbox{\boldmath$f$}(\mbox{\boldmath$x$},t)\in\mbox{\boldmath$X$}_{e}\;|\;\textrm{ess
sup}_{t\in[0,T]}\frac{d^{4}}{d\mbox{\boldmath$x$}^{4}}f_{i}\in
C^{0}(\bar{\Omega}),\;i=1,2,3\\},$

where $C^{0}(\Omega)$ is the space of continuous function.

The unique solvability of scheme eq. 2.3-eq. 2.4 follows from the equivalent
form of eq. 2.3:

$(\frac{3}{2{k}}I-\alpha\Delta_{h})\tilde{\mbox{\boldmath$m$}}_{h}^{n+2}(\mbox{\boldmath$x$})=\mbox{\boldmath$q$}_{h}^{n+2}(\mbox{\boldmath$x$}),\quad\mbox{\boldmath$x$}\in\Omega_{h}^{0},$

where $\tilde{\mbox{\boldmath$m$}}_{h}^{n+2}\in\mbox{\boldmath$X$}$ and
$\mbox{\boldmath$q$}_{h}^{n+2}:=\frac{2\mbox{\boldmath$m$}_{h}^{n+1}-\frac{1}{2}\mbox{\boldmath$m$}_{h}^{n}}{{k}}-\hat{\mbox{\boldmath$m$}}_{h}^{n+2}\times\Delta_{h}\hat{\mbox{\boldmath$m$}}_{h}^{n+2}+\alpha|{\mathcal{A}}_{h}\nabla_{h}\hat{\mbox{\boldmath$m$}}_{h}^{n+1}|^{2}\hat{\mbox{\boldmath$m$}}^{n+2}$.
The left hand side corresponds to a positive-definite symmetric matrix, and
the unique solvability of the proposed scheme eq. 2.3-eq. 2.4, as well as the
Algorithm 2.1, is obvious. With the fast discrete Cosine transform, the above
linear system can be very efficiently solved.

The theoretical results concerning the convergence analysis is stated below.

######  Theorem 2.1.

Assume that the exact solution of eq. 2.1 has the regularity
$\mbox{\boldmath$m$}_{e}\in C^{3}([0,T];[C^{0}(\bar{}\Omega)]^{3})\cap
C^{2}([0,T];[C^{2}(\bar{\Omega})]^{3})\cap
L^{\infty}([0,T];[C^{4}(\bar{\Omega})]^{3})$. Denote
${\mbox{\boldmath$m$}}_{h}^{n}$ ($n\geq 0$) as the numerical solution obtained
from eq. 2.3-eq. 2.4 with the initial error satisfying
$\|\mathcal{P}_{h}\mbox{\boldmath$m$}_{e}(\cdot,t_{p})-\mbox{\boldmath$m$}_{h}^{p}\|_{2}+\|\nabla_{h}(\mathcal{P}_{h}\mbox{\boldmath$m$}_{e}(\cdot,t_{p})-\mbox{\boldmath$m$}_{h}^{p})\|_{2}=\mathcal{O}({k}^{2}+h^{2}),\,p=0,1$.
In addition, we assume the technical assumption $\alpha>3$, and
$\mathcal{C}_{1}h\leq{k}\leq\mathcal{C}_{2}h$, with
$\mathcal{C}_{1},\mathcal{C}_{2}$ being the positive constants. Then the
following convergence result holds for $2\leq
n\leq\left\lfloor\frac{T}{k}\right\rfloor$ as $h,{k}\to 0^{+}$:

(2.6)
$\displaystyle\|\mathcal{P}_{h}\mbox{\boldmath$m$}_{e}(\cdot,t_{n})-\mbox{\boldmath$m$}_{h}^{n}\|_{2}+\Big{(}{k}\sum_{p=1}^{n}\|\nabla_{h}(\mathcal{P}_{h}\mbox{\boldmath$m$}_{e}(\cdot,t_{p})-\tilde{\mbox{\boldmath$m$}}_{h}^{p})\|_{2}^{2}\Big{)}^{\frac{1}{2}}$
$\displaystyle\leq\mathcal{C}({k}^{2}+h^{2}),$

in which the constant $\mathcal{C}>0$ is independent of ${k}$ and $h$.

###  2.4. A few preliminary estimates

In this section, some preliminary inequalities are derived, which will be
useful in the error analysis presented in the next section. In addition, we
have to build a stability estimate of the projection step in the numerical
algorithm.

The proof of the standard inverse inequality and discrete Gronwall inequality
can be obtained in existing textbooks and references; we just cite the results
here. In the sequel, for simplicity of notation, we will use the uniform
constant $\mathcal{C}$ to denote all the controllable constants.

######  Lemma 2.1.

(Inverse inequality) [chen16, chenW20a, Ciarlet1978] . For each vector-valued
grid function $\mbox{\boldmath$e$}_{h}\in X$, we have

(2.7)
$\displaystyle\|{\mbox{\boldmath$e$}}_{h}^{n}\|_{\infty}\leq\gamma{h}^{-1/2}(\|{\mbox{\boldmath$e$}}_{h}^{n}\|_{2}+\|\nabla_{h}\mbox{\boldmath$e$}_{h}^{n}\|_{2}),\quad\|\nabla_{h}{\mbox{\boldmath$e$}}_{h}^{n}\|_{4}\leq\gamma{h}^{-3/4}\|\nabla_{h}{\mbox{\boldmath$e$}}_{h}^{n}\|_{2},$

in which constant $\gamma$ depends on $\Omega$, as well as the form of the
discrete $\|\cdot\|_{2}$ norm.

The following estimate will be utilized in the convergence analysis. A rough
version has been provided in a recent article [jingrun2019analysis]; here we
give a further refined estimate.

######  Lemma 2.2 (Discrete gradient acting on cross product).

For grid functions
$\mbox{\boldmath$f$}_{h},\mbox{\boldmath$g$}_{h},\mbox{\boldmath$F$}_{h}\in\mbox{\boldmath$X$}$,
we have for any $\delta>0$

(2.8)
$\displaystyle\left\langle\mbox{\boldmath$f$}_{h}\times\Delta_{h}\mbox{\boldmath$g$}_{h},\mbox{\boldmath$F$}_{h}\right\rangle$
$\displaystyle=\left\langle\mbox{\boldmath$F$}_{h}\times\mbox{\boldmath$f$}_{h},\Delta_{h}\mbox{\boldmath$g$}_{h}\right\rangle,$
(2.9)
$\displaystyle\|\nabla_{h}({\mbox{\boldmath$f$}}_{h}\times{\mbox{\boldmath$g$}}_{h})\|_{2}^{2}$
$\displaystyle\leq(1+\delta)\|{\mbox{\boldmath$f$}}_{h}\|_{\infty}^{2}\cdot\|\nabla_{h}{\mbox{\boldmath$g$}}_{h}\|_{2}^{2}+{\mathcal{C}}_{\delta}\|{\mbox{\boldmath$g$}}_{h}\|_{4}^{2}\cdot\|\nabla_{h}\mbox{\boldmath$f$}_{h}\|_{4}^{2}.$

###### Proof.

Equality (2.8) has been proved in [jingrun2019analysis], so that we only focus
on the proof of (2.9).

At each numerical mesh cell, from $(\hat{x}_{i},\hat{y}_{j},\hat{z}_{\ell})$
to $(\hat{x}_{i+1},\hat{y}_{j},\hat{z}_{\ell})$, the following expansion
identity is valid:

(2.10) $\displaystyle
D_{x}(\mbox{\boldmath$f$}_{h}\times\mbox{\boldmath$g$}_{h})_{i,j,\ell}=$
$\displaystyle(A_{x}\mbox{\boldmath$f$}_{h})_{i,j,\ell}\times(D_{x}\mbox{\boldmath$g$}_{h})_{i,j,\ell}+(A_{x}\mbox{\boldmath$g$}_{h})_{i,j,\ell}\times(D_{x}\mbox{\boldmath$f$}_{h})_{i,j,\ell},$
$\displaystyle\mbox{with}\,\,\,(A_{x}\mbox{\boldmath$f$}_{h})_{i,j,\ell}=$
$\displaystyle\frac{1}{2}((\mbox{\boldmath$f$}_{h})_{i,j,\ell}+(\mbox{\boldmath$f$}_{h})_{i+1,j,\ell}).$

In turn, we get the following expansion, over each mesh cell:

(2.11)
$D_{x}(\mbox{\boldmath$f$}_{h}\times\mbox{\boldmath$g$}_{h})=(A_{x}\mbox{\boldmath$f$}_{h})\times(D_{x}\mbox{\boldmath$g$}_{h})+(D_{x}\mbox{\boldmath$f$}_{h})\times(A_{x}\mbox{\boldmath$g$}_{h}).$

Subsequently, a careful application of discrete Hölder inequality reveals that

(2.12)
$\displaystyle\|(A_{x}\mbox{\boldmath$f$}_{h})\times(D_{x}\mbox{\boldmath$g$}_{h})\|_{2}\leq$
$\displaystyle\|A_{x}\mbox{\boldmath$f$}_{h}\|_{\infty}\cdot\|D_{x}\mbox{\boldmath$g$}_{h}\|_{2}\leq\|\mbox{\boldmath$f$}_{h}\|_{\infty}\cdot\|D_{x}\mbox{\boldmath$g$}_{h}\|_{2},$
(2.13)
$\displaystyle\|(D_{x}\mbox{\boldmath$f$}_{h})\times(A_{x}\mbox{\boldmath$g$}_{h})\|_{2}\leq$
$\displaystyle\|D_{x}\mbox{\boldmath$f$}_{h}\|_{4}\cdot\|A_{x}\mbox{\boldmath$g$}_{h}\|_{4}\leq\|D_{x}\mbox{\boldmath$f$}_{h}\|_{4}\cdot\|\mbox{\boldmath$g$}_{h}\|_{4},$

in which the fact that
$\|A_{x}\mbox{\boldmath$f$}_{h}\|_{\infty}\leq\|\mbox{\boldmath$f$}_{h}\|_{\infty}$,
$\|A_{x}\mbox{\boldmath$g$}_{h}\|_{4}\leq\|\mbox{\boldmath$g$}_{h}\|_{4}$, has
been applied. Then we get

(2.14)
$\|D_{x}(\mbox{\boldmath$f$}_{h}\times\mbox{\boldmath$g$}_{h})\|_{2}\leq\|\mbox{\boldmath$f$}_{h}\|_{\infty}\cdot\|D_{x}\mbox{\boldmath$g$}_{h}\|_{2}+\|D_{x}\mbox{\boldmath$f$}_{h}\|_{4}\cdot\|\mbox{\boldmath$g$}_{h}\|_{4}.$

The corresponding estimates in the $y$ and $z$ directions can be similarly
derived, and the technical details are skipped for the sake of brevity. A
combination of (2.14) and its counterparts in $y$ and $z$ directions leads to

(2.15)
$\displaystyle\|\nabla_{h}(\mbox{\boldmath$f$}_{h}\times\mbox{\boldmath$g$}_{h})\|_{2}^{2}\leq$
$\displaystyle\|\mbox{\boldmath$f$}_{h}\|_{\infty}^{2}\cdot\|\nabla_{h}\mbox{\boldmath$g$}_{h}\|_{2}+3\|\nabla_{h}\mbox{\boldmath$f$}_{h}\|_{4}^{2}\cdot\|\mbox{\boldmath$g$}_{h}\|_{4}^{2}$
$\displaystyle+6\|\mbox{\boldmath$f$}_{h}\|_{\infty}\cdot\|\nabla_{h}\mbox{\boldmath$f$}_{h}\|_{4}\cdot\|\mbox{\boldmath$g$}_{h}\|_{4}$
$\displaystyle\leq$
$\displaystyle(1+\delta)\|\mbox{\boldmath$f$}_{h}\|_{\infty}^{2}\cdot\|\nabla_{h}\mbox{\boldmath$g$}_{h}\|_{2}+(3+9\delta^{-1})\|\nabla_{h}\mbox{\boldmath$f$}_{h}\|_{4}^{2}\cdot\|\mbox{\boldmath$g$}_{h}\|_{4}^{2},$

for any $\delta>0$, in which the Cauchy inequality has been applied in the
last step. Therefore, the nonlinear cross product estimate (2.9) has been
derived. This finishes the proof of Lemma 2.2. ∎

The following discrete Sobolev inequality has been derived in the existing
works [guan17a, guan14a], for the discrete grid function with periodic
boundary condition; an extension to the discrete homogeneous Neumann boundary
condition can be made in a similar manner.

######  Lemma 2.3 (Discrete Sobolev inequality).

[guan17a, guan14a] For a grid function
$\mbox{\boldmath$f$}_{h}\in\mbox{\boldmath$X$}$, we have the following
discrete Sobolev inequality:

(2.16)
$\displaystyle\|\mbox{\boldmath$f$}_{h}\|_{4}\leq\mathcal{C}\|\mbox{\boldmath$f$}_{h}\|_{2}^{\frac{1}{4}}\cdot\|\mbox{\boldmath$f$}_{h}\|_{H_{h}^{1}}^{\frac{3}{4}}\leq\mathcal{C}(\|\mbox{\boldmath$f$}_{h}\|_{2}+\|\mbox{\boldmath$f$}_{h}\|_{2}^{\frac{1}{4}}\cdot\|\nabla_{h}\mbox{\boldmath$f$}_{h}\|_{2}^{\frac{3}{4}}),$

in which the positive constant $\mathcal{C}$ only depends on the domain
$\Omega$.

The following stability estimates of the point-wise projection (2.4) are
crucial for the error analysis, and the proof could be found in Appendix A.

######  Lemma 2.4.

Assume the continuous vector function
$\mbox{\boldmath$m$}_{e}\in[C(\Omega)]^{3}$ satisfies a regularity requirement
$\|\mbox{\boldmath$m$}_{e}\|_{W^{1,\infty}}\leq C^{*}$ (with $C^{*}$ being a
positive constant) and the point-wise constraint
$|\mbox{\boldmath$m$}_{e}|=1$. Denoting
$\underline{\bm{m}}_{h}={\mathcal{P}}_{h}\mbox{\boldmath$m$}_{e}\in X$, for
any grid function $\tilde{\mbox{\boldmath$m$}}_{h}\in X$, we define the
projected grid function $\mbox{\boldmath$m$}_{h}\in X$ as
$\mbox{\boldmath$m$}_{h}=\frac{\tilde{\mbox{\boldmath$m$}}_{h}}{|\tilde{\mbox{\boldmath$m$}}_{h}|}$,
and introduce the error functions as
$\mbox{\boldmath$e$}_{h}(\mbox{\boldmath$x$})=\underline{\bm{m}}_{h}(\mbox{\boldmath$x$})-\mbox{\boldmath$m$}_{h}(\mbox{\boldmath$x$})$,
$\tilde{\mbox{\boldmath$e$}}_{h}(\mbox{\boldmath$x$})=\underline{\bm{m}}_{h}(\mbox{\boldmath$x$})-\tilde{\mbox{\boldmath$m$}}_{h}(\mbox{\boldmath$x$})$
($\mbox{\boldmath$x$}\in\Omega_{h}$). Under the a-priori assumptions on
$\tilde{\mbox{\boldmath$e$}}_{h}$ or equivalently on the profile
$\tilde{\mbox{\boldmath$m$}}_{h}$:

(2.17) $\|\tilde{\mbox{\boldmath$e$}}_{h}\|_{2}\leq
2{k}^{\frac{15}{8}},\quad\|\nabla_{h}\tilde{\mbox{\boldmath$e$}}_{h}\|_{2}\leq\frac{1}{2}{k}^{\frac{11}{8}},$

the following estimates hold for sufficiently small ${k}$ and $h$ satisfying
$C_{1}h\leq{k}\leq C_{2}h$ (with $C_{1},C_{2}$ being two positive constants)

(2.18)
$\displaystyle\|\tilde{\mbox{\boldmath$e$}}_{h}\|_{2}^{2}\geq(1-{k}^{\frac{5}{4}})\|\mbox{\boldmath$e$}_{h}\|_{2}^{2}+(1-{k}^{\frac{1}{4}})\|\tilde{\mbox{\boldmath$e$}}_{h}-\mbox{\boldmath$e$}_{h}\|_{2}^{2},$
(2.19)
$\displaystyle\|\nabla_{h}\mbox{\boldmath$e$}_{h}\|_{2}^{2}\leq(1+\delta)\|\nabla_{h}\tilde{\mbox{\boldmath$e$}}_{h}\|_{2}^{2}+\mathcal{C}_{\delta}\|\tilde{\mbox{\boldmath$e$}}_{h}\|_{2}^{2},\quad\mbox{for
any $\delta>0$ }.$

In addition, for the analysis of the BDF2 temporal stencil at the projection
stage, a further refined error estimate is needed.

######  Lemma 2.5.

Consider $\underline{\bm{m}}_{h}^{(1)},\underline{\bm{m}}_{h}^{(2)}\in X$
($|\mbox{\boldmath$m$}_{h}^{(q)}(\mbox{\boldmath$x$})|=1$,
$\mbox{\boldmath$x$}\in\Omega_{h}$, $q=1,2$) with the $W_{h}^{1,\infty}$ bound
$\|\underline{\bm{m}}_{h}^{(q)}\|_{\infty}+\|\nabla_{h}\underline{\bm{m}}_{h}^{(q)}\|_{\infty}\leq
C^{*}$ ($q=1,2$). For any grid functions
$\tilde{\mbox{\boldmath$m$}}_{h}^{(1)},\tilde{\mbox{\boldmath$m$}}_{h}^{(2)}\in
X$, we define the projected grid functions
$\mbox{\boldmath$m$}_{h}^{(q)}=\frac{\tilde{\mbox{\boldmath$m$}}_{h}^{(q)}}{|\tilde{\mbox{\boldmath$m$}}_{h}^{(q)}|}$,
and introduce the error functions as
$\mbox{\boldmath$e$}_{h}^{(q)}(\mbox{\boldmath$x$})=\underline{\bm{m}}_{h}^{(q)}(\mbox{\boldmath$x$})-\mbox{\boldmath$m$}_{h}^{(q)}(\mbox{\boldmath$x$})\in
X$,
$\tilde{\mbox{\boldmath$e$}}_{h}^{(q)}(\mbox{\boldmath$x$})=\underline{\bm{m}}_{h}^{(q)}(\mbox{\boldmath$x$})-\tilde{\mbox{\boldmath$m$}}_{h}^{(q)}(\mbox{\boldmath$x$})\in
X$, $q=1,2$. Under the a-priori assumptions for
$\tilde{\mbox{\boldmath$e$}}_{h}^{(q)}$ ($q=1,2$):

(2.20) $\|\tilde{\mbox{\boldmath$e$}}_{h}^{(q)}\|_{2}\leq
2{k}^{\frac{15}{8}},\quad\|\nabla_{h}\tilde{\mbox{\boldmath$e$}}^{(q)}\|_{2}\leq\frac{1}{2}{k}^{\frac{11}{8}},\,\,\,q=1,2,$

and the assumptions for $\underline{\bm{m}}_{h}^{(q)}$ ($q=1,2$):

(2.21)
$\|\underline{\bm{m}}_{h}^{(1)}-\underline{\bm{m}}_{h}^{(2)}\|_{\infty}\leq\frac{1}{4}{k}^{\frac{7}{8}},$

the following estimate is valid for sufficiently small ${k}$ and $h$
satisfying $C_{1}h\leq{k}\leq C_{2}h$:

(2.22)
$\Big{|}\langle\tilde{\mbox{\boldmath$e$}}_{h}^{(1)}-\mbox{\boldmath$e$}_{h}^{(1)},\mbox{\boldmath$e$}_{h}^{(2)}\rangle\Big{|}\leq{k}^{\frac{5}{4}}\|\mbox{\boldmath$e$}_{h}^{(2)}\|_{2}^{2}+{k}^{\frac{1}{4}}\|\tilde{\mbox{\boldmath$e$}}_{h}^{(1)}-\mbox{\boldmath$e$}_{h}^{(1)}\|_{2}^{2}.$

We leave the proof of Lemma 2.5 to Appendix B. Lemmas 2.4 and 2.5 essentially
establish the stability of the projection step (2.4) under the assumptions
that the previous numerical solution at $t^{n}$ and $t^{n+1}$ are sufficiently
close to the exact solution.

##  3\. The optimal rate convergence analysis: Proof of Theorem 2.1

Denote
$\underline{\mbox{\boldmath$m$}}(\mbox{\boldmath$x$},t)=\mathcal{P}_{h}\mbox{\boldmath$m$}_{e}(\mbox{\boldmath$x$},t)\in\mbox{\boldmath$X$}$
($\mbox{\boldmath$x$}\in\Omega_{h}$) and
$\underline{\mbox{\boldmath$m$}}_{h}^{n}(\mbox{\boldmath$x$})=\underline{\mbox{\boldmath$m$}}_{h}(\mbox{\boldmath$x$},t^{n})$
($n\geq 0$). Around the boundary section $z=0$, we set
$\hat{z}_{0}=-\frac{1}{2}h$, $\hat{z}_{1}=\frac{1}{2}h$, and we can extend the
profile $\underline{\mbox{\boldmath$m$}}$ to the numerical “ghost” points,
according to the extrapolation formula eq. 2.2:

(3.1)
$\underline{\mbox{\boldmath$m$}}_{i,j,0}=\underline{\mbox{\boldmath$m$}}_{i,j,1},\quad\underline{\mbox{\boldmath$m$}}_{i,j,N_{z}+1}=\underline{\mbox{\boldmath$m$}}_{i,j,N_{z}},$

and the extrapolation for other boundaries can be formulated in the same
manner. The proof of such an extrapolation yields a higher order
$\mathcal{O}(h^{5})$ approximation, instead of the standard
$\mathcal{O}(h^{3})$ accuracy has been applied in [jingrun2019analysis]. Also
see the related works [STWW2003, Wang2000, Wang2004] in the existing
literature.

Performing a careful Taylor expansion for the exact solution around the
boundary section $z=0$, combined with the mesh point values:
$\hat{z}_{0}=-\frac{1}{2}h$, $\hat{z}_{1}=\frac{1}{2}h$, we get

$\displaystyle\mbox{\boldmath$m$}_{e}(\hat{x}_{i},\hat{y}_{j},\hat{z}_{0})$
$\displaystyle=\mbox{\boldmath$m$}_{e}(\hat{x}_{i},\hat{y}_{j},\hat{z}_{1})-h\partial_{z}\mbox{\boldmath$m$}_{e}(\hat{x}_{i},\hat{y}_{j},0)-\frac{h^{3}}{24}\partial_{z}^{3}\mbox{\boldmath$m$}_{e}(\hat{x}_{i},\hat{y}_{j},0)+\mathcal{O}(h^{5})$
(3.2)
$\displaystyle=\mbox{\boldmath$m$}_{e}(\hat{x}_{i},\hat{y}_{j},\hat{z}_{1})-\frac{h^{3}}{24}\partial_{z}^{3}\mbox{\boldmath$m$}_{e}(\hat{x}_{i},\hat{y}_{j},0)+\mathcal{O}(h^{5}),$

in which the homogenous boundary condition has been applied in the second
step. It remains to determine
$\partial_{z}^{3}\mbox{\boldmath$m$}_{e}(\hat{x}_{i},\hat{y}_{j},0)$, for
which we use information from the rewritten PDE (2.1) and its derivatives.
Applying $\partial_{z}$ to the first evolutionary equation in (2.1) along the
boundary section $\Gamma_{z}:z=0$ gives

(3.3) $\displaystyle\begin{aligned} &(m_{1})_{zt}-2\alpha(m_{1}(\nabla
m_{1}\cdot\nabla(m_{1})_{z}+\nabla m_{2}\cdot\nabla(m_{2})_{z}+\nabla
m_{3}\cdot\nabla(m_{3})_{z}))\\\
&-\alpha|\nabla\mbox{\boldmath$m$}_{e}|^{2}(m_{1})_{z}-\alpha((m_{1})_{zxx}+(m_{1})_{zyy}+\partial_{z}^{3}m_{1})\\\
=&(m_{3})_{z}\Delta
m_{2}+m_{3}((m_{2})_{zxx}+(m_{2})_{zyy}+\partial_{z}^{3}m_{2})\\\
&-(m_{2})_{z}\Delta
m_{3}-m_{2}((m_{3})_{zxx}+(m_{3})_{zyy}+\partial_{z}^{3}m_{3}),\quad\mbox{on}\,\,\Gamma_{z}.\end{aligned}$

The first, third terms, and the first two parts in the fourth term on the
left-hand side of (3.3) disappear, due to the homogeneous Neumann boundary
condition for $m_{1}$. For the second term on the left hand side, we observe
that

(3.4) $\nabla
m_{1}\cdot\nabla(m_{1})_{z}=(m_{1})_{x}\cdot(m_{1})_{zx}+(m_{1})_{y}\cdot(m_{1})_{zy}+(m_{1})_{z}\cdot(m_{1})_{zz}=0,\quad\mbox{on}\,\,\Gamma_{z},$

since $(m_{1})_{z}=0$ on the boundary section. Similar derivations could be
made to the two other terms on the left hand side:

(3.5) $\nabla m_{2}\cdot\nabla(m_{2})_{z}=0,\,\,\,\nabla
m_{3}\cdot\nabla(m_{3})_{z}=0,\quad\mbox{on}\,\,\Gamma_{z}.$

Meanwhile, on the right hand side of (3.3), we see that the first and third
terms, as well as the first two parts in the second and fourth terms,
disappear, which comes from the homogeneous Neumann boundary condition for
$m_{2}$ and $m_{3}$. Then we arrive at

(3.6)
$\alpha\partial_{z}^{3}m_{1}=-m_{3}\partial_{z}^{3}m_{2}+m_{2}\partial_{z}^{3}m_{3},\quad\mbox{on}\,\,\Gamma_{z}.$

Similarly, we are able to derive the following equalities:

(3.7)
$\displaystyle\alpha\partial_{z}^{3}m_{2}=m_{1}\partial_{z}^{3}m_{3}-m_{3}\partial_{z}^{3}m_{1},\quad\mbox{on}\,\,\Gamma_{z}$
$\displaystyle\alpha\partial_{z}^{3}m_{3}=m_{2}\partial_{z}^{3}m_{1}-m_{1}\partial_{z}^{3}m_{2},\quad\mbox{on}\,\,\Gamma_{z}.$

In turn, for any $\alpha>0$, we observe that the matrix
$\left(\begin{array}[]{ccc}\alpha&m_{3}&-m_{2}\\\ -m_{3}&\alpha&m_{1}\\\
m_{2}&-m_{1}&\alpha\end{array}\right)$ has a positive determinant, so that the
linear system (3.6)-(3.7) has only one trivial solution:

(3.8)
$\partial_{z}^{3}m_{1}=\partial_{z}^{3}m_{2}=\partial_{z}^{3}m_{3}=0,\quad\mbox{on}\,\,\Gamma_{z}.$

As a result, an ${\mathcal{O}}(h^{5})$ consistency accuracy for the symmetric
extrapolation is obtained:

(3.9)
$\displaystyle\mbox{\boldmath$m$}_{e}(\hat{x}_{i},\hat{y}_{j},\hat{z}_{0})=\mbox{\boldmath$m$}_{e}(\hat{x}_{i},\hat{y}_{j},\hat{z}_{1})+\mathcal{O}(h^{5}),\;\underline{\mbox{\boldmath$m$}}(\hat{x}_{i},\hat{y}_{j},\hat{z}_{0})=\underline{\mbox{\boldmath$m$}}(\hat{x}_{i},\hat{y}_{j},\hat{z}_{1})+\mathcal{O}(h^{5}).$

In other words, the extrapolation formula eq. 3.1 is indeed
$\mathcal{O}(h^{5})$ accurate.

Subsequently, a detailed calculation of Taylor expansion, in both time and
space, leads to the following truncation error estimate:

(3.10) $\displaystyle\begin{aligned}
&\frac{\frac{3}{2}\underline{\mbox{\boldmath$m$}}_{h}^{n+2}(\mbox{\boldmath$x$})-2\underline{\mbox{\boldmath$m$}}_{h}^{n+1}(\mbox{\boldmath$x$})+\frac{1}{2}\underline{\mbox{\boldmath$m$}}_{h}^{n}(\mbox{\boldmath$x$})}{{k}}\\\
&=-\hat{\underline{\mbox{\boldmath$m$}}}_{h}^{n+2}\times\Delta_{h}\hat{\underline{\mbox{\boldmath$m$}}}_{h}^{n+2}+\tau^{n+2}+\alpha\Delta_{h}\underline{\mbox{\boldmath$m$}}_{h}^{n+2}+\alpha|{\mathcal{A}}_{h}\nabla_{h}\hat{\underline{\mbox{\boldmath$m$}}}_{h}^{n+2}|^{2}\hat{\underline{\mbox{\boldmath$m$}}}_{h}^{n+2},\,\mbox{\boldmath$x$}\in\Omega_{h}^{0},\end{aligned}$

where
$\hat{\underline{\mbox{\boldmath$m$}}}_{h}^{n+2}=2\underline{\mbox{\boldmath$m$}}_{h}^{n+1}-\underline{\mbox{\boldmath$m$}}_{h}^{n}\in\mbox{\boldmath$X$}$,
$\tau^{n+2}\in\mbox{\boldmath$X$}$ and
$\|\tau^{n+2}\|_{2}\leq\mathcal{C}({k}^{2}+h^{2})$. Introducing the numerical
error functions
$\tilde{\mbox{\boldmath$e$}}_{h}^{n}=\underline{\mbox{\boldmath$m$}}_{h}^{n}-\tilde{\mbox{\boldmath$m$}}_{h}^{n}\in\mbox{\boldmath$X$}$,
${\mbox{\boldmath$e$}}_{h}^{n}=\underline{\mbox{\boldmath$m$}}_{h}^{n}-\mbox{\boldmath$m$}_{h}^{n}\in\mbox{\boldmath$X$}$,
and subtracting eq. 2.3-eq. 2.4 from the consistency estimate eq. 3.10, we
have the error evolutionary equation at the interior points
$\mbox{\boldmath$x$}\in\Omega_{h}^{0}$, for $0\leq n\leq\lfloor
T/{k}\rfloor-2$:

(3.11) $\displaystyle\begin{aligned}
\frac{\frac{3}{2}\tilde{\mbox{\boldmath$e$}}_{h}^{n+2}-2\mbox{\boldmath$e$}_{h}^{n+1}+\frac{1}{2}\mbox{\boldmath$e$}_{h}^{n}}{{k}}&=-\hat{\mbox{\boldmath$m$}}_{h}^{n+2}\times\Delta_{h}\hat{\mbox{\boldmath$e$}}_{h}^{n+2}-\left(2{\mbox{\boldmath$e$}}_{h}^{n+1}-{\mbox{\boldmath$e$}}_{h}^{n}\right)\times\Delta_{h}\hat{\underline{\bm{m}}}_{h}^{n+2}\\\
&+\alpha\Delta_{h}\tilde{\mbox{\boldmath$e$}}_{h}^{n+2}+\alpha|{\mathcal{A}}_{h}\nabla_{h}\hat{\underline{\mbox{\boldmath$m$}}}_{h}^{n+2}|^{2}\hat{\mbox{\boldmath$e$}}_{h}^{n+2}+\tau^{n+2}\\\
&+\alpha\Big{(}{\mathcal{A}}_{h}\nabla_{h}(\hat{\underline{\mbox{\boldmath$m$}}}_{h}^{n+2}+\hat{\mbox{\boldmath$m$}}_{h}^{n+2})\cdot{\mathcal{A}}_{h}\nabla_{h}\hat{\mbox{\boldmath$e$}}_{h}^{n+2}\Big{)}\hat{\mbox{\boldmath$m$}}_{h}^{n+2}\end{aligned}$

with
$\hat{\mbox{\boldmath$e$}}_{h}^{n+2}=2\mbox{\boldmath$e$}_{h}^{n+1}-\mbox{\boldmath$e$}_{h}^{n}\in\mbox{\boldmath$X$}$.

Before proceeding into the formal error estimate, we state the bound for the
exact solution ${\underline{\bm{m}}}$ and the numerical solution
$\mbox{\boldmath$m$}_{h}$. Since the exact solution
$\mbox{\boldmath$m$}_{e}\in L^{\infty}([0,T];[C^{4}(\bar{\Omega})]^{3})$, the
following bound is available, for some positive constant $\mathcal{C}$:

(3.12)
$\displaystyle\|\nabla_{h}^{r}\underline{\bm{m}}_{h}(\cdot,t)\|_{4},\,\,\|\nabla_{h}^{r}\underline{\bm{m}}_{h}(\cdot,t)\|_{\infty}\leq\mathcal{C},\quad
r=0,1,2,3,\,0\leq t\leq T.$

In addition, we make the following a-priori assumption for the numerical error
function:

(3.13)
$\|{\mbox{\boldmath$e$}}_{h}^{n}\|_{2}\leq{k}^{\frac{31}{16}},\,\,\,\|\tilde{\mbox{\boldmath$e$}}_{h}^{n}\|_{2}\leq
2{k}^{\frac{15}{8}},\,\,\,\|\nabla_{h}\tilde{\mbox{\boldmath$e$}}_{h}^{n}\|_{2}\leq\frac{1}{2}{k}^{\frac{11}{8}},\,\,\,\mbox{$n\leq
q+1$}.$

Such an assumption will be recovered by the convergence analysis at time step
$t^{q+2}$. Based on this a-priori assumption, we see that (2.17) is satisfied,
so that we are able to apply Lemma 2.4 and the estimate (2.19) to get

(3.14)
$\displaystyle\|\nabla_{h}\mbox{\boldmath$e$}_{h}^{n}\|_{2}\leq\frac{5}{4}\|\nabla_{h}\tilde{\mbox{\boldmath$e$}}_{h}^{n}\|_{2}+\mathcal{C}\|\tilde{\mbox{\boldmath$e$}}_{h}^{n}\|_{2}\leq\frac{5}{8}{k}^{\frac{11}{8}}+{\mathcal{C}}{k}^{\frac{15}{8}}\leq{k}^{\frac{11}{8}},\,\,\,\mbox{$n\leq
q+1$}.$

In turn, an application of inverse inequality implies the $\|\cdot\|_{\infty}$
and $\|\cdot\|_{W_{h}^{1,4}}$ bounds of the numerical error function
$\mbox{\boldmath$e$}_{h}^{n}$ ($n\leq q+1$):

(3.15) $\displaystyle\|\mbox{\boldmath$e$}_{h}^{n}\|_{\infty}$
$\displaystyle\leq\frac{C\|\mbox{\boldmath$e$}_{h}^{n}\|_{2}}{h^{\frac{3}{2}}}\leq\frac{C\cdot{k}^{\frac{31}{16}}}{h^{\frac{3}{2}}}\leq
C{k}^{\frac{7}{16}}\leq{k}^{\frac{3}{8}},$
$\displaystyle\|\nabla_{h}\mbox{\boldmath$e$}_{h}^{n}\|_{4}$
$\displaystyle\leq\frac{C\|\nabla_{h}\mbox{\boldmath$e$}_{h}^{n}\|_{2}}{h^{\frac{3}{4}}}\leq\frac{C\cdot{k}^{\frac{11}{8}}}{h^{\frac{3}{4}}}\leq
C{k}^{\frac{5}{8}}\leq{k}^{\frac{1}{2}}\leq\frac{1}{3}.$

Subsequently, the triangle inequality yields the desired $W_{h}^{1,4}$ bound
for the numerical solutions $\mbox{\boldmath$m$}_{h}^{n}$ and
$\tilde{\mbox{\boldmath$m$}}_{h}^{n}$ ($n\leq q+1$):

(3.16) $\displaystyle\|\nabla_{h}{\mbox{\boldmath$m$}}_{h}^{n}\|_{4}$
$\displaystyle=\|\nabla_{h}\underline{\bm{m}}_{h}^{n}-\nabla_{h}{\mbox{\boldmath$e$}}_{h}^{n}\|_{4}\leq\|\nabla_{h}\underline{\bm{m}}_{h}^{n}\|_{4}+\|\nabla_{h}{\mbox{\boldmath$e$}}_{h}^{n}\|_{4}\leq\mathcal{C}+\frac{1}{3}.$

Furthermore, we need a sharper $\|\cdot\|_{\infty}$ bound for
$\hat{\mbox{\boldmath$m$}}_{h}^{n+2}=2\mbox{\boldmath$m$}_{h}^{n+1}-\mbox{\boldmath$m$}_{h}^{n}$,
which will be needed in the later error analysis. The following extrapolation
estimate is valid, due to the $C^{3}([0,T];[C^{0}(\Omega)]^{3})$ regularity of
the exact solution $\mbox{\boldmath$m$}(\cdot,t)$:

(3.17)
$\underline{\bm{m}}_{h}^{n+2}=2\underline{\bm{m}}_{h}^{n+1}-\underline{\bm{m}}_{h}^{n}+{\mathcal{O}}({k}^{2}).$

Meanwhile, since $|\mbox{\boldmath$m$}(\mbox{\boldmath$x$},t)|\equiv 1$
($\mbox{\boldmath$x$}\in\Omega$), we conclude that

(3.18)
$\|2\underline{\bm{m}}_{h}^{n+1}-\underline{\bm{m}}_{h}^{n}\|_{\infty}\leq
1+{\mathcal{C}}{k}^{2},\quad n\leq q+1.$

Its combination with the a-priori assumption that
$\|{\mbox{\boldmath$e$}}_{h}^{n}\|_{\infty}\leq{k}+h$, for $n\leq q+1$, (as
given by (3.13)), implies that

(3.19) $\displaystyle\|\hat{\mbox{\boldmath$m$}}_{h}^{n+2}\|_{\infty}=$
$\displaystyle\|2\mbox{\boldmath$m$}_{h}^{n+1}-\mbox{\boldmath$m$}_{h}^{n}\|_{\infty}\leq\|2\underline{\bm{m}}_{h}^{n+1}-\underline{\bm{m}}_{h}^{n}\|_{\infty}+\|2\mbox{\boldmath$e$}_{h}^{n+1}-\mbox{\boldmath$e$}_{h}^{n}\|_{\infty}$
$\displaystyle\leq$ $\displaystyle
1+{\mathcal{C}}{k}^{2}+3{k}^{\frac{3}{8}}\leq\alpha_{1}:=\Big{(}\frac{3+\alpha}{6}\Big{)}^{\frac{1}{2}},$

provided that ${k}$ and $h$ are sufficiently small. In addition, we denote
$\gamma_{0}:=\alpha-3>0$, so that $\alpha_{1}^{2}=1+\frac{1}{6}\gamma_{0}$.

Next, we perform a discrete $\ell^{2}$ error estimate at $t^{q+2}$ using the
mathematical induction. Taking a discrete inner product with the numerical
error equation eq. 3.11 by
$\tilde{\mbox{\boldmath$e$}}_{h}^{n+2}\in\mbox{\boldmath$X$}$ ($n\leq q+1$)
gives that

(3.20) $\displaystyle\begin{aligned}
R.H.S.=&\left\langle-\left(2{\mbox{\boldmath$m$}}_{h}^{n+1}-{\mbox{\boldmath$m$}}_{h}^{n}\right)\times\Delta_{h}\hat{\mbox{\boldmath$e$}}_{h}^{n+2},\tilde{\mbox{\boldmath$e$}}_{h}^{n+2}\right\rangle\\\
&-\left\langle\left(2{\mbox{\boldmath$e$}}_{h}^{n+1}-{\mbox{\boldmath$e$}}_{h}^{n}\right)\times\Delta_{h}\hat{\underline{\mbox{\boldmath$m$}}}_{h}^{n+2},\tilde{\mbox{\boldmath$e$}}_{h}^{n+2}\right\rangle+\left\langle\tau^{n+2},\tilde{\mbox{\boldmath$e$}}_{h}^{n+2}\right\rangle\\\
&-\alpha\|\nabla_{h}\tilde{\mbox{\boldmath$e$}}_{h}^{n+2}\|_{2}^{2}+\alpha\langle|{\mathcal{A}}_{h}\nabla_{h}\hat{\underline{\mbox{\boldmath$m$}}}_{h}^{n+2}|^{2}\hat{\mbox{\boldmath$e$}}_{h}^{n+2},\tilde{\mbox{\boldmath$e$}}_{h}^{n+2}\rangle\\\
&+\alpha\Big{\langle}\Big{(}{\mathcal{A}}_{h}\nabla_{h}(\hat{\underline{\mbox{\boldmath$m$}}}_{h}^{n+2}+\hat{\mbox{\boldmath$m$}}_{h}^{n+2})\cdot{\mathcal{A}}_{h}\nabla_{h}\hat{\mbox{\boldmath$e$}}_{h}^{n+2}\Big{)}\hat{\mbox{\boldmath$m$}}_{h}^{n+2},\tilde{\mbox{\boldmath$e$}}_{h}^{n+2}\Big{\rangle}\\\
=:&\tilde{I}_{1}+\tilde{I}_{2}+\tilde{I}_{3}-\alpha\|\nabla_{h}\tilde{\mbox{\boldmath$e$}}_{h}^{n+2}\|_{2}^{2}+\tilde{I}_{5}+\tilde{I}_{6}.\end{aligned}$

Then all the terms are accordingly analyzed. For the term $\tilde{I}_{1}$, a
combination of the summation by parts formula by and Cauchy inequality results
in

(3.21) $\displaystyle\begin{aligned}
\tilde{I}_{1}=&\left\langle-\hat{\mbox{\boldmath$m$}}_{h}^{n+2}\times\Delta_{h}\hat{\mbox{\boldmath$e$}}_{h}^{n+2},\tilde{\mbox{\boldmath$e$}}_{h}^{n+2}\right\rangle=\left\langle\tilde{\mbox{\boldmath$e$}}_{h}^{n+2}\times\hat{\mbox{\boldmath$m$}}_{h}^{n+2},-\Delta_{h}\hat{\mbox{\boldmath$e$}}_{h}^{n+2}\right\rangle\\\
=&\left\langle\nabla_{h}\Big{[}\tilde{\mbox{\boldmath$e$}}_{h}^{n+2}\times\hat{\mbox{\boldmath$m$}}_{h}^{n+2}\Big{]},\nabla_{h}\hat{\mbox{\boldmath$e$}}_{h}^{n+2}\right\rangle\leq\frac{3}{2}\Big{\|}\nabla_{h}\Big{(}\tilde{\mbox{\boldmath$e$}}_{h}^{n+2}\times\hat{\mbox{\boldmath$m$}}_{h}^{n+2}\Big{)}\Big{\|}_{2}^{2}+\frac{1}{6}\|\nabla_{h}\hat{\mbox{\boldmath$e$}}_{h}^{n+2}\|_{2}^{2}.\end{aligned}$

Meanwhile, an application of the cross product gradient estimate (2.9) implies
that, for any $\delta>0$, the following inequality is valid:

(3.22)
$\displaystyle\Big{\|}\nabla_{h}\Big{(}\tilde{\mbox{\boldmath$e$}}_{h}^{n+2}\times\hat{\mbox{\boldmath$m$}}_{h}^{n+2}\Big{)}\Big{\|}_{2}^{2}$
$\displaystyle\leq$
$\displaystyle(1+\delta)\|\hat{\mbox{\boldmath$m$}}_{h}^{n+2}\|_{\infty}^{2}\cdot\|\nabla_{h}\tilde{\mbox{\boldmath$e$}}_{h}^{n+2}\|_{2}^{2}+{\mathcal{C}}_{\delta}\|\tilde{\mbox{\boldmath$e$}}_{h}^{n+2}\|_{4}^{2}\cdot\|\nabla_{h}\hat{\mbox{\boldmath$m$}}_{h}^{n+2}\|_{4}^{2}$
$\displaystyle\leq$
$\displaystyle(1+\delta)\alpha_{1}^{2}\|\nabla_{h}\tilde{\mbox{\boldmath$e$}}_{h}^{n+2}\|_{2}^{2}+{\mathcal{C}}_{\delta}\|\tilde{\mbox{\boldmath$e$}}_{h}^{n+2}\|_{4}^{2},$

where $\mathcal{C}_{\delta}$ is a positive constant dependent on $\delta$, in
which the a-priori bound estimates (3.15) and (3.19) have been applied. The
term $\|\nabla_{h}\hat{\mbox{\boldmath$e$}}_{h}^{n+2}\|_{2}^{2}$ can be
analyzed as follows:

(3.23)
$\displaystyle\|\nabla_{h}\hat{\mbox{\boldmath$e$}}_{h}^{n+2}\|_{2}^{2}=$
$\displaystyle\|\nabla_{h}(2\mbox{\boldmath$e$}_{h}^{n+1}-\mbox{\boldmath$e$}_{h}^{n})\|_{2}^{2}\leq
6\|\nabla_{h}\mbox{\boldmath$e$}_{h}^{n+1}\|_{2}^{2}+3\|\nabla_{h}\mbox{\boldmath$e$}_{h}^{n}\|_{2}^{2}$
$\displaystyle\leq$
$\displaystyle(1+\delta)(6\|\nabla_{h}\tilde{\mbox{\boldmath$e$}}_{h}^{n+1}\|_{2}^{2}+3\|\nabla_{h}\tilde{\mbox{\boldmath$e$}}_{h}^{n}\|_{2}^{2})+{\mathcal{C}}_{\delta}(\|\tilde{\mbox{\boldmath$e$}}_{h}^{n+1}\|_{2}^{2}+\|\tilde{\mbox{\boldmath$e$}}_{h}^{n}\|_{2}^{2}),$

in which the estimate (2.19) (in Lemma 2.4) has been repeatedly applied, due
to the a-priori assumption (3.13). Combining (3.22), (3.23) and (3.21), we get

(3.24) $\displaystyle\begin{aligned}
\tilde{I}_{1}\leq\,&\frac{3}{2}\Big{\|}\nabla_{h}\Big{(}\tilde{\mbox{\boldmath$e$}}_{h}^{n+2}\times\hat{\mbox{\boldmath$m$}}_{h}^{n+2}\Big{)}\Big{\|}_{2}^{2}+\frac{1}{6}\|\nabla_{h}\hat{\mbox{\boldmath$e$}}_{h}^{n+2}\|_{2}^{2}\\\
\leq\,&(1+\delta)\Big{(}\frac{3}{2}\alpha_{1}^{2}\|\nabla_{h}\tilde{\mbox{\boldmath$e$}}_{h}^{n+2}\|_{2}^{2}+\|\nabla_{h}\tilde{\mbox{\boldmath$e$}}_{h}^{n+1}\|_{2}^{2}+\frac{1}{2}\|\nabla_{h}\tilde{\mbox{\boldmath$e$}}_{h}^{n}\|_{2}^{2}\Big{)}\\\
&+{\mathcal{C}}_{\delta}(\|\tilde{\mbox{\boldmath$e$}}_{h}^{n+2}\|_{4}^{2}+\|\tilde{\mbox{\boldmath$e$}}_{h}^{n+1}\|_{2}^{2}+\|\tilde{\mbox{\boldmath$e$}}_{h}^{n}\|_{2}^{2}).\end{aligned}$

For the term $\tilde{I}_{2}$, by the preliminary estimate (3.12) for the exact
solution, we have

(3.25) $\displaystyle\begin{aligned}
\tilde{I}_{2}=\,&-\left\langle\hat{\mbox{\boldmath$e$}}_{h}^{n+2}\times\Delta_{h}\hat{\underline{\bm{m}}}_{h}^{n+2},\tilde{\mbox{\boldmath$e$}}_{h}^{n+2}\right\rangle\leq\,\frac{1}{2}\big{[}\|\tilde{\mbox{\boldmath$e$}}_{h}^{n+2}\|_{2}^{2}+\|\hat{\mbox{\boldmath$e$}}_{h}^{n+2}\|_{2}^{2}\cdot\|\Delta_{h}\hat{\underline{\bm{m}}}_{h}^{n+2}\|_{\infty}^{2}\big{]}\\\
\leq\,&\mathcal{C}(\|\tilde{\mbox{\boldmath$e$}}_{h}^{n+2}\|_{2}^{2}+\|\mbox{\boldmath$e$}_{h}^{n+1}\|_{2}^{2}+\|\mbox{\boldmath$e$}_{h}^{n}\|_{2}^{2}).\end{aligned}$

For the term $\tilde{I}_{3}$, an application of Cauchy inequality gives

(3.26)
$\displaystyle\tilde{I}_{3}=\left\langle\tau^{n+2},\tilde{\mbox{\boldmath$e$}}_{h}^{n+2}\right\rangle\leq\mathcal{C}\|\tilde{\mbox{\boldmath$e$}}_{h}^{n+2}\|_{2}^{2}+\mathcal{C}({k}^{4}+h^{4}).$

In terms of $\tilde{I}_{5}$, based on the $W_{h}^{1,\infty}$ bound (3.12) for
the exact solution, an application of discrete Hölder inequality gives

(3.27)
$\displaystyle\Big{\|}|{\mathcal{A}}_{h}\nabla_{h}\hat{\underline{\mbox{\boldmath$m$}}}_{h}^{n+2}|^{2}\hat{\mbox{\boldmath$e$}}_{h}^{n+2}\Big{\|}_{2}\leq\,$
$\displaystyle\|\nabla_{h}\hat{\underline{\mbox{\boldmath$m$}}}_{h}^{n+2}\|_{\infty}^{2}\cdot\|\hat{\mbox{\boldmath$e$}}_{h}^{n+2}\|_{2}\leq{\mathcal{C}}\|\hat{\mbox{\boldmath$e$}}_{h}^{n+2}\|_{2},$

and

(3.28) $\displaystyle\begin{aligned}
\tilde{I}_{5}=\,&\alpha\langle|{\mathcal{A}}_{h}\nabla_{h}\hat{\underline{\mbox{\boldmath$m$}}}_{h}^{n+2}|^{2}\hat{\mbox{\boldmath$e$}}_{h}^{n+2},\tilde{\mbox{\boldmath$e$}}_{h}^{n+2}\rangle\leq\alpha\Big{\|}|{\mathcal{A}}_{h}\nabla_{h}\hat{\underline{\mbox{\boldmath$m$}}}_{h}^{n+2}|^{2}\hat{\mbox{\boldmath$e$}}_{h}^{n+2}\Big{\|}_{2}\cdot\|\tilde{\mbox{\boldmath$e$}}_{h}^{n+2}\|_{2}\\\
\leq\,&{\mathcal{C}}\alpha\|\hat{\mbox{\boldmath$e$}}_{h}^{n+2}\|_{2}\cdot\|\tilde{\mbox{\boldmath$e$}}_{h}^{n+2}\|_{2}\leq{\mathcal{C}}(\|\mbox{\boldmath$e$}_{h}^{n+1}\|_{2}^{2}+\|\mbox{\boldmath$e$}_{h}^{n}\|_{2}^{2}+\|\tilde{\mbox{\boldmath$e$}}_{h}^{n+2}\|_{2}^{2}).\end{aligned}$

For the term $\tilde{I}_{6}$, an application of discrete Hölder inequality
gives

(3.29) $\displaystyle\begin{aligned}
\tilde{I}_{6}=\,&\alpha\Big{\langle}\Big{(}{\mathcal{A}}_{h}\nabla_{h}(\hat{\underline{\mbox{\boldmath$m$}}}_{h}^{n+2}+\hat{\mbox{\boldmath$m$}}_{h}^{n+2})\cdot{\mathcal{A}}_{h}\nabla_{h}\hat{\mbox{\boldmath$e$}}_{h}^{n+2}\Big{)}\hat{\mbox{\boldmath$m$}}_{h}^{n+2},\tilde{\mbox{\boldmath$e$}}_{h}^{n+2}\Big{\rangle}\\\
\leq\,&\alpha\Big{(}\|\nabla_{h}\hat{\underline{\mbox{\boldmath$m$}}}_{h}^{n+2}\|_{4}+\|\nabla_{h}\hat{\mbox{\boldmath$m$}}_{h}^{n+2}\|_{4}\Big{)}\cdot\|\nabla_{h}\hat{\mbox{\boldmath$e$}}_{h}^{n+2}\|_{2}\cdot\|\hat{\mbox{\boldmath$m$}}_{h}^{n+2}\|_{\infty}\cdot\|\tilde{\mbox{\boldmath$e$}}_{h}^{n+2}\|_{4}\\\
\leq\,&{\mathcal{C}}\alpha\|\nabla_{h}\hat{\mbox{\boldmath$e$}}_{h}^{n+2}\|_{2}\cdot\|\tilde{\mbox{\boldmath$e$}}_{h}^{n+2}\|_{4}\leq\mathcal{C}\alpha^{2}\gamma_{0}^{-1}\|\tilde{\mbox{\boldmath$e$}}_{h}^{n+2}\|_{4}^{2}+\frac{\gamma_{0}}{72}\|\nabla_{h}\hat{\mbox{\boldmath$e$}}_{h}^{n+2}\|_{2}^{2}\\\
\leq\,&\mathcal{C}\|\tilde{\mbox{\boldmath$e$}}_{h}^{n+2}\|_{4}^{2}+\frac{(1+\delta)\gamma_{0}}{24}(2\|\nabla_{h}\tilde{\mbox{\boldmath$e$}}_{h}^{n+1}\|_{2}^{2}+\|\nabla_{h}\tilde{\mbox{\boldmath$e$}}_{h}^{n}\|_{2}^{2})+{\mathcal{C}}_{\delta}(\|\tilde{\mbox{\boldmath$e$}}_{h}^{n+1}\|_{2}^{2}+\|\tilde{\mbox{\boldmath$e$}}_{h}^{n}\|_{2}^{2}),\end{aligned}$

in which the $W_{h}^{1,4}$ bounds (3.12), (3.16), for the exact and numerical
solutions, as well as the preliminary error estimate (3.23), have been
applied. On the other hand, the inner product of the left hand side of eq.
3.11 with $\tilde{\mbox{\boldmath$e$}}_{h}^{n+2}$ turns out to be

$\displaystyle L.H.S.=$
$\displaystyle\frac{1}{4{k}}\big{(}\|\tilde{\mbox{\boldmath$e$}}_{h}^{n+2}\|_{2}^{2}-\|\mbox{\boldmath$e$}_{h}^{n+1}\|_{2}^{2}+\|2\tilde{\mbox{\boldmath$e$}}_{h}^{n+2}-\mbox{\boldmath$e$}_{h}^{n+1}\|_{2}^{2}-\|2\mbox{\boldmath$e$}_{h}^{n+1}-\mbox{\boldmath$e$}_{h}^{n}\|_{2}^{2}$
$\displaystyle+\|\tilde{\mbox{\boldmath$e$}}_{h}^{n+2}-2\mbox{\boldmath$e$}_{h}^{n+1}+\mbox{\boldmath$e$}_{h}^{n}\|_{2}^{2}\big{)}.$

Its combination with eqs. 3.24, 3.25, 3.26, 3.28, and 3.29 and eq. 3.20 leads
to

(3.30) $\displaystyle\begin{aligned}
&\frac{1}{4{k}}(\|\tilde{\mbox{\boldmath$e$}}_{h}^{n+2}\|_{2}^{2}-\|\mbox{\boldmath$e$}_{h}^{n+1}\|_{2}^{2}+\|2\tilde{\mbox{\boldmath$e$}}_{h}^{n+2}-\mbox{\boldmath$e$}_{h}^{n+1}\|_{2}^{2}-\|2\mbox{\boldmath$e$}_{h}^{n+1}-\mbox{\boldmath$e$}_{h}^{n}\|_{2}^{2})\\\
&+\alpha_{2}\|\nabla_{h}\tilde{\mbox{\boldmath$e$}}_{h}^{n+2}\|_{2}^{2}-(1+\delta)(1+\frac{\gamma_{0}}{12})\Big{(}\|\nabla_{h}\tilde{\mbox{\boldmath$e$}}_{h}^{n+1}\|_{2}^{2}+\frac{1}{2}\|\nabla_{h}\tilde{\mbox{\boldmath$e$}}_{h}^{n}\|_{2}^{2}\Big{)}\\\
&\leq\mathcal{C}_{\delta}\Big{(}\sum\limits_{p=n}^{n+2}\|\tilde{\mbox{\boldmath$e$}}_{h}^{p}\|_{2}^{2}+\|\mbox{\boldmath$e$}_{h}^{n+1}\|_{2}^{2}+\|\mbox{\boldmath$e$}_{h}^{n}\|_{2}^{2}+\|\tilde{\mbox{\boldmath$e$}}_{h}^{n+2}\|_{4}^{2}\Big{)}+\mathcal{C}({k}^{4}+h^{4}),\end{aligned}$

where $\alpha_{2}=(\alpha-\frac{3}{2}\alpha_{1}^{2}(1+\delta))$. Meanwhile,
for the $\|\cdot\|_{4}$ error estimate for
$\tilde{\mbox{\boldmath$e$}}_{h}^{n+2}$, an application of the discrete
Sobolev inequality (2.16) (in Lemma 2.3) gives

(3.31)
$\displaystyle{\mathcal{C}}_{\delta}\|\tilde{\mbox{\boldmath$e$}}_{h}^{n+2}\|_{4}^{2}\leq$
$\displaystyle\mathcal{C}_{\delta}(\|\tilde{\mbox{\boldmath$e$}}_{h}^{n+2}\|_{2}^{2}+\|\tilde{\mbox{\boldmath$e$}}_{h}^{n+2}\|_{2}^{\frac{1}{2}}\cdot\|\nabla_{h}\tilde{\mbox{\boldmath$e$}}_{h}^{n+2}\|_{2}^{\frac{3}{2}})$
$\displaystyle\leq$
$\displaystyle\mathcal{C}_{\delta}\|\tilde{\mbox{\boldmath$e$}}_{h}^{n+2}\|_{2}^{2}+\frac{\gamma_{0}}{12}\|\nabla_{h}\tilde{\mbox{\boldmath$e$}}_{h}^{n+2}\|_{2}^{2},$

in which the Young’s inequality has been applied. Then we get

(3.32)
$\displaystyle\frac{1}{4{k}}\big{(}\|\tilde{\mbox{\boldmath$e$}}_{h}^{n+2}\|_{2}^{2}-\|\mbox{\boldmath$e$}_{h}^{n+1}\|_{2}^{2}+\|2\tilde{\mbox{\boldmath$e$}}_{h}^{n+2}-\mbox{\boldmath$e$}_{h}^{n+1}\|_{2}^{2}-\|2\mbox{\boldmath$e$}_{h}^{n+1}-\mbox{\boldmath$e$}_{h}^{n}\|_{2}^{2})$
$\displaystyle+\Big{(}\alpha_{2}-\frac{\gamma_{0}}{12}\Big{)}\|\nabla_{h}\tilde{\mbox{\boldmath$e$}}_{h}^{n+2}\|_{2}^{2}-(1+\delta)(1+\frac{\gamma_{0}}{12})\Big{(}\|\nabla_{h}\tilde{\mbox{\boldmath$e$}}_{h}^{n+1}\|_{2}^{2}+\frac{1}{2}\|\nabla_{h}\tilde{\mbox{\boldmath$e$}}_{h}^{n}\|_{2}^{2}\Big{)}$
$\displaystyle\leq$
$\displaystyle\mathcal{C}_{\delta}\Big{(}\|\tilde{\mbox{\boldmath$e$}}_{h}^{n+2}\|_{2}^{2}+\|\tilde{\mbox{\boldmath$e$}}_{h}^{n+1}\|_{2}^{2}+\|\tilde{\mbox{\boldmath$e$}}_{h}^{n}\|_{2}^{2}+\|\mbox{\boldmath$e$}_{h}^{n+1}\|_{2}^{2}+\|\mbox{\boldmath$e$}_{h}^{n}\|_{2}^{2}\Big{)}+\mathcal{C}({k}^{4}+h^{4}).$

In particular, we observe that

(3.33)
$\displaystyle\big{(}\alpha-\frac{3}{2}\alpha_{1}^{2}(1+\delta)-\frac{\gamma_{0}}{12}\Big{)}-\frac{3}{2}(1+\delta)(1+\frac{\gamma_{0}}{12})$
$\displaystyle\geq$
$\displaystyle\alpha-3\alpha_{1}^{2}(1+\delta)\geq(3+\gamma_{0})-3(1+\frac{1}{6}\gamma_{0})(1+\delta)\geq\frac{1}{4}\gamma_{0},$

if $\delta>0$ is chosen with $(1+\frac{1}{6}\gamma_{0})(1+\delta)\leq
1+\frac{1}{4}\gamma_{0}$.

Moreover, an application of the a-priori assumption (3.13) into (3.32) yields

(3.34)
$\displaystyle(\frac{1}{4{k}}-{\mathcal{C}}_{\delta})\|\tilde{\mbox{\boldmath$e$}}_{h}^{n+2}\|_{2}^{2}+\Big{(}\frac{3}{2}(1+\delta)(1+\frac{\gamma_{0}}{12})+\frac{\gamma_{0}}{4}\Big{)}\|\nabla_{h}\tilde{\mbox{\boldmath$e$}}_{h}^{n+2}\|_{2}^{2}$
$\displaystyle\leq$
$\displaystyle(1+\delta)(1+\frac{\gamma_{0}}{12})\Big{(}\|\nabla_{h}\tilde{\mbox{\boldmath$e$}}_{h}^{n+1}\|_{2}^{2}+\frac{1}{2}\|\nabla_{h}\tilde{\mbox{\boldmath$e$}}_{h}^{n}\|_{2}^{2}\Big{)}+\frac{1}{4{k}}\big{(}\|\mbox{\boldmath$e$}_{h}^{n+1}\|_{2}^{2}+\|2\hat{\mbox{\boldmath$e$}}_{h}^{n+2}\|_{2}^{2})$
$\displaystyle+\mathcal{C}_{\delta}\Big{(}\|\tilde{\mbox{\boldmath$e$}}_{h}^{n+1}\|_{2}^{2}+\|\tilde{\mbox{\boldmath$e$}}_{h}^{n}\|_{2}^{2}+\|\mbox{\boldmath$e$}_{h}^{n+1}\|_{2}^{2}+\|\mbox{\boldmath$e$}_{h}^{n}\|_{2}^{2}\Big{)}+\mathcal{C}({k}^{4}+h^{4})$
$\displaystyle\leq$
$\displaystyle\frac{3}{2}(1+\delta)(1+\frac{\gamma_{0}}{12})\cdot\frac{1}{4}{k}^{\frac{11}{4}}+{\mathcal{C}}({k}^{\frac{23}{8}}+{k}^{\frac{15}{4}}+{k}^{4}+h^{4})$
$\displaystyle\leq$
$\displaystyle\Big{(}\frac{3}{2}(1+\delta)(1+\frac{\gamma_{0}}{12})+\frac{\gamma_{0}}{8}\Big{)}\cdot\frac{1}{4}{k}^{\frac{11}{4}},$

provided that ${k}$ and $h$ are sufficiently small, and under the linear
refinement requirement, $C_{1}h\leq{k}\leq C_{2}h$. As a matter of fact, we
can choose $\gamma_{0}>0$ and $\delta>0$ and make ${k}$ sufficiently small so
that

(3.35)
$\frac{3}{2}(1+\delta)(1+\frac{\gamma_{0}}{12})+\frac{\gamma_{0}}{8}\leq
2,\quad\frac{1}{4{k}}-{\mathcal{C}}_{\delta}\geq\frac{1}{6{k}}.$

As a result, by (3.34), we arrive at

(3.36)
$\displaystyle\frac{1}{6{k}}\|\tilde{\mbox{\boldmath$e$}}_{h}^{q+2}\|_{2}^{2}\leq\frac{1}{2}{k}^{\frac{11}{4}},\quad\|\nabla_{h}\tilde{\mbox{\boldmath$e$}}_{h}^{q+2}\|_{2}^{2}\leq\frac{1}{4}{k}^{\frac{11}{4}},$
$\displaystyle\mbox{i.e.},\quad\|\tilde{\mbox{\boldmath$e$}}_{h}^{q+2}\|_{2}\leq\sqrt{3}{k}^{\frac{15}{8}}\leq
2{k}^{\frac{15}{8}},\quad\|\nabla_{h}\tilde{\mbox{\boldmath$e$}}_{h}^{q+2}\|_{2}\leq\frac{1}{2}{k}^{\frac{11}{8}},$

so that the a-priori assumption is valid for
$\tilde{\mbox{\boldmath$m$}}^{q+2}$.

With the recovery of the a-priori estimate (3.36) at time step $t^{q+2}$, we
are able to apply estimates (2.18), (2.22) (in Lemmas 2.4 and 2.5),
respectively:

(3.37)
$\displaystyle\|\tilde{\mbox{\boldmath$e$}}_{h}^{n+2}\|_{2}^{2}\geq(1-{k}^{\frac{5}{4}})\|\mbox{\boldmath$e$}_{h}^{n+2}\|_{2}^{2}+(1-{k}^{\frac{1}{4}})\|\tilde{\mbox{\boldmath$e$}}_{h}^{n+2}-\mbox{\boldmath$e$}_{h}^{n+2}\|_{2}^{2},$
(3.38)
$\displaystyle\Big{|}\langle\tilde{\mbox{\boldmath$e$}}_{h}^{n+2}-\mbox{\boldmath$e$}_{h}^{n+2},\mbox{\boldmath$e$}_{h}^{p}\rangle\Big{|}\leq{k}^{\frac{5}{4}}\|\mbox{\boldmath$e$}_{h}^{p}\|_{2}^{2}+{k}^{\frac{1}{4}}\|\tilde{\mbox{\boldmath$e$}}_{h}^{n+2}-\mbox{\boldmath$e$}_{h}^{n+2}\|_{2}^{2},\,\,p=n+1,n+2.$

Moreover, the following inequality becomes available for $n\leq q$:

(3.39)
$\displaystyle\|2\tilde{\mbox{\boldmath$e$}}_{h}^{n+2}-\mbox{\boldmath$e$}_{h}^{n+1}\|_{2}^{2}=$
$\displaystyle\|\hat{\mbox{\boldmath$e$}}_{h}^{n+3}\|_{2}^{2}+4\langle\tilde{\mbox{\boldmath$e$}}_{h}^{n+2}-\mbox{\boldmath$e$}_{h}^{n+2},\hat{\mbox{\boldmath$e$}}_{h}^{n+3}\rangle+4\|\tilde{\mbox{\boldmath$e$}}_{h}^{n+2}-\mbox{\boldmath$e$}_{h}^{n+2}\|_{2}^{2}$
$\displaystyle\geq$
$\displaystyle\|2\mbox{\boldmath$e$}_{h}^{n+2}-\mbox{\boldmath$e$}_{h}^{n+1}\|_{2}^{2}-8{k}^{\frac{5}{4}}\|\mbox{\boldmath$e$}_{h}^{n+2}\|_{2}^{2}-8{k}^{\frac{1}{4}}\|\tilde{\mbox{\boldmath$e$}}_{h}^{n+2}-\mbox{\boldmath$e$}_{h}^{n+2}\|_{2}^{2}$
$\displaystyle-4{k}^{\frac{5}{4}}\|\mbox{\boldmath$e$}_{h}^{n+1}\|_{2}^{2}-4{k}^{\frac{1}{4}}\|\tilde{\mbox{\boldmath$e$}}_{h}^{n+2}-\mbox{\boldmath$e$}_{h}^{n+2}\|_{2}^{2}+4\|\tilde{\mbox{\boldmath$e$}}_{h}^{n+2}-\mbox{\boldmath$e$}_{h}^{n+2}\|_{2}^{2}$
$\displaystyle\geq$
$\displaystyle\|2\mbox{\boldmath$e$}_{h}^{n+2}-\mbox{\boldmath$e$}_{h}^{n+1}\|_{2}^{2}-{k}(\|\mbox{\boldmath$e$}_{h}^{n+2}\|_{2}^{2}+\|\mbox{\boldmath$e$}_{h}^{n+1}\|_{2}^{2})$
$\displaystyle+(4-12{k}^{\frac{1}{4}})\|\tilde{\mbox{\boldmath$e$}}_{h}^{n+2}-\mbox{\boldmath$e$}_{h}^{n+2}\|_{2}^{2}.$

Going back (3.32), we arrive at the following estimate, for $n\leq q$:

(3.40)
$\displaystyle\frac{1}{4{k}}\big{(}\|\mbox{\boldmath$e$}_{h}^{n+2}\|_{2}^{2}-\|\mbox{\boldmath$e$}_{h}^{n+1}\|_{2}^{2}+\|2\mbox{\boldmath$e$}_{h}^{n+2}-\mbox{\boldmath$e$}_{h}^{n+1}\|_{2}^{2}-\|2\mbox{\boldmath$e$}_{h}^{n+1}-\mbox{\boldmath$e$}_{h}^{n}\|_{2}^{2})$
$\displaystyle+\frac{5-13{k}^{\frac{1}{4}}}{4{k}}\|\tilde{\mbox{\boldmath$e$}}_{h}^{n+2}-\mbox{\boldmath$e$}_{h}^{n+2}\|_{2}^{2}+\Big{(}\alpha_{2}-\frac{\gamma_{0}}{12}\Big{)}\|\nabla_{h}\tilde{\mbox{\boldmath$e$}}_{h}^{n+2}\|_{2}^{2}$
$\displaystyle-(1+\delta)(1+\frac{\gamma_{0}}{12})\Big{(}\|\nabla_{h}\tilde{\mbox{\boldmath$e$}}_{h}^{n+1}\|_{2}^{2}+\frac{1}{2}\|\nabla_{h}\tilde{\mbox{\boldmath$e$}}_{h}^{n}\|_{2}^{2}\Big{)}$
$\displaystyle\leq$
$\displaystyle\mathcal{C}_{\delta}\sum_{p=n}^{n+2}(\|\tilde{\mbox{\boldmath$e$}}_{h}^{p}\|_{2}^{2}+\|\mbox{\boldmath$e$}_{h}^{p}\|_{2}^{2}\Big{)}+\mathcal{C}({k}^{4}+h^{4}).$

Meanwhile, for the terms $\|\tilde{\mbox{\boldmath$e$}}_{h}^{p}\|_{2}^{2}$,
$p=n,n+1,n+2$, an application of Cauchy inequality gives

(3.41) $\|\tilde{\mbox{\boldmath$e$}}_{h}^{p}\|_{2}^{2}\leq
2(\|\mbox{\boldmath$e$}_{h}^{p}\|_{2}^{2}+\|\tilde{\mbox{\boldmath$e$}}_{h}^{p}-\mbox{\boldmath$e$}_{h}^{p}\|_{2}^{2}),\quad
p=n,n+1,n+2.$

Its substitution into (3.40) leads to the following inequality for $n\leq q$:

(3.42)
$\displaystyle\frac{1}{4{k}}\big{(}\|\mbox{\boldmath$e$}_{h}^{n+2}\|_{2}^{2}-\|\mbox{\boldmath$e$}_{h}^{n+1}\|_{2}^{2}+\|2\mbox{\boldmath$e$}_{h}^{n+2}-\mbox{\boldmath$e$}_{h}^{n+1}\|_{2}^{2}-\|2\mbox{\boldmath$e$}_{h}^{n+1}-\mbox{\boldmath$e$}_{h}^{n}\|_{2}^{2})$
$\displaystyle+\frac{1}{4{k}}\|\tilde{\mbox{\boldmath$e$}}_{h}^{n+2}-\mbox{\boldmath$e$}_{h}^{n+2}\|_{2}^{2}-{\mathcal{C}}_{\delta}(\|\tilde{\mbox{\boldmath$e$}}_{h}^{n+1}-\mbox{\boldmath$e$}_{h}^{n+1}\|_{2}^{2}+\|\tilde{\mbox{\boldmath$e$}}_{h}^{n}-\mbox{\boldmath$e$}_{h}^{n}\|_{2}^{2})$
$\displaystyle+\Big{(}\alpha_{2}-\frac{\gamma_{0}}{12}\Big{)}\|\nabla_{h}\tilde{\mbox{\boldmath$e$}}_{h}^{n+2}\|_{2}^{2}-(1+\delta)(1+\frac{\gamma_{0}}{12})\Big{(}\|\nabla_{h}\tilde{\mbox{\boldmath$e$}}_{h}^{n+1}\|_{2}^{2}+\frac{1}{2}\|\nabla_{h}\tilde{\mbox{\boldmath$e$}}_{h}^{n}\|_{2}^{2}\Big{)}$
$\displaystyle\leq$
$\displaystyle\mathcal{C}_{\delta}\Big{(}\|\mbox{\boldmath$e$}_{h}^{n+2}\|_{2}^{2}+\|\mbox{\boldmath$e$}_{h}^{n+1}\|_{2}^{2}+\|\mbox{\boldmath$e$}_{h}^{n}\|_{2}^{2}\Big{)}+\mathcal{C}({k}^{4}+h^{4}).$

In turn, an application of discrete Gronwall inequality [Girault1986],
combined with the fact (3.33), yields the desired convergence estimate at
$t^{q+2}$ :

(3.43)
$\displaystyle\|\mbox{\boldmath$e$}_{h}^{n}\|_{2}^{2}+\gamma_{0}{k}\sum_{p=0}^{n}\|\nabla_{h}\tilde{\mbox{\boldmath$e$}}_{h}^{p}\|_{2}^{2}\leq\mathcal{C}Te^{\mathcal{C}T}({k}^{4}+h^{4}),\quad\forall\,n\leq
q+2\leq\left\lfloor\frac{T}{k}\right\rfloor,$
$\displaystyle\mbox{i.e.,}\quad\|\mbox{\boldmath$e$}_{h}^{n}\|_{2}+\Big{(}\gamma_{0}{k}\sum_{p=0}^{n}\|\nabla_{h}\tilde{\mbox{\boldmath$e$}}_{h}^{p}\|_{2}^{2}\Big{)}^{\frac{1}{2}}\leq\mathcal{C}({k}^{2}+h^{2}).$

The convergence estimate (2.6) has been proved at $t^{q+2}$. In addition, we
see that the a-priori assumption (3.13) has also been validated at the next
time step $t^{q+2}$, provided that ${k}$ and $h$ are sufficiently small. By
mathematical induction, this completes the proof of Theorem 2.1. $\hfill\Box$

######  Remark 3.1.

The condition $\alpha>3$ is a relatively strong constraint. In fact, such a
condition is used in the estimate (3.24) for $\tilde{I}_{1}$, since we need
$\alpha>3$ to control these Laplace terms, due to the explicit treatment of
$\Delta_{h}\hat{{\mbox{\boldmath$m$}}}_{h}^{n+2}$. Meanwhile, such an
inequality only stands for a theoretical difficulty, and the practical
computations may not need that large value of $\alpha$. In most practical
simulation examples, a value of $\alpha>1$ would be sufficient to ensure the
numerical stability of the proposed numerical scheme eq. 2.3-eq. 2.4.

In addition, the explicit treatment of the Laplace term, namely
$\Delta_{h}\hat{\mbox{\boldmath$m$}}_{h}^{n+2}=\Delta_{h}(2\mbox{\boldmath$m$}_{h}^{n+1}-\mbox{\boldmath$m$}_{h}^{n})$,
will greatly improve the numerical efficiency, since only a constant-
coefficient Poisson solver is needed at each step. This crucial fact enables
one to produce very robust simulation results at a much-reduced computational
cost.

######  Remark 3.2.

In a recent work [Xie2021BDF2Linear], a rough stability estimate for the
projection step, namely, $\|\mbox{\boldmath$e$}_{h}\|_{2}\leq
2\|\tilde{\mbox{\boldmath$e$}}_{h}\|_{2}+{\mathcal{O}}(h^{2})$,
$\|\nabla_{h}\mbox{\boldmath$e$}_{h}\|_{2}\leq{\mathcal{C}}(\|\mbox{\boldmath$e$}_{h}\|_{2}+\|\nabla_{h}\tilde{\mbox{\boldmath$e$}}_{h}\|_{2})+{\mathcal{O}}(h^{2})$,
has been proved. These inequalities are sufficient to establish the stability
and convergence analysis for a semi-implicit numerical scheme, if the BDF2
temporal stencil is formulated as
$\frac{1}{k}(\frac{3}{2}\tilde{\mbox{\boldmath$m$}}_{h}^{n+2}-2\tilde{\mbox{\boldmath$m$}}_{h}^{n+1}+\frac{1}{2}\tilde{\mbox{\boldmath$m$}}_{h}^{n})$.
However, for the BDF2 temporal stencil formulated as
$\frac{1}{k}(\frac{3}{2}\tilde{\mbox{\boldmath$m$}}_{h}^{n+2}-2\mbox{\boldmath$m$}_{h}^{n+1}+\frac{1}{2}\mbox{\boldmath$m$}_{h}^{n})$
(as given by Algorithm 2.2), such a stability estimate are not sufficient to
derive the stability and convergence analysis, due to the singular coefficient
$\frac{1}{k}$. Instead, a much more refined stability estimate, as given by
(2.18), (2.19) (in Lemma 2.4), is needed to pass through the convergence
analysis. The proof of these two refined inequalities has to be based on a
more precise geometric analysis of the corresponding vectors, and the details
will be presented in Appendix A.

Extensive numerical experiments have demonstrated a better stability property
of the temporal stencil in Algorithm 2.2 than that of Algorithm 2.1, for
physical models with large damping parameters. For the theoretical analysis of
Algorithm 2.2, the refined stability estimate (2.18), (2.19) has played a
crucial role.

##  4\. Conclusions

In this paper, we have presented an optimal rate convergence analysis and
error estimate for a second-order accurate, linear numerical scheme to the LL
equation. The second-order backward differentiation formula is applied in the
temporal discretization, the linear diffusion term is treated implicitly,
while the nonlinear terms are updated by a fully explicit extrapolation
formula. Afterward, a point-wise projection is applied to normalize the
magnetization vector. In turn, only a linear system independent of the updated
magnetization needs to be solved at each time step, which has greatly improved
the computational efficiency, and many great advantages of this numerical
scheme have been reported in the numerical simulation with large damping
parameters. The error estimate has been theoretically established in the
discrete $\ell^{\infty}(0,T;\ell^{2})\cap\ell^{2}(0,T;H_{h}^{1})$ norm, under
suitable regularity assumptions and reasonable ratio between the time step-
size and the spatial mesh-size. The key difficulty of the theoretical analysis
is associated with the fact that the projection step is highly nonlinear and
non-convex. To overcome this subtle difficulty, we build a stability estimate
for the projection step, which plays a crucial role in the derivation of the
convergence analysis for the numerical scheme.

## Appendix A. Proof of 2.4

###### Proof.

First of all, an $\|\cdot\|_{\infty}$ bound for the numerical error
$\tilde{\mbox{\boldmath$e$}}_{h}$ can be derived, by the a-priori estimate
(2.17):

(4.1) $\|\tilde{\mbox{\boldmath$e$}}_{h}\|_{\infty}\leq\gamma
h^{-\frac{1}{2}}(\|\tilde{\mbox{\boldmath$e$}}_{h}\|_{2}+\|\nabla_{h}\tilde{\mbox{\boldmath$e$}}_{h}\|_{2})\leq\gamma
h^{-\frac{1}{2}}\cdot{k}^{\frac{11}{8}}\leq\frac{1}{4}{k}^{\frac{3}{4}},$

provided that ${k}$ and $h$ are sufficiently small, and under the linear
refinement requirement $C_{1}h\leq{k}\leq C_{2}h$. Notice that the inverse
inequality (2.7) has been applied in the first step. In turn, we observe that
the following bounds are available for the numerical profile
$\tilde{\mbox{\boldmath$m$}}_{h}$:

(4.2) $\displaystyle
1-\frac{1}{4}{k}^{\frac{3}{4}}\leq|\tilde{\mbox{\boldmath$m$}}_{h}|\leq
1+\frac{1}{4}{k}^{\frac{3}{4}},\quad\mbox{at a point-wise level},$ (4.3)
$\displaystyle\|\nabla_{h}\tilde{\mbox{\boldmath$e$}}_{h}\|_{4}\leq\gamma
h^{-\frac{3}{4}}\|\nabla_{h}\tilde{\mbox{\boldmath$e$}}_{h}\|_{2}\leq\gamma
h^{-\frac{3}{4}}\cdot\frac{1}{2}{k}^{\frac{11}{8}}\leq{k}^{\frac{1}{2}},$
(4.4)
$\displaystyle\|\nabla_{h}\tilde{\mbox{\boldmath$m$}}_{h}\|_{4}\leq\|\nabla_{h}\underline{\bm{m}}_{h}\|_{4}+\|\nabla_{h}\tilde{\mbox{\boldmath$e$}}_{h}\|_{4}\leq
C^{*}+1:=M.$

Again, the inverse inequality (2.7) has been applied in the derivation.

A careful calculation indicates that

(4.5) $\displaystyle{\mbox{\boldmath$e$}}_{h}=$
$\displaystyle\underline{\bm{m}}_{h}-{\mbox{\boldmath$m$}}_{h}=\underline{\bm{m}}_{h}-\tilde{\mbox{\boldmath$m$}}_{h}+\tilde{\mbox{\boldmath$m$}}_{h}-\frac{\tilde{\mbox{\boldmath$m$}}_{h}}{|\tilde{\mbox{\boldmath$m$}}_{h}|}=\tilde{\mbox{\boldmath$e$}}_{h}+\frac{\tilde{\mbox{\boldmath$m$}}_{h}}{|\tilde{\mbox{\boldmath$m$}}_{h}|}(|\tilde{\mbox{\boldmath$m$}}_{h}|-1),\quad\mbox{so
that}$ $\displaystyle\tilde{\mbox{\boldmath$e$}}_{h}=$
$\displaystyle\mbox{\boldmath$e$}_{h}+\tilde{\mbox{\boldmath$e$}}_{h,c},\quad\tilde{\mbox{\boldmath$e$}}_{h,c}:=\frac{\tilde{\mbox{\boldmath$m$}}_{h}}{|\tilde{\mbox{\boldmath$m$}}_{h}|}(1-|\tilde{\mbox{\boldmath$m$}}_{h}|).$

Meanwhile, at a fixed grid point
$(\hat{x}_{i},\hat{y}_{j},\hat{z}_{\ell})\in\Omega_{h}^{0}$, we look at the
triangle formed by the vectors: $\tilde{\mbox{\boldmath$m$}}_{h}$,
$\tilde{\mbox{\boldmath$e$}}_{h}$ and $\underline{\bm{m}}_{h}$. In particular,
we denote the angle between $\tilde{\mbox{\boldmath$m$}}_{h}$ and
$\underline{\bm{m}}_{h}$ as $\theta$. Since the lengths of these three vectors
have the following estimates:

$1-\frac{1}{4}{k}^{\frac{3}{4}}\leq|\tilde{\mbox{\boldmath$m$}}_{h}|\leq
1+\frac{1}{4}{k}^{\frac{3}{4}},\,\,\,|\tilde{\mbox{\boldmath$e$}}_{h}|\leq\frac{1}{4}{k}^{\frac{3}{4}},\,\,\,|\underline{\bm{m}}_{h}|\equiv
1,$

a careful application of Sine law indicates that

(4.6) $0\leq\sin\theta\leq\frac{1}{4}{k}^{\frac{3}{4}}.$

And also, we look at the triangle formed by the vectors:
$\mbox{\boldmath$m$}_{h}$, $\mbox{\boldmath$e$}_{h}$ and
$\underline{\bm{m}}_{h}$. Two sides have equal lengths:
$|\mbox{\boldmath$m$}_{h}|=|\underline{\bm{m}}_{h}|=1$, and the angle between
$\mbox{\boldmath$m$}_{h}$ and $\underline{\bm{m}}_{h}$ is exactly $\theta$,
since
$\mbox{\boldmath$m$}_{h}=\frac{\tilde{\mbox{\boldmath$m$}}_{h}}{|\tilde{\mbox{\boldmath$m$}}_{h}|}$
is in the same direction as $\tilde{\mbox{\boldmath$m$}}_{h}$. In turn, the
angle between vectors $\mbox{\boldmath$m$}_{h}$ and $\mbox{\boldmath$e$}_{h}$
is given by $\varphi^{\prime}=\frac{\pi}{2}-\frac{\theta}{2}$. Because of this
fact, the angle between vectors $\tilde{\mbox{\boldmath$m$}}_{h}$ and
$\mbox{\boldmath$e$}_{h}$ is the same as
$\varphi^{\prime}=\frac{\pi}{2}-\frac{\theta}{2}$,

Subsequently, we denote the angle between vectors $\tilde{e}_{h,c}$ and
$\mbox{\boldmath$e$}_{h}$ as $\varphi$. In fact, this angle has the following
representation

(4.7)
$\displaystyle\varphi=\varphi^{\prime}=\frac{\pi}{2}-\frac{\theta}{2},\quad\mbox{if}\,\,\,|\tilde{\mbox{\boldmath$m$}}_{h}|\leq
1;\quad\varphi=\pi-\varphi^{\prime}=\frac{\pi}{2}+\frac{\theta}{2},\quad\mbox{if}\,\,\,|\tilde{\mbox{\boldmath$m$}}_{h}|>1.$

In either case, the following estimate is valid:

(4.8)
$|\cos\varphi|=|\sin\frac{\theta}{2}|\leq\sin\theta\leq\frac{1}{4}{k}^{\frac{3}{4}}.$

Consequently, the definition of the point-wise inner product implies the
following estimate

(4.9) $\displaystyle|\tilde{\mbox{\boldmath$e$}}_{h}|^{2}=$
$\displaystyle|\mbox{\boldmath$e$}_{h}+\tilde{\mbox{\boldmath$e$}}_{h,c}|^{2}=|\mbox{\boldmath$e$}_{h}|^{2}+|\tilde{\mbox{\boldmath$e$}}_{h,c}|^{2}+2\mbox{\boldmath$e$}_{h}\cdot\tilde{\mbox{\boldmath$e$}}_{h,c}$
$\displaystyle=$
$\displaystyle|\mbox{\boldmath$e$}_{h}|^{2}+|\tilde{\mbox{\boldmath$e$}}_{h,c}|^{2}+2|\mbox{\boldmath$e$}_{h}|\cdot|\tilde{\mbox{\boldmath$e$}}_{h,c}|\cdot\cos\varphi$
$\displaystyle\geq$
$\displaystyle|\mbox{\boldmath$e$}_{h}|^{2}+|\tilde{\mbox{\boldmath$e$}}_{h,c}|^{2}-2|\mbox{\boldmath$e$}_{h}|\cdot|\tilde{\mbox{\boldmath$e$}}_{h,c}|\cdot\frac{1}{4}{k}^{\frac{3}{4}}=|\mbox{\boldmath$e$}_{h}|^{2}+|\tilde{\mbox{\boldmath$e$}}_{h,c}|^{2}-\frac{1}{2}{k}^{\frac{3}{4}}|\mbox{\boldmath$e$}_{h}|\cdot|\tilde{\mbox{\boldmath$e$}}_{h,c}|$
$\displaystyle\geq$
$\displaystyle|\mbox{\boldmath$e$}_{h}|^{2}+|\tilde{\mbox{\boldmath$e$}}_{h,c}|^{2}-({k}^{\frac{5}{4}}|\mbox{\boldmath$e$}_{h}|^{2}+\frac{1}{4}{k}^{\frac{1}{4}}|\tilde{\mbox{\boldmath$e$}}_{h,c}|^{2})$
$\displaystyle\geq$
$\displaystyle(1-{k}^{\frac{5}{4}})|\mbox{\boldmath$e$}_{h}|^{2}+(1-{k}^{\frac{1}{4}})|\tilde{\mbox{\boldmath$e$}}_{h}-\mbox{\boldmath$e$}_{h}|^{2}.$

Notice that this estimate is valid at a point-wise level, for any fixed mesh
point $(i,j,\ell)$. Therefore, a summation in space leads to the first
inequality (2.18).

To derive the second inequality (2.19), we will focus on the $D_{x}$ part in
the discrete gradient; the analysis for the $D_{y}$ and $D_{z}$ parts can be
performed in a similar manner. We begin with the following expansion:

(4.10) $\displaystyle D_{x}\mbox{\boldmath$e$}_{h}$
$\displaystyle=D_{x}{\underline{\bm{m}}}_{h}-D_{x}\frac{\tilde{\mbox{\boldmath$m$}}_{h}}{|\tilde{\mbox{\boldmath$m$}}_{h}|}=D_{x}\left[\frac{{\underline{\bm{m}}}_{h}}{|\tilde{\mbox{\boldmath$m$}}_{h}|}-\frac{\tilde{\mbox{\boldmath$m$}}_{h}}{|\tilde{\mbox{\boldmath$m$}}_{h}|}\right]+D_{x}\left[\underline{\bm{m}}_{h}-\frac{{\underline{\bm{m}}}_{h}}{|\tilde{\mbox{\boldmath$m$}}_{h}|}\right]$
$\displaystyle=D_{x}\frac{\tilde{\mbox{\boldmath$e$}}_{h}}{|\tilde{\mbox{\boldmath$m$}}_{h}|}-D_{x}\left[\frac{{\underline{\bm{m}}}_{h}}{|\tilde{\mbox{\boldmath$m$}}_{h}|}(1-|\tilde{\mbox{\boldmath$m$}}_{h}|)\right]$
$\displaystyle=D_{x}\frac{\tilde{\mbox{\boldmath$e$}}_{h}}{|\tilde{\mbox{\boldmath$m$}}_{h}|}-D_{x}\left[\frac{{\underline{\bm{m}}}_{h}}{|\tilde{\mbox{\boldmath$m$}}_{h}|}\frac{(\underline{\bm{m}}_{h}+\tilde{\mbox{\boldmath$m$}}_{h})\cdot(\underline{\bm{m}}_{h}-\tilde{\mbox{\boldmath$m$}}_{h})}{1+|\tilde{\mbox{\boldmath$m$}}_{h}|}\right]$
$\displaystyle=D_{x}\frac{\tilde{\mbox{\boldmath$e$}}_{h}}{|\tilde{\mbox{\boldmath$m$}}_{h}|}-D_{x}\left[\underline{\bm{m}}_{h}\frac{(\underline{\bm{m}}_{h}+\tilde{\mbox{\boldmath$m$}}_{h})\cdot\tilde{\mbox{\boldmath$e$}}_{h}}{|\tilde{\mbox{\boldmath$m$}}_{h}|+|\tilde{\mbox{\boldmath$m$}}_{h}|^{2}}\right]$
$\displaystyle=D_{x}\frac{\tilde{\mbox{\boldmath$e$}}_{h}}{|\tilde{\mbox{\boldmath$m$}}_{h}|}-D_{x}\left[\underline{\bm{m}}_{h}\frac{(2\underline{\bm{m}}_{h}-\tilde{\mbox{\boldmath$e$}}_{h})\cdot\tilde{\mbox{\boldmath$e$}}_{h}}{|\tilde{\mbox{\boldmath$m$}}_{h}|+|\tilde{\mbox{\boldmath$m$}}_{h}|^{2}}\right]$
$\displaystyle=D_{x}\frac{\tilde{\mbox{\boldmath$e$}}_{h}}{|\tilde{\mbox{\boldmath$m$}}_{h}|}-D_{x}\Big{(}\underline{\bm{m}}_{h}\frac{2\underline{\bm{m}}_{h}\cdot\tilde{\mbox{\boldmath$e$}}_{h}}{|\tilde{\mbox{\boldmath$m$}}_{h}|+|\tilde{\mbox{\boldmath$m$}}_{h}|^{2}}\Big{)}+D_{x}\Big{(}\underline{\bm{m}}_{h}\frac{|\tilde{\mbox{\boldmath$e$}}_{h}|^{2}}{|\tilde{\mbox{\boldmath$m$}}_{h}|+|\tilde{\mbox{\boldmath$m$}}_{h}|^{2}}\Big{)},$

in which the identity
$\tilde{\mbox{\boldmath$m$}}_{h}=\underline{\bm{m}}_{h}-\tilde{\mbox{\boldmath$e$}}_{h}$
has been applied at the last two steps. Meanwhile, at each numerical mesh
cell, from $(i,j,\ell)$ to $(i+1,j,\ell)$, the identity (2.10) is valid, so
that the following expansions are able to be made at the center location
$(i+1/2,j,\ell)$ over the numerical cell:

(4.11) $\displaystyle
D_{x}\frac{\tilde{\mbox{\boldmath$e$}}_{h}}{|\tilde{\mbox{\boldmath$m$}}_{h}|}-D_{x}\Big{(}\underline{\bm{m}}_{h}\frac{2\underline{\bm{m}}_{h}\cdot\tilde{\mbox{\boldmath$e$}}_{h}}{|\tilde{\mbox{\boldmath$m$}}_{h}|+|\tilde{\mbox{\boldmath$m$}}_{h}|^{2}}\Big{)}=J_{1}+J_{2}+J_{3}+J_{4}+J_{5},$
$\displaystyle J_{1}=$ $\displaystyle
A_{x}\Big{(}\frac{1}{|\tilde{\mbox{\boldmath$m$}}_{h}|}\Big{)}D_{x}\tilde{\mbox{\boldmath$e$}}_{h}-A_{x}\underline{\bm{m}}_{h}A_{x}\Big{(}\frac{2}{|\tilde{\mbox{\boldmath$m$}}_{h}|+|\tilde{\mbox{\boldmath$m$}}_{h}|^{2}}\Big{)}(A_{x}\underline{\bm{m}}_{h}\cdot
D_{x}\tilde{\mbox{\boldmath$e$}}_{h}),$ $\displaystyle J_{2}=$
$\displaystyle-\Big{(}A_{x}^{(2)}\Big{(}\frac{1}{|\tilde{\mbox{\boldmath$m$}}_{h}|^{2}}\Big{)}D_{x}|\tilde{\mbox{\boldmath$m$}}_{h}|\Big{)}\tilde{\mbox{\boldmath$e$}}_{h},\quad
J_{3}=-2(D_{x}\underline{\bm{m}}_{h})A_{x}\Big{(}\frac{\underline{\bm{m}}_{h}\cdot\tilde{\mbox{\boldmath$e$}}_{h}}{|\tilde{\mbox{\boldmath$m$}}_{h}|+|\tilde{\mbox{\boldmath$m$}}_{h}|^{2}}\Big{)},$
$\displaystyle J_{4}=$
$\displaystyle(A_{x}\underline{\bm{m}}_{h})A_{x}(\underline{\bm{m}}_{h}\cdot\tilde{\mbox{\boldmath$e$}}_{h})A_{x}^{(2)}\Big{(}\frac{2}{(|\tilde{\mbox{\boldmath$m$}}_{h}|+|\tilde{\mbox{\boldmath$m$}}_{h}|^{2})^{2}}\Big{)}$
$\displaystyle\cdot(D_{x}|\tilde{\mbox{\boldmath$m$}}_{h}|+2A_{x}\tilde{\mbox{\boldmath$m$}}_{h}\cdot
D_{x}\tilde{\mbox{\boldmath$m$}}_{h}),$ $\displaystyle J_{5}=$ $\displaystyle-
A_{x}\underline{\bm{m}}_{h}A_{x}\Big{(}\frac{2}{|\tilde{\mbox{\boldmath$m$}}_{h}|+|\tilde{\mbox{\boldmath$m$}}_{h}|^{2}}\Big{)}(D_{x}\underline{\bm{m}}_{h}\cdot
A_{x}\tilde{\mbox{\boldmath$e$}}_{h}),$

in which the nonlinear average operator $A_{x}^{(2)}$ is introduced as

(4.12)
$A_{x}^{(2)}\Big{(}\frac{1}{(f_{h})^{2}}\Big{)}_{i,j,\ell}=\frac{1}{(f_{h})_{i,j,\ell}(f_{h})_{i+1,j,\ell}},\,\,\mbox{for
scalar grid function $f_{h}$}.$

Moreover, by the point-wise a-priori estimate (4.2) for
$\tilde{\mbox{\boldmath$m$}}_{h}$, the following bounds are available:

(4.13) $\displaystyle\frac{1}{1+\frac{1}{4}{k}^{\frac{3}{4}}}\leq
A_{x}\Big{(}\frac{1}{|\tilde{\mbox{\boldmath$m$}}_{h}|}\Big{)}\leq\frac{1}{1-\frac{1}{4}{k}^{\frac{3}{4}}},$
$\displaystyle 1-\frac{3}{2}{k}^{\frac{3}{4}}\leq
A_{x}\Big{(}\frac{2}{|\tilde{\mbox{\boldmath$m$}}_{h}|+|\tilde{\mbox{\boldmath$m$}}_{h}|^{2}}\Big{)}\leq
1+\frac{3}{2}{k}^{\frac{3}{4}},$
$\displaystyle\frac{1}{(1+\frac{1}{4}{k}^{\frac{3}{4}})^{2}}\leq
A_{x}^{(2)}\Big{(}\frac{1}{|\tilde{\mbox{\boldmath$m$}}_{h}|^{2}}\Big{)}\leq\frac{1}{(1-\frac{1}{4}{k}^{\frac{3}{4}})^{2}},$
$\displaystyle\frac{1}{2(1+\frac{1}{4}{k}^{\frac{3}{4}})^{4}}\leq
A_{x}^{(2)}\Big{(}\frac{2}{(|\tilde{\mbox{\boldmath$m$}}_{h}|+|\tilde{\mbox{\boldmath$m$}}_{h}|^{2})^{2}}\Big{)}\leq\frac{1}{(1-\frac{1}{4}{k}^{\frac{3}{4}})^{4}}.$

For the term $J_{1}$, we make the following decomposition
$J_{1}=J_{11}+J_{12}+J_{13}$ to facilitate the analysis

(4.14) $\displaystyle J_{11}=$ $\displaystyle
D_{x}\tilde{\mbox{\boldmath$e$}}_{h}-A_{x}\underline{\bm{m}}_{h}(A_{x}\underline{\bm{m}}_{h}\cdot
D_{x}\tilde{\mbox{\boldmath$e$}}_{h}),\quad
J_{12}=\Big{(}A_{x}\Big{(}\frac{1}{|\tilde{\mbox{\boldmath$m$}}_{h}|}\Big{)}-1\Big{)}D_{x}\tilde{\mbox{\boldmath$e$}}_{h},$
$\displaystyle J_{13}=$ $\displaystyle-
A_{x}\underline{\bm{m}}_{h}\Big{(}A_{x}\Big{(}\frac{2}{|\tilde{\mbox{\boldmath$m$}}_{h}|+|\tilde{\mbox{\boldmath$m$}}_{h}|^{2}}\Big{)}-1\Big{)}(A_{x}\underline{\bm{m}}_{h}\cdot
D_{x}\tilde{\mbox{\boldmath$e$}}_{h}).$

For the quantity $J_{11}$, we see that

(4.15) $\displaystyle|J_{11}|^{2}=$
$\displaystyle|D_{x}\tilde{\mbox{\boldmath$e$}}_{h}|^{2}+|A_{x}\underline{\bm{m}}_{h}|^{2}(A_{x}\underline{\bm{m}}_{h}\cdot
D_{x}\tilde{\mbox{\boldmath$e$}}_{h})^{2}-2(A_{x}\underline{\bm{m}}_{h}\cdot
D_{x}\tilde{\mbox{\boldmath$e$}}_{h})^{2}$ $\displaystyle\leq$
$\displaystyle|D_{x}\tilde{\mbox{\boldmath$e$}}_{h}|^{2}-(A_{x}\underline{\bm{m}}_{h}\cdot
D_{x}\tilde{\mbox{\boldmath$e$}}_{h})^{2},\quad\mbox{so
that}\,\,\,|J_{11}|\leq|D_{x}\tilde{\mbox{\boldmath$e$}}_{h}|,$

in which the last step comes from the fact that
$|A_{x}\underline{\bm{m}}_{h}|^{2}\leq 1$. For the quantity $J_{12}$, the
a-priori bounds (4.13) imply that

(4.16)
$J_{12}=\Big{|}A_{x}\Big{(}\frac{1}{|\tilde{\mbox{\boldmath$m$}}_{h}|}\Big{)}-1\Big{|}\cdot|D_{x}\tilde{\mbox{\boldmath$e$}}_{h}|\leq\frac{1}{2}{k}^{\frac{3}{4}}|D_{x}\tilde{\mbox{\boldmath$e$}}_{h}|.$

Similarly, the quantity $J_{13}$ can be analyzed as

(4.17) $\displaystyle J_{13}=$
$\displaystyle|A_{x}\underline{\bm{m}}_{h}|\cdot\Big{|}A_{x}\Big{(}\frac{2}{|\tilde{\mbox{\boldmath$m$}}_{h}|+|\tilde{\mbox{\boldmath$m$}}_{h}|^{2}}\Big{)}-1\Big{|}\cdot|A_{x}\underline{\bm{m}}_{h}|\cdot|D_{x}\tilde{\mbox{\boldmath$e$}}_{h}|$
$\displaystyle\leq$
$\displaystyle(1+Ch^{2})^{2}(\frac{1}{(1-\frac{1}{4}{k}^{\frac{3}{4}})^{2}}|D_{x}\tilde{\mbox{\boldmath$e$}}_{h}|\leq(1+\frac{3}{4}{k}^{\frac{3}{4}})|D_{x}\tilde{\mbox{\boldmath$e$}}_{h}|.$

Consequently, a combination of (4.15)-(4.17) leads to

(4.18)
$|J_{1}|\leq|J_{11}|+|J_{12}|+|J_{13}|\leq(1+\frac{5}{4}{k}^{\frac{3}{4}})|D_{x}\tilde{\mbox{\boldmath$e$}}_{h}|,$

which turns out to be a point-wise inequality. In turn, a summation in space
implies that

(4.19)
$\|J_{1}\|_{2}\leq(1+\frac{5}{4}{k}^{\frac{3}{4}})\|D_{x}\tilde{\mbox{\boldmath$e$}}_{h}\|_{2}.$

For the term $J_{2}$, we make use of the a-priori bounds (4.13), as well as
the fact that
$|D_{x}|\tilde{\mbox{\boldmath$m$}}_{h}||\leq|D_{x}\tilde{\mbox{\boldmath$m$}}_{h}|$,
so that an application of discrete Hölder inequality gives

(4.20) $\displaystyle\|J_{2}\|_{2}\leq$
$\displaystyle\Big{\|}A_{x}^{(2)}\Big{(}\frac{1}{|\tilde{\mbox{\boldmath$m$}}_{h}|^{2}}\Big{)}\Big{\|}_{\infty}\cdot\|D_{x}\tilde{\mbox{\boldmath$m$}}_{h}\|_{4}\cdot\|\tilde{\mbox{\boldmath$e$}}_{h}\|_{4}\leq(1+\frac{3}{4}{k}^{\frac{3}{4}})M\|\tilde{\mbox{\boldmath$e$}}_{h}\|_{4},$

in which the a-priori $W_{h}^{1,4}$ bound (4.4) for the numerical profile
$\tilde{\mbox{\boldmath$m$}}_{h}$ has been applied. The other terms in the
expansion (4.11) can be analyzed in a similar manner.

(4.21) $\displaystyle\|J_{3}\|_{2}\leq$
$\displaystyle\|D_{x}\underline{\bm{m}}_{h}\|_{\infty}\cdot\max\Big{(}\frac{2}{|\tilde{\mbox{\boldmath$m$}}_{h}|+|\tilde{\mbox{\boldmath$m$}}_{h}|^{2}}\Big{)}\cdot\|\underline{\bm{m}}_{h}\|_{\infty}\cdot\|\tilde{\mbox{\boldmath$e$}}_{h}\|_{2}$
$\displaystyle\leq$ $\displaystyle C^{*}(1+\frac{3}{4}{k}^{\frac{3}{4}})\cdot
1\cdot\|\tilde{\mbox{\boldmath$e$}}_{h}\|_{2}\leq(1+\frac{3}{4}{k}^{\frac{3}{4}})M\|\tilde{\mbox{\boldmath$e$}}_{h}\|_{2},$
$\displaystyle\|J_{4}\|_{2}\leq$
$\displaystyle\|\underline{\bm{m}}_{h}\|_{\infty}^{2}\cdot\|\tilde{\mbox{\boldmath$e$}}_{h}\|_{4}\cdot\max\Big{(}\frac{2}{(|\tilde{\mbox{\boldmath$m$}}_{h}|+|\tilde{\mbox{\boldmath$m$}}_{h}|^{2})^{2}}\Big{)}$
$\displaystyle\cdot(1+2\|\tilde{\mbox{\boldmath$m$}}_{h}\|_{\infty})\|D_{x}\tilde{\mbox{\boldmath$m$}}_{h}\|_{4}\leq\frac{3}{2}(1+\frac{3}{2}{k}^{\frac{3}{4}})M\|\tilde{\mbox{\boldmath$e$}}_{h}\|_{4},$
$\displaystyle\|J_{5}\|_{2}\leq$
$\displaystyle\|\underline{\bm{m}}_{h}\|_{\infty}\cdot\max\Big{(}\frac{2}{|\tilde{\mbox{\boldmath$m$}}_{h}|+|\tilde{\mbox{\boldmath$m$}}_{h}|^{2}}\Big{)}\cdot\|D_{x}\underline{\bm{m}}_{h}\|_{\infty}\cdot\|\tilde{\mbox{\boldmath$e$}}_{h}\|_{2}$
$\displaystyle\leq$
$\displaystyle(1+\frac{3}{2}{k}^{\frac{3}{4}})C^{*}\|\tilde{\mbox{\boldmath$e$}}_{h}\|_{2}\leq(1+\frac{3}{2}{k}^{\frac{3}{4}})M\|\tilde{\mbox{\boldmath$e$}}_{h}\|_{2}.$

Therefore, a substitution of (4.19), (4.20) and (4.21) into (4.11) leads to

(4.22)
$\displaystyle\Big{\|}D_{x}\frac{\tilde{\mbox{\boldmath$e$}}_{h}}{|\tilde{\mbox{\boldmath$m$}}_{h}|}-D_{x}\Big{(}\underline{\bm{m}}_{h}\frac{2\underline{\bm{m}}_{h}\cdot\tilde{\mbox{\boldmath$e$}}_{h}}{|\tilde{\mbox{\boldmath$m$}}_{h}|+|\tilde{\mbox{\boldmath$m$}}_{h}|^{2}}\Big{)}\Big{\|}_{2}$
$\displaystyle\leq$
$\displaystyle(1+\frac{5}{4}{k}^{\frac{3}{4}})\|D_{x}\tilde{\mbox{\boldmath$e$}}_{h}\|_{2}+3M(\|\tilde{\mbox{\boldmath$e$}}_{h}\|_{2}+\|\tilde{\mbox{\boldmath$e$}}_{h}\|_{4}).$

The analysis for the last term on the right hand side of (4.10) can be
similarly carried out; some technical details are skipped for the sake of
brevity.

(4.23)
$\displaystyle\Big{\|}D_{x}\Big{(}\underline{\bm{m}}_{h}\frac{|\tilde{\mbox{\boldmath$e$}}_{h}|^{2}}{|\tilde{\mbox{\boldmath$m$}}_{h}|+|\tilde{\mbox{\boldmath$m$}}_{h}|^{2}}\Big{)}\Big{\|}_{2}$
$\displaystyle\leq$
$\displaystyle{\mathcal{C}}\Big{(}\|\underline{\bm{m}}_{h}\|_{\infty}\cdot\max\Big{(}\frac{2}{|\tilde{\mbox{\boldmath$m$}}_{h}|+|\tilde{\mbox{\boldmath$m$}}_{h}|^{2}}\Big{)}\cdot\|\tilde{\mbox{\boldmath$e$}}_{h}\|_{\infty}\cdot\|D_{x}\tilde{\mbox{\boldmath$e$}}_{h}\|_{2}$
$\displaystyle+\|D_{x}\underline{\bm{m}}_{h}\|_{\infty}\cdot\max\Big{(}\frac{1}{|\tilde{\mbox{\boldmath$m$}}_{h}|+|\tilde{\mbox{\boldmath$m$}}_{h}|^{2}}\Big{)}\cdot\|\tilde{\mbox{\boldmath$e$}}_{h}\|_{\infty}\cdot\|\tilde{\mbox{\boldmath$e$}}_{h}\|_{2}$
$\displaystyle+\|\underline{\bm{m}}_{h}\|_{\infty}\cdot\max\Big{(}\frac{1}{(|\tilde{\mbox{\boldmath$m$}}_{h}|+|\tilde{\mbox{\boldmath$m$}}_{h}|^{2})^{2}}\Big{)}\cdot(1+2\|\tilde{\mbox{\boldmath$m$}}_{h}\|_{\infty})\|D_{x}\tilde{\mbox{\boldmath$m$}}_{h}\|_{4}$
$\displaystyle\cdot\|\tilde{\mbox{\boldmath$e$}}_{h}\|_{\infty}\cdot\|\tilde{\mbox{\boldmath$e$}}_{h}\|_{4}\Big{)}$
$\displaystyle\leq$
$\displaystyle{\mathcal{C}}{k}^{\frac{3}{4}}(\|D_{x}\tilde{\mbox{\boldmath$e$}}_{h}\|_{2}+\|\tilde{\mbox{\boldmath$e$}}_{h}\|_{2}+\|\tilde{\mbox{\boldmath$e$}}_{h}\|_{4})\leq{k}^{\frac{5}{8}}(\|D_{x}\tilde{\mbox{\boldmath$e$}}_{h}\|_{2}+\|\tilde{\mbox{\boldmath$e$}}_{h}\|_{2}+\|\tilde{\mbox{\boldmath$e$}}_{h}\|_{4}).$

Finally, a substitution of (4.22) and (4.23) into (4.10) yields

(4.24) $\displaystyle\|D_{x}\mbox{\boldmath$e$}_{h}\|_{2}\leq$
$\displaystyle(1+2{k}^{\frac{5}{8}})\|D_{x}\tilde{\mbox{\boldmath$e$}}_{h}\|_{2}+(3M+1)(\|\tilde{\mbox{\boldmath$e$}}_{h}\|_{2}+\|\tilde{\mbox{\boldmath$e$}}_{h}\|_{4}),\quad\mbox{so
that}$ $\displaystyle\|D_{x}\mbox{\boldmath$e$}_{h}\|_{2}^{2}\leq$
$\displaystyle(1+5{k}^{\frac{5}{8}})\|D_{x}\tilde{\mbox{\boldmath$e$}}_{h}\|_{2}^{2}+(3M+1)^{2}(\|\tilde{\mbox{\boldmath$e$}}_{h}\|_{2}+\|\tilde{\mbox{\boldmath$e$}}_{h}\|_{4})^{2}$
$\displaystyle+4(3M+1)\|D_{x}\tilde{\mbox{\boldmath$e$}}_{h}\|_{2}(\|\tilde{\mbox{\boldmath$e$}}_{h}\|_{2}+\|\tilde{\mbox{\boldmath$e$}}_{h}\|_{4})$
$\displaystyle\leq$
$\displaystyle(1+5{k}^{\frac{5}{8}})\|D_{x}\tilde{\mbox{\boldmath$e$}}_{h}\|_{2}^{2}+2(3M+1)^{2}(\|\tilde{\mbox{\boldmath$e$}}_{h}\|_{2}^{2}+\|\tilde{\mbox{\boldmath$e$}}_{h}\|_{4}^{2})$
$\displaystyle+\frac{\delta}{4}\|D_{x}\tilde{\mbox{\boldmath$e$}}_{h}\|_{2}^{2}+16(3M+1)^{2}\delta^{-1}(\|\tilde{\mbox{\boldmath$e$}}_{h}\|_{2}+\|\tilde{\mbox{\boldmath$e$}}_{h}\|_{4})^{2}$
$\displaystyle\leq$
$\displaystyle(1+\frac{\delta}{2})\|D_{x}\tilde{\mbox{\boldmath$e$}}_{h}\|_{2}^{2}+(32\delta^{-1}+2)(3M+1)^{2}(\|\tilde{\mbox{\boldmath$e$}}_{h}\|_{2}^{2}+\|\tilde{\mbox{\boldmath$e$}}_{h}\|_{4}^{2}),$

for any $\delta>0$, provided that ${k}$ is sufficiently small. Similar
estimates can be derived for the gradient in the $y$ and $z$ directions; the
technical details are skipped for the sake of brevity.

(4.25) $\displaystyle\|D_{y}\mbox{\boldmath$e$}_{h}\|_{2}^{2}\leq$
$\displaystyle(1+\frac{\delta}{2})\|D_{y}\tilde{\mbox{\boldmath$e$}}_{h}\|_{2}^{2}+(32\delta^{-1}+2)(3M+1)^{2}(\|\tilde{\mbox{\boldmath$e$}}_{h}\|_{2}^{2}+\|\tilde{\mbox{\boldmath$e$}}_{h}\|_{4}^{2}),$
$\displaystyle\|D_{z}\mbox{\boldmath$e$}_{h}\|_{2}^{2}\leq$
$\displaystyle(1+\frac{\delta}{2})\|D_{z}\tilde{\mbox{\boldmath$e$}}_{h}\|_{2}^{2}+(32\delta^{-1}+2)(3M+1)^{2}(\|\tilde{\mbox{\boldmath$e$}}_{h}\|_{2}^{2}+\|\tilde{\mbox{\boldmath$e$}}_{h}\|_{4}^{2}).$

Then we arrive at

(4.26)
$\|\nabla_{h}\mbox{\boldmath$e$}_{h}\|_{2}^{2}\leq(1+\frac{\delta}{2})\|\nabla_{h}\tilde{\mbox{\boldmath$e$}}_{h}\|_{2}^{2}+3(32\delta^{-1}+2)(3M+1)^{2}(\|\tilde{\mbox{\boldmath$e$}}_{h}\|_{2}^{2}+\|\tilde{\mbox{\boldmath$e$}}_{h}\|_{4}^{2}).$

Meanwhile, by the discrete Sobolev inequality (2.16) (in Lemma 2.3), we get

(4.27)
$\displaystyle\|\tilde{\mbox{\boldmath$e$}}_{h}\|_{4}\leq\mathcal{C}(\|\tilde{\mbox{\boldmath$e$}}_{h}\|_{2}+\|\tilde{\mbox{\boldmath$e$}}_{h}\|_{2}^{\frac{1}{4}}\cdot\|\nabla_{h}\tilde{\mbox{\boldmath$e$}}_{h}\|_{2}^{\frac{3}{4}}),$

so that

(4.28) $\displaystyle
3(32\delta^{-1}+2)(3M+1)^{2}\|\tilde{\mbox{\boldmath$e$}}_{h}\|_{4}^{2}$
$\displaystyle\leq$ $\displaystyle
3(32\delta^{-1}+2)(3M+1)^{2}\mathcal{C}(\|\tilde{\mbox{\boldmath$e$}}_{h}\|_{2}^{2}+\|\tilde{\mbox{\boldmath$e$}}_{h}\|_{2}^{\frac{1}{2}}\cdot\|\nabla_{h}\tilde{\mbox{\boldmath$e$}}_{h}\|_{2}^{\frac{3}{2}})$
$\displaystyle\leq$
$\displaystyle{\mathcal{C}}_{\delta}\|\tilde{\mbox{\boldmath$e$}}_{h}\|_{2}^{2}+\frac{\delta}{2}\|\nabla_{h}\tilde{\mbox{\boldmath$e$}}_{h}\|_{2}^{2},\quad\forall\delta>0,$

in which the Young’s inequality has been applied in the last step. Going back
(4.26), we obtain

(4.29) $\displaystyle\|\nabla_{h}\mbox{\boldmath$e$}_{h}\|_{2}\leq$
$\displaystyle(1+\frac{\delta}{2})\|\nabla_{h}\tilde{\mbox{\boldmath$e$}}_{h}\|_{2}+\frac{\delta}{2}\|\nabla_{h}\tilde{\mbox{\boldmath$e$}}_{h}\|_{2}+{\mathcal{C}}_{\delta}\|\tilde{\mbox{\boldmath$e$}}_{h}\|_{2}\leq(1+\delta)\|\nabla_{h}\tilde{\mbox{\boldmath$e$}}_{h}\|_{2}+{\mathcal{C}}_{\delta}\|\tilde{\mbox{\boldmath$e$}}_{h}\|_{2},$

provided that $k$ is sufficiently small. This finishes the proof of Lemma 2.4.
∎

## Appendix B. Proof of 2.5

###### Proof.

Here $\underline{\bm{m}}_{h}^{(1)}$ and $\underline{\bm{m}}_{h}^{(2)}$ serve
as the exact solution at different time steps and
$\tilde{\mbox{\boldmath$m$}}_{h}^{(1)},\tilde{\mbox{\boldmath$m$}}_{h}^{(2)}\in
X$ are the corresponding numerical solutions.

The $\|\cdot\|_{\infty}$ bounds for the error functions
$\tilde{\mbox{\boldmath$e$}}_{h}^{(q)}$ can be derived, with the help of the
a-priori estimate (2.20):

(4.30)
$\displaystyle\|\tilde{\mbox{\boldmath$e$}}_{h}^{(q)}\|_{\infty}\leq\gamma
h^{-\frac{1}{2}}(\|\tilde{\mbox{\boldmath$e$}}_{h}^{(q)}\|_{2}+\|\nabla_{h}\tilde{\mbox{\boldmath$e$}}_{h}^{(q)}\|_{2})\leq\gamma
h^{-\frac{1}{2}}\cdot{k}^{\frac{11}{8}}\leq\frac{1}{4}{k}^{\frac{3}{4}},$
$\displaystyle\mbox{so
that}\,\,1-\frac{1}{4}{k}^{\frac{3}{4}}\leq|\tilde{\mbox{\boldmath$m$}}_{h}^{(q)}|\leq
1+\frac{1}{4}{k}^{\frac{3}{4}},\,\,\mbox{at a point-wise level, \, $q=1,2$}.$

For the numerical error functions at different time steps, the decomposition
(4.5) is still valid:

(4.31)
$\tilde{\mbox{\boldmath$e$}}_{h}^{(q)}=\mbox{\boldmath$e$}_{h}^{(q)}+\tilde{\mbox{\boldmath$e$}}_{h,c}^{(q)},\quad\tilde{\mbox{\boldmath$e$}}_{h,c}^{(q)}:=\frac{\tilde{\mbox{\boldmath$m$}}_{h}^{(q)}}{|\tilde{\mbox{\boldmath$m$}}_{h}^{(q)}|}(1-|\tilde{\mbox{\boldmath$m$}}_{h}^{(q)}|).$

At a fixed grid point $(\hat{x}_{i},\hat{y}_{j},\hat{z}_{k})$, we look at the
triangle formed by the vectors: $\tilde{\mbox{\boldmath$m$}}_{h}^{(q)}$,
$\tilde{\mbox{\boldmath$e$}}_{h}^{(q)}$, $\underline{\bm{m}}_{h}^{(q)}$, and
denote the angle between $\tilde{\mbox{\boldmath$m$}}_{h}^{(q))}$ and
$\underline{\bm{m}}_{h}^{(q)}$ as $\theta_{q}$. In turn, the estimate (4.6) is
laid for each $\theta_{q}$:

(4.32) $\displaystyle
1-\frac{1}{4}{k}^{\frac{3}{4}}\leq|\tilde{\mbox{\boldmath$m$}}_{h}^{(q)}|\leq
1+\frac{1}{4}{k}^{\frac{3}{4}},\,\,\,|\tilde{\mbox{\boldmath$e$}}_{h}|\leq\frac{1}{4}{k}^{\frac{3}{4}},\,\,\,|\underline{\bm{m}}_{h}^{(q)}|\equiv
1,$ $\displaystyle 0\leq\sin\theta_{j}\leq\frac{1}{4}{k}^{\frac{3}{4}},\quad
q=1,2.$

Similarly, in the triangle formed by the vectors:
$\mbox{\boldmath$m$}_{h}^{(q)}$, $\mbox{\boldmath$e$}_{h}^{(q)}$,
$\underline{\bm{m}}_{h}^{(q)}$, two sides have equal lengths:
$|\mbox{\boldmath$m$}_{h}^{(q)}|=|\underline{\bm{m}}_{h}^{(q)}|=1$, and the
angle between $\mbox{\boldmath$m$}_{h}^{(q)}$ and
$\underline{\bm{m}}_{h}^{(q)}$ is exactly $\theta_{q}$. In turn, the angle
between vectors $\mbox{\boldmath$m$}_{h}^{(q)}$ and
$\mbox{\boldmath$e$}_{h}^{(q)}$, as well as the angle between vectors
$\underline{\bm{m}}_{h}^{(q)}$ and $\mbox{\boldmath$e$}_{h}^{(q)}$ is given by
$\varphi^{(q)}=\frac{\pi}{2}-\frac{\theta_{q}}{2}$.

Meanwhile, we denote the angle between $\underline{\bm{m}}_{h}^{(1)}$ and
$\underline{\bm{m}}_{h}^{(2)}$ as $\theta^{*}$. By the a-priori assumption and
the fact that
$|\underline{\bm{m}}_{h}^{(1)}|=|\underline{\bm{m}}_{h}^{(2)}|=1$, we have an
estimate for $\theta^{*}$:

(4.33) $2\sin\frac{\theta^{*}}{2}\leq\frac{1}{4}{k}^{\frac{7}{8}}.$

Subsequently, we denote the angle between
$\tilde{\mbox{\boldmath$m$}}_{h}^{(1)}$ and $\mbox{\boldmath$e$}_{h}^{(2)}$ as
$\phi^{*}$. By the above analyses, we see that

(4.34)
$\phi^{*}=\varphi^{(2)}\pm(\theta_{1}+\theta^{*})=\frac{\pi}{2}-\frac{\theta_{2}}{2}\pm(\theta_{1}+\theta^{*}).$

In either case,
$\phi^{*}=\frac{\pi}{2}-\frac{\theta_{2}}{2}+(\theta_{1}+\theta^{*})$ or
$\phi^{*}=\frac{\pi}{2}-\frac{\theta_{2}}{2}-(\theta_{1}+\theta^{*})$, the
following estimate is valid:

(4.35) $\displaystyle|\cos\phi^{*}|=$
$\displaystyle|\sin(\frac{\theta_{2}}{2}\pm(\theta_{1}+\theta^{*}))|\leq\sin(\frac{\theta_{2}}{2}+\theta_{1}+\theta^{*})$
$\displaystyle\leq$
$\displaystyle\sin\frac{\theta_{2}}{2}+\sin\theta_{1}+\sin\theta^{*}\leq\sin\theta_{2}+\sin\theta_{1}+\sin\theta^{*}$
$\displaystyle\leq$
$\displaystyle\frac{1}{4}{k}^{\frac{3}{4}}+\frac{1}{4}{k}^{\frac{3}{4}}+\frac{1}{4}{k}^{\frac{7}{8}}\leq\frac{3}{4}{k}^{\frac{3}{4}},$

in which we have used the triangular inequality,
$\sin(s_{1}+s_{2}+s_{3})\leq\sin s_{1}+\sin s_{2}+\sin s_{3}$, provided that
$s_{1},s_{2},s_{3}>0$ are sufficiently small. As a consequence, the definition
of the point-wise inner product implies the following estimate

(4.36)
$\displaystyle\Big{|}(\tilde{\mbox{\boldmath$e$}}_{h}^{(1)}-\mbox{\boldmath$e$}_{h}^{(1)})\cdot\mbox{\boldmath$e$}_{h}^{(2)}\Big{|}=$
$\displaystyle|\tilde{\mbox{\boldmath$e$}}_{h}^{(1)}-\mbox{\boldmath$e$}_{h}^{(1)}|\cdot|\tilde{\mbox{\boldmath$e$}}_{h}^{(2)}|\cdot|\cos\phi^{*}|$
$\displaystyle\leq$
$\displaystyle|\tilde{\mbox{\boldmath$e$}}_{h}^{(1)}-\mbox{\boldmath$e$}_{h}^{(1)}|\cdot|\tilde{\mbox{\boldmath$e$}}_{h}^{(2)}|\cdot\frac{3}{4}{k}^{\frac{3}{4}}$
$\displaystyle\leq$
$\displaystyle{k}^{\frac{1}{4}}|\tilde{\mbox{\boldmath$e$}}_{h}^{(1)}-\mbox{\boldmath$e$}_{h}^{(1)}|^{2}+{k}^{\frac{4}{5}}|\tilde{\mbox{\boldmath$e$}}_{h}^{(2)}|^{2}.$

Again, this estimate is valid at a point-wise level for
$\mbox{\boldmath$x$}\in\Omega_{h}^{0}$. Therefore, a summation in space leads
to the first inequality (2.22). ∎

## Acknowledgments

This work is supported in part by the grants NSFC 11971021 (J. Chen), NSF
DMS-2012669 (C. Wang), NSFC 11771036 and 12171041 (Y. Cai).

## References