#  Unconditionally optimal convergence of an energy-conserving and linearly
implicit scheme for nonlinear wave equations††thanks: This work is supported
in part by the NSFC (Grant Nos. 11771128, 11871106, 11871092, 11926356) and
NSAF (Grant No. U1930402).

Waixiang Cao  School of Mathematical Sciences, Beijing Normal University,
Beijing 100875, China<EMAIL_ADDRESS>Dongfang Li School of Mathematics and
Statistics, Huazhong University of Science and Technology, Wuhan 430074, China
(dfli@hust.edu.cn); and Hubei Key Laboratory of Engineering Modeling and
Scientific Computing, Huazhong University of Science and Technology, Wuhan
430074, China Zhimin Zhang Beijing Computational Science Research Center,
Beijing 100193, China (zmzhang@csrc.ac.cn); and Department of Mathematics,
Wayne State University, Detroit, MI 48202, USA<EMAIL_ADDRESS>
###### Abstract

In this paper, we present and analyze an energy-conserving and linearly
implicit scheme for solving the nonlinear wave equations. Optimal error
estimates in time and superconvergent error estimates in space are established
without time-step dependent on the spatial mesh size. The key is to estimate
directly the solution bounds in the $H^{2}$-norm for both the nonlinear wave
equation and the corresponding fully discrete scheme, while the previous
investigations rely on the temporal-spatial error splitting approach.
Numerical examples are presented to confirm energy-conserving properties,
unconditional convergence, and optimal error estimates, respectively, of the
proposed fully discrete schemes.

Keywords: scalar auxiliary variable, wave equations, stability, error
estimate, superconvergence

AMS subject classifications: 65M60,65M12, 65N12

##  1 Introduction

We present an energy-conserving and linearly implicit scheme as well as the
unconditionally optimal error estimates for solving the following wave
equation

$\displaystyle\begin{aligned} &u_{tt}=\triangle u-\lambda u-F^{\prime}(u),\ \
&&({\bf x},t)\in\Omega\times(0,T],\\\ &u({\bf x},0)=u_{0}({\bf x}),\ \
u_{t}({\bf x},0)=u_{1}({\bf x}),\ \ &&{\bf x}\in\Omega\end{aligned}$ (1.1)

with the periodic boundary condition, where $\lambda\geq 0$ is a constant,
$\Omega$ is a polygonal or polyhedral domain in ${\mathbb{R}}^{d}(d=2,3)$,
$u_{0}$ and $u_{1}$ are sufficiently smooth, and $F\in C^{2}(\mathbb{R})$ is
the nonlinear potential. For simplicity, we assume that $\Omega$ is a
rectangular or cubic domain. Nonlinear wave equations are widely used to model
plenty of complicated natural phenomena in variety of scientific fields [3,
12, 13, 34]. In the past several decades, it has been one of the hot spots in
the numerical analysis of different schemes for the equations [4, 9, 16, 33,
35].

There are many papers that consider error analysis of the fully discrete
schemes for the nonlinear problems under the following assumption (e.g., [8,
10, 11, 17]),

$|F^{\prime}(u^{n})-F^{\prime}(u_{h}^{n})|\leq L|u^{n}-u_{h}^{n}|,$ (1.2)

where $u^{n}$ and $u_{h}^{n}$ are respectively the theoretical and numerical
solutions, and $L>0$ is the Lipschitz coefficient. A classical model
satisfying (1.2) is the Sine-Gorden equation, whose nonlinear term is
$\sin(u)$. However, as pointed in [16], the assumption (1.2) is not the
typical behavior of the general nonlinear wave equations and thus its range of
the actual applicability is limited.

In order to deal with the non-Lipschitz nonlinearity, one common way is to
impose a priori boundedness of the numerical approximations $u_{h}^{n}$. In
classical finite element analysis, the numerical solutions in the maximum norm
are usually estimated by

$\displaystyle\|u^{n}_{h}\|_{L^{\infty}}$ $\displaystyle\leq$
$\displaystyle\|R_{h}u^{n}\|_{L^{\infty}}+\|R_{h}u^{n}-u^{n}_{h}\|_{L^{\infty}}$
(1.3) $\displaystyle\leq$
$\displaystyle\|R_{h}u^{n}\|_{L^{\infty}}+Ch^{-d/2}\|R_{h}u^{n}-u^{n}_{h}\|_{L^{2}}$
$\displaystyle\leq$
$\displaystyle\|R_{h}u^{n}\|_{L^{\infty}}+Ch^{-d/2}(\tau^{p}+h^{r+1}),$

where $R_{h}$ is the projection operator, $r+1$ and $p$ are convergence orders
in spatial and temporal directions, respectively. Consequently, a time-step
restriction $\tau=O(h^{\frac{d}{2p}})$ is needed in (1.3)(e.g., [7, 14, 15,
26, 27]). Such time-step restriction appears widely in the numerical analysis
but is not always necessary in actual applications.

Unconditional convergence means that the established error bound is valid
without the above mentioned time-step restriction. To achieve the
unconditional convergence, a temporal-spatial error splitting approach is
presented recently [19, 20, 21, 23, 24]. The key idea of the approach is to
introduce a time discrete system, whose solution is denoted by $U^{n}$. Then,
one can get the following error estimates

$\|U^{n}\|_{H^{2}}\leq
C~{}~{}~{}\textrm{and}~{}~{}~{}\|R_{h}U^{n}-u_{h}^{n}\|_{L^{2}}\leq Ch^{2}.$

The boundedness of the numerical solutions is obtained by

$\displaystyle\|u^{n}_{h}\|_{L^{\infty}}$ $\displaystyle\leq$
$\displaystyle\|R_{h}U^{n}-u^{n}_{h}\|_{L^{\infty}}+\|R_{h}U^{n}-{\cal
I}_{h}U^{n}\|_{L^{\infty}}+\|{\cal I}_{h}U^{n}\|_{L^{\infty}}$ (1.4)
$\displaystyle\leq$ $\displaystyle
Ch^{-d/2}(\|R_{h}U^{n}-u^{n}_{h}\|_{L^{2}}+\|R_{h}U^{n}-{\cal
I}_{h}U^{n}\|_{L^{2}})+\|U^{n}\|_{L^{\infty}}$ $\displaystyle\leq$
$\displaystyle Ch^{2-d/2}+\|U^{n}\|_{H^{2}}.$

Here ${\cal I}_{h}U^{n}$ denotes the interponation function of $U^{n}$. It
implies the numerical solutions is bounded if the temporal and spatial step-
sizes are sufficiently small, respectively. Then, the error estimates can be
proved following the usual way. In spite of the interesting and instructive
work, additional error estimates in different norms are required in the proof,
and so far, most unconditional convergence results are focused on nonlinear
parabolic problems.

Nonlinear wave equations (1.1) have several remarkable features. First, the
models (1.1) are energy-conserving, i.e.,

$E(t)=\int_{\Omega}u_{t}^{2}+|\nabla u|^{2}+\lambda|u|^{2}+2F(u)dx=E(0).$

Second, the typical nonlinear terms are non-Lipschitz continuous. Third, the
solutions have different regularities. A natural question is whether we can
develop some effective unconditional convergence numerical schemes for
nonlinear wave equations, taking all the remarkable features into account.

In the present paper, we present an energy-conserving and linearly implicit
scheme for solving the nonlinear wave equations (1.1). The scheme is of order
$2$ in the temporal directions and no additional initial iterations are
required. The construction of the scheme is based on recently-developed scalar
auxiliary variable (SAV) approach combined with finite element methods,
classical Crank-Nicolson methods, and extrapolation approximation. We show
that our fully-discrete schemes conserve the energy and converge without any
temporal step-size restrictions dependent on the spatial step-size. Unlike the
previous temporal-spatial error splitting approach, we estimate the solution
directly in the following procedure: (1) obtain the bounds in the $H^{2}$-norm
of the solutions for both the nonlinear wave equations and the corresponding
fully discrete schemes; (2) establish the bound for numerical solutions by
applying the embedding inequality; (3) obtain the unconditionally optimal
error estimates in time and superconvergent error estimates in space.

We remark that the key to construct the energy-conserving schemes is the SAV
idea, which has been applied successfully to the gradient flows [1, 28, 29,
30]. Very recently, the idea was introduced to develop energy-conserving
schemesfor the conservative laws [5, 6, 22, 25]. However, much attention has
been paid to the stability and energy-conserving properties, and no
unconditional convergence results of fully discrete SAV schemes for nonlinear
wave equations are found in the literature. This is the main motivation and
contribution of the present study.

The rest of paper is organized as follows. In Section 2, we propose a fully
discrete scheme for the nonlinear wave equation (1.1). In Section 3, we
present a detailed proof to show the energy-conserving properties and
unconditional convergence for the temporal discretization. Error estimates for
the fully discrete solution is established in Section 4, where we prove that
the approximation error is unconditionally optimal in time and superconvergent
in space (under the $H^{1}$ norm). In Section 5, we present several numerical
examples to confirm the theoretical results. Finally, conclusions are
presented in Section 6.

##  2 The linearly implicit method

In this section, we present a fully discrete numerical scheme, which preserves
the discrete energy.

Suppose $E_{1}(u)=\int_{\Omega}F(u)d{\bf x}\geq-c_{0}$ for some $c_{0}>0$,
i.e., it is bounded from below, and let $C_{0}>c_{0}$ so that
$E_{1}(u)+C_{0}>0$. We introduce the following scalar auxiliary variable (SAV)

$r(t)=\sqrt{E(u)},\ \ E(u)=\int_{\Omega}F(u)d{\bf x}+C_{0},$

and rewrite (1.1) as

$\displaystyle\begin{aligned} &u_{t}=v,&\\\ &v_{t}=\triangle u-\lambda
u-\frac{r}{\sqrt{E(u)}}f(u),&\\\
&r_{t}=\frac{1}{2\sqrt{E(u)}}\int_{\Omega}f(u)u_{t}d{\bf x},&\end{aligned}$
(2.1)

where $f(u)=F^{\prime}(u)$.

Let ${\mathcal{T}_{h}}$ be the usual regular triangulation of the polygonal
domain $\Omega$. Denote by $h_{T}$ the meshsize of $\mathcal{T}_{h}$, where
$h_{T}$ is the diameter of the element ${T}\in\mathcal{T}_{h}$, and
$h=\max_{T\in\mathcal{T}_{h}}$. Let $V_{h}$ be the classical finite-
dimensional subspace of $H^{1}(\Omega)$, which consists of the usual
continuous piecewise polynomials of degree $k$ ($k\geq 1$) on
$\mathcal{T}_{h}$. That is,

$V_{h}=\\{v\in C^{0}(\Omega):v|_{T}\in{\mathbb{P}}_{k},\ \forall T\in{\cal
T}_{h}\\},$

where ${\mathbb{P}}_{k}$ denotes the space of polynomials of degree no more
than $k$.

Let $\tau=\frac{T}{N}$ with $N$ a given integer and $t_{n}=n\tau$,
$n=0,1,\cdots,N$. Denote

$u^{n}=u({{\bf x}},t_{n}),\ \ v^{n}=v({{\bf x}},t_{n}),\ \ r^{n}=r(t_{n}).$

For any sequence of the functions $\\{f^{n}\\}_{n=0}^{N}$, we define for all
$n,n=0,\ldots N-1$

$D_{\tau}f^{n+1}:=\frac{Df^{n+1}}{\tau}=\frac{f^{n+1}-f^{n}}{\tau},\ \
\tilde{f}^{n+\frac{1}{2}}:=\frac{1}{2}(3f^{n}-f^{n-1}),\ \
\hat{f}^{n+\frac{1}{2}}:=\frac{f^{n+1}+f^{n}}{2}.$

Note that for $n=0$, we denote by $\tilde{f}^{\frac{1}{2}}=f^{0}$.

To design an energy-conserving and linearly implicit numerical scheme, which
is easy to implement and efficient, we consider the following fully discrete
Crank-Nicolson Galerkin SAV method: Find $u^{n+1}_{h}\in V_{h},v^{n+1}_{h}\in
V_{h},r_{h}^{n+1}\in\mathbb{R}$ for $n=0,\ldots N-1$ such that for all
$(w_{h},\zeta_{h})\in V_{h}\times V_{h}$

$\displaystyle\begin{split}&\big{(}D_{\tau}u_{h}^{n+1},w_{h}\big{)}=\big{(}\hat{v}_{h}^{n+\frac{1}{2}},w_{h}\big{)},&&\\\
&\big{(}D_{\tau}v_{h}^{n+1},\zeta_{h}\big{)}=-\big{(}\nabla\hat{u}_{h}^{n+\frac{1}{2}},\nabla\zeta_{h}\big{)}-\big{(}\lambda\hat{u}_{h}^{n+\frac{1}{2}},\zeta_{h}\big{)}-\big{(}\frac{\hat{r}_{h}^{n+\frac{1}{2}}}{\sqrt{E(\tilde{u}_{h}^{n+\frac{1}{2}})}}f(\tilde{u}_{h}^{n+\frac{1}{2}}),\zeta_{h}\big{)},&&\\\
&{r_{h}^{n+1}-r_{h}^{n}}=\frac{1}{2\sqrt{E(\tilde{u}_{h}^{n+\frac{1}{2}})}}\int_{\Omega}f(\tilde{u}_{h}^{n+\frac{1}{2}})({u_{h}^{n+1}-u_{h}^{n}})d{\bf
x},&&\end{split}$ (2.2)

where $(u,v)=\int_{\Omega}u({\bf x})v({\bf x})d{\bf x}$,
$f(\tilde{u}_{h}^{\frac{1}{2}})=f(u_{h}^{0})$, and initial values are chosen
as

${(u_{h}^{0},v_{h}^{0},r_{h}^{0})=(R_{h}u_{0},R_{h}u_{1},\sqrt{E(u_{0})})}.$

Here $R_{h}u_{0}$ is the Ritz projection of $u_{0}$, which will be defined
later.

Equivalently, we rewrite the above scheme (2.2) into the following linear form

$((4I-\tau^{2}{\triangle_{h}}{+}\tau^{2}\lambda
I)u_{h}^{n+1},w_{h})+\frac{\tau^{2}}{2}(u_{h}^{n+1},b_{1})\big{(}b_{1},w_{h}\big{)}=(g,w_{h})+\frac{\tau^{2}}{2}(u_{h}^{n},b_{1})\big{(}b_{1},w_{h}\big{)}$
(2.3)

for all $w_{h}\in V_{h}$, where $(\triangle_{h}u_{h},v_{h}):=-(\nabla
u_{h},\nabla v_{h})$, and

$b_{1}=\frac{f(\tilde{u}_{h}^{n+\frac{1}{2}})}{\sqrt{E(\tilde{u}_{h}^{n+\frac{1}{2}})}},\
\ g=(4I{+}\tau^{2}{\triangle_{h}}-\tau^{2}\lambda I)u_{h}^{n}+4\tau
v_{h}^{n}-2\tau^{2}r^{n}_{h}b_{1}.$

Choosing $w_{h}$ in (2.3) to be the basis function of $V_{h}$ leads to a
linear equation of the form

$A{\bf u}^{n+1}+({\bf u}^{n+1},{\bf b}_{1}){\bf b}_{2}={\bf g}$

for some matrix $A$ and vectors ${\bf b}_{1},{\bf b}_{2},{\bf g}$. By taking
the inner product with ${\bf b}_{1}$ in the above equation, we obtain $({\bf
u}^{n+1},{\bf b}_{1})$ and then derive ${\bf u}^{n+1}$. Hence the scheme is
easy to implement and very efficient. We also refer to [29, 30] for more
detailed information.

##  3 Unconditionally energy preservation and convergence for the temporal
discretization

In this section, we shall prove that the Galerkin SAV method (2.2) preserves
the energy unconditionally. Moreover, we establish the convergence analysis of
the SAV approach with minimum assumptions.

We begin with the energy preservation property of the Galerkin SAV approach.
We define the energy

$E^{n}=\sqrt{\frac{1}{2}\big{(}\|v_{h}^{n}\|^{2}+\|\nabla
u_{h}^{n}\|^{2}+\lambda\|u_{h}^{n}\|^{2}\big{)}+(r_{h}^{n})^{2}},\ \ 1\leq
n\leq N.$

Here $\|u\|^{2}=(u,u)=\|u||_{L^{2}}^{2}$. Taking
$(w_{h},\zeta_{h})=(v_{h}^{n+1}-v_{h}^{n},u_{h}^{n+1}-u_{h}^{n})$ and
multiplying the third equation of (2.2) by ${r_{h}^{n+1}+r_{h}^{n}}$, we
derive

$\displaystyle\frac{1}{2}\big{(}\|v_{h}^{n+1}\|^{2}-\|v_{h}^{n}\|^{2}+\|\nabla
u_{h}^{n+1}\|^{2}-\|\nabla
u_{h}^{n}\|^{2}+\lambda\|u_{h}^{n+1}\|^{2}-\lambda\|u_{h}^{n}\|^{2}\big{)}+(r_{h}^{n+1})^{2}-(r_{h}^{n})^{2}=0.$

Consequently,

$E^{n+1}=E^{n}=E^{0},\ \ \forall n\geq 1.$

Now we consider a time-discrete system of equation:

$\displaystyle\begin{split}&D_{\tau}{U^{n+1}}=\hat{V}^{n+\frac{1}{2}},&\\\
&D_{\tau}{V^{n+1}}=\triangle\hat{U}^{n+\frac{1}{2}}-\lambda\hat{U}^{n+\frac{1}{2}}-\frac{\hat{R}^{n+\frac{1}{2}}}{\sqrt{E(\tilde{U}^{n+\frac{1}{2}})}}f(\tilde{U}^{n+\frac{1}{2}}),&\\\
&{R^{n+1}-R^{n}}=\frac{1}{2\sqrt{E(\tilde{U}^{n+\frac{1}{2}})}}\int_{\Omega}f(\tilde{U}^{n+\frac{1}{2}})({U^{n+1}-U^{n}})d{\bf
x},&\end{split}$ (3.1)

subject to the periodic boundary condition and the following initial
conditions

$\displaystyle U^{0}({\bf x})=u_{0}({\bf x}),\ \ V^{0}({\bf x})=u_{1}({\bf
x}).$

As we may observe, the numerical solution $(u_{h}^{n},v_{h}^{n},r_{h}^{n})$
can be viewed as the Galerkin approximation of the above time-discrete system
of equation. To study the convergence of the temporal discretization (3.1), we
need some preliminaries.

First, for the simplicity of notations, throughout this paper, we denote by
$C$ a generic positive constant, which depend solely upon the physical
parameters of the problem and independent of $\tau,h,n$, and it is not
necessary to be the same at every appearance. We adopt the usual notations for
Sobolev spaces, e.g., $W^{m,p}(I)$ on sub-domain $I\in\Omega$ equipped with
the norm $\|\cdot\|_{W^{m,p},I}$ and semi-norm $|\cdot|_{W^{m,p},I}$. We omit
the index $I$ when $I=\Omega$. Especially, when $p=2$, we set
$W^{m,p}(I)=H^{m}(I)$ and $\|\cdot\|_{W^{m,p},I}=\|\cdot\|_{H^{m},I}$ and
$|\cdot|_{W^{m,p},I}=|\cdot|_{H^{m},I}.$ The notation $\alpha\lesssim\beta$
implies that $\alpha$ is bounded by $\beta$ multiplied by a constant
independent of $\tau,h,n$.

Second, we would like to present a Gronwall-type inequalities, which play
important role in our later convergence analysis and error estimates.

###### Lemma 1

( [18] ) Let $\tau$, $B$ and $a_{k}$, $b_{k}$, $c_{k}$, $\gamma_{k}$, for
integers $k>0$, be nonnegative numbers such that

$a_{n}+\tau\sum_{k=0}^{n}b_{k}\leq\tau\sum_{k=0}^{n}\gamma_{k}a_{k}+\tau\sum_{k=0}^{n}c_{k}+B,~{}~{}~{}~{}~{}~{}\mbox{for}~{}~{}n\geq
0.$

Suppose that $\tau\gamma_{k}<1$, for all $k$, and set
$\sigma_{k}=(1-\tau\gamma_{k})^{-1}$. Then,

$a_{n}+\tau\sum_{k=0}^{n}b_{k}\leq\Big{(}\tau\sum_{k=0}^{n}c_{k}+B\Big{)}\exp\Big{(}\tau\sum_{k=0}^{n}\gamma_{k}\sigma_{k}\Big{)}.$

###### Lemma 2

( [2] ) Let $I=[a,b]$ and $\alpha(t),\beta(t),u(t)\in C^{0}(I)$. Suppose
$\beta(t)\geq 0$ and

$u(t)\leq\alpha(t)+\int_{a}^{t}\beta(s)u(s)ds,\ \forall t\in I.$

Then

$u(t)\leq\alpha(t)+\int_{a}^{t}\alpha(s)\beta(s)e^{\int_{s}^{t}\beta(r)dr}ds,\
\ \forall t\in I.$

Now we are ready to study the convergence of the solution of (3.1). Taking the
inner product of the first two equations with $V^{n+1}-V^{n},U^{n+1}-U^{n}$
and multiplying the third equation of (3.1) by ${R^{n+1}+R^{n}}$, we derive

$\displaystyle\\!\frac{1}{2}\big{(}\|V^{n+1}\|^{2}\\!-\\!\|V^{n}\|^{2}\\!+\\!\|\nabla
U^{n+1}\|^{2}-\|\nabla
U^{n}\|^{2}\\!+\\!\lambda\|U^{n+1}\|^{2}\\!-\\!\lambda\|U^{n}\|^{2}\big{)}\\!+\\!(R^{n+1})^{2}\\!-\\!(R^{n})^{2}=0,$

which indicates that

$\|V^{n}\|+\|U^{n}\|_{H^{1}}+|R^{n}|\lesssim 1,\ \ \forall 1\leq n\leq N.$

As pointed out in [28], energy stable is not sufficient for the convergence
which typically needs bounds in higher norms. Following the idea in [28], our
convergence analysis is along this line: we first start from the energy
preservation to derive the error bounds in higher norms, (i.e., the $H^{2}$
estimates) for the solution $U^{n}$, and thus get the $L^{\infty}$ for $U^{n}$
thanks to the embedding theory, and then we use the bounds in $H^{2}$ norms to
show that the numerical solution $U^{n}$ converges to the exact solution
$u^{n}$ in some suitable norms as $\tau$ tends to zero. To this end, we need
to the bounds in $H^{2}$ norm of the PDE system (1.1). The error bounds for
the solution of (3.1) is similar to those of the PDE system.

Note that most of the convergence and error analysis for linearly implicit are
based on the so called Lipschitz assumption, i.e.,

$|F^{\prime}(u_{1})-F^{\prime}(u_{2})|\leq L|u_{1}-u_{2}|,\ \ \forall
u_{1},u_{2}.$ (3.2)

The above assumption greatly limits its range of applicability. Following the
basic idea of [28], we adopt the following assumption instead of the Lipschitz
assumption in our convergence analysis:

$\displaystyle|f^{\prime}(x)|<C(|x|^{p}+1),\ \ p\geq 0\ {\rm if}\ n=1,2;\
0<p<4\ {\rm if}\ n=3,$ (3.3)
$\displaystyle|f^{\prime\prime}(x)|<C(|x|^{p}+1),\ \ p\geq 0\ {\rm if}\
n=1,2;\ 0<p<3\ {\rm if}\ n=3.$ (3.4)

It has been proved in [28] that if $f(u)$ satisfies the conditions
(3.3)-(3.4), there holds for some $\sigma$, where $0\leq\sigma<1$, such that

$\displaystyle\|f^{\prime\prime}(u)\|_{L^{\infty}}+\|f^{\prime}(u)\|_{L^{\infty}}\leq
C(1+\|\nabla\triangle u\|^{\sigma})\leq\epsilon\|\nabla\triangle
u\|+C_{\epsilon}$ (3.5)

for any $\epsilon>0$ with $C_{\epsilon}$ a constant depending on $\epsilon$.

We present the following estimates for the exact solution of (1.1).

###### Proposition 1

Assume that $u$ is the solution of (1.1), and $u_{0}\in H^{3},u_{1}\in H^{2}$
and (3.3)-(3.4) holds. Then for any $T>0$,

$(\|\triangle u\|+\|\triangle u_{t}\|+\|\nabla\triangle u\|)(T)\lesssim 1.$

Proof. First, multiplying $u_{t}$ on both sides of (1.1) and using the
integration by parts yields

$\frac{d}{dt}(\|u_{t}\|^{2}+\|\nabla
u\|^{2}+\lambda\|u\|^{2}+{2}\int_{\Omega}F(u))=0,$

which indicates that

$\|u_{t}\|^{2}+\|u\|_{H^{1}}^{2}+\int_{\Omega}F(u)\lesssim 1.$ (3.6)

On the other hand, we multiply $\triangle^{2}u_{t}$ on both sides of (1.1) and
again use the integration by parts to obtain

$\displaystyle\begin{split}\frac{1}{2}\frac{d}{dt}\left(\|\triangle
u_{t}\|^{2}+\|\nabla\triangle u\|^{2}+\lambda\|\triangle
u\|^{2}\right)&={-(\triangle f(u),\triangle u_{t}).}&\\\ \end{split}$ (3.7)

Integrating with respect to time between $0$ to $t$ and using Cauchy-Schwarz
inequality yields

$\displaystyle\left(\|\triangle u_{t}\|^{2}+\|\nabla\triangle
u\|^{2}+\lambda\|\triangle u\|^{2}\right)(t)$ $\displaystyle\leq$
$\displaystyle\left(\|\triangle u_{t}\|^{2}+\|\nabla\triangle
u\|^{2}+\lambda\|\triangle u\|^{2}\right)(0)$ (3.8) $\displaystyle+$
$\displaystyle\int_{0}^{t}(\|\triangle u_{t}\|^{2}+\|\triangle
f(u)\|^{2})(t)dt.$

By (3.5) and the identity

$\triangle f(u)=f^{\prime}(u)\triangle u+f^{\prime\prime}(u)|\nabla u|^{2},$

we have for all $0\leq\delta<1$,

$\displaystyle\|\triangle f(u)\|$ $\displaystyle\leq$
$\displaystyle\|f^{\prime\prime}(u)\|_{L^{\infty}}\|\nabla
u\|_{L^{4}}^{2}+\|f^{\prime}(u)\|_{L^{\infty}}\|\triangle u\|$
$\displaystyle\leq$ $\displaystyle
C(\|f^{\prime\prime}(u)\|_{L^{\infty}}+\|f^{\prime}(u)\|_{L^{\infty}})(\|\nabla
u\|_{L^{4}}^{2}+\|\triangle u\|)$ $\displaystyle\leq$ $\displaystyle
C(1+\|\nabla\triangle u\|^{\delta})({\|\nabla u\|_{L^{4}}^{2}}+\|\triangle
u\|).$

As for the term $\|\nabla u\|_{L^{4}}$, we use the Sobolev embedding theory
and interpolation inequality about the spaces $H^{s}$ (see, e.g., [31]) and
then obtain

$\|\nabla u\|_{L^{4}}\leq C\|\nabla u\|_{H^{d/4}}\leq C\|\nabla
u\|^{1-d/8}\|\nabla\triangle u\|^{d/8}\leq C\|\nabla\triangle u\|^{d/8}.$

Moreover, by using the integration by parts and (3.6),

$\|\triangle u\|^{2}\leq C\|\nabla\triangle u\|\|\nabla u\|\leq
C\|\nabla\triangle u\|.$

Consequently,

$\|\triangle f(u)\|^{2}\leq C(1+\|\nabla\triangle u\|^{2\delta})(\|\nabla
u\|_{L^{4}}^{4}+\|\triangle u\|^{2})\leq C(1+\|\nabla\triangle u\|^{2}).$

Substituting the above inequality into (3.8) gives

$\displaystyle\begin{split}\left(\|\triangle u_{t}\|^{2}+\|\nabla\triangle
u\|^{2}+\lambda\|\triangle u\|^{2}\right)(t)&\lesssim\left(\|\triangle
u_{t}\|^{2}+\|\nabla\triangle u\|^{2}+\lambda\|\triangle
u\|^{2}\right)(0)+1+&\\\ &+\int_{0}^{t}\left(\|\triangle
u_{t}\|^{2}+\|\nabla\triangle u\|^{2}\right)(t)dt.&\end{split}$ (3.9)

By Gronwall inequality in Lemma 2, there holds

$\displaystyle\left(\|\triangle u_{t}\|^{2}+\|\nabla\triangle
u\|^{2}+\lambda\|\triangle u\|^{2}\right)(t)$
$\displaystyle\lesssim\left(\|\triangle u_{t}\|^{2}+\|\nabla\triangle
u\|^{2}+\lambda\|\triangle u\|^{2}\right)(0)+1\lesssim 1.$

This finishes our proof.

Similar to the proof in the above Proposition, we also have the following
$H^{2}$ estimates for the solution of (3.1).

###### Proposition 2

Assume that $(U^{n},V^{n},R^{n})$ are the solutions of (3.1) and (3.3)-(3.4)
holds. Then

$\|\triangle U^{n}\|+\|\triangle V^{n}\|^{2}+\|\nabla\triangle
U^{n}\|^{2}\lesssim 1.$ (3.10)

Proof. First, we have from the first equation of (3.1) that

$D_{\tau}\nabla\triangle U^{n+1}=\nabla\triangle\hat{V}^{n+{\frac{1}{2}}}.$

Multiplying the above equation with $\nabla(V^{n+1}-V^{n})$ and the second
equation of (3.1) with $\triangle^{2}(U^{n+1}-U^{n})$, and then using the
integration by parts, we obtain

$\displaystyle\|\nabla\triangle U^{n+1}\|^{2}-\|\nabla\triangle
U^{n}\|^{2}+\|\triangle V^{n+1}\|^{2}-\|\triangle
V^{n}\|^{2}+\lambda\|\triangle U^{n+1}\|^{2}-\lambda\|\triangle U^{n}\|^{2}$
$\displaystyle=\frac{{2}\hat{R}^{n+\frac{1}{2}}}{\sqrt{E(\tilde{U}^{n+\frac{1}{2}})}}(\nabla
f(\tilde{U}^{n+\frac{1}{2}}),\nabla\triangle(U^{n+1}-U^{n})).$

Denoting
$f_{1}(\tilde{U}^{n+\frac{1}{2}})=\frac{2\hat{R}^{n+\frac{1}{2}}f(\tilde{U}^{n+\frac{1}{2}})}{\sqrt{E(\tilde{U}^{n+\frac{1}{2}})}}$
and summing up the above equation for all $n$ from $0$ to $m$ yields

$\displaystyle\begin{split}&\|\nabla\triangle U^{m+1}\|^{2}-\|\nabla\triangle
U^{0}\|^{2}+\|\triangle V^{m+1}\|^{2}-\|\triangle
V^{0}\|^{2}+\lambda\|\triangle U^{m+1}\|^{2}-\lambda\|\triangle
U^{0}\|^{2}&\\\ &=\\!(\nabla
f_{1}(\tilde{U}^{m\\!+\\!\frac{1}{2}}),\\!\nabla\triangle\\!U^{m\\!+\\!1}\\!)\\!-\\!(\nabla
f_{1}\\!(\tilde{U}^{\frac{1}{2}}\\!),\\!\nabla\triangle\\!U^{0}\\!)\\!+\\!\sum_{n=1}^{m}\\!(\nabla\\!(f_{1}(\tilde{U}^{n\\!-\\!\frac{1}{2}})\\!-\\!f_{1}(\tilde{U}^{n+\frac{1}{2}})\\!),\nabla\triangle
U^{n})&\\\ &=(\nabla f_{1}(U^{m+1}),\nabla\triangle U^{m+1})-(\nabla
f_{1}(\tilde{U}^{\frac{1}{2}}),\nabla\triangle U^{0})+I_{1}-I,&\end{split}$
(3.11)

where
$f_{1}(U^{m+1})=\frac{2\hat{R}^{m+\frac{1}{2}}f(U^{m+1})}{\sqrt{E(\tilde{U}^{m+\frac{1}{2}})}}$
, and

$I=(\nabla f_{1}(U^{m+1})-\nabla
f_{1}(\tilde{U}^{m+\frac{1}{2}}),\nabla\triangle U^{m+1}),\ \ \
I_{1}=\sum_{n=1}^{m}(\nabla(f_{1}(\tilde{U}^{n-\frac{1}{2}})-f_{1}(\tilde{U}^{n+\frac{1}{2}})),\nabla\triangle
U^{n}).$

Since $R^{n}$ is bounded and $E(U)$ is bounded from below, we have from (3.5)
that,

$\|\nabla f_{1}(U^{m+1})\|\leq
C(1+\|f^{\prime}(U^{m+1})\|_{L^{\infty}})\leq\epsilon\|\nabla\triangle
U^{m+1}\|+C_{\epsilon}.$ (3.12)

Consequently,

$|(\nabla f_{1}(U^{m+1}),\nabla\triangle U^{m+1})-(\nabla
f_{1}(\tilde{U}^{\frac{1}{2}}),\nabla\triangle
U^{0})|\leq\frac{1}{4}\|\nabla\triangle
U^{m+1}\|^{2}+\frac{1}{4}\|\nabla\triangle U^{0}\|^{2}+C.$ (3.13)

On the other hand, we note that

$\displaystyle\begin{split}&\nabla f_{1}(\tilde{U}^{n-\frac{1}{2}})-\nabla
f_{1}(\tilde{U}^{n+\frac{1}{2}})&\\\
&=f^{\prime}_{1}(\tilde{U}^{n-\frac{1}{2}})\nabla\tilde{U}^{n-\frac{1}{2}}-f^{\prime}_{1}(\tilde{U}^{n+\frac{1}{2}})\nabla\tilde{U}^{n+\frac{1}{2}}\\\
&=f^{\prime}_{1}(\tilde{U}^{n-\frac{1}{2}})(\nabla\tilde{U}^{n-\frac{1}{2}}-\nabla\tilde{U}^{n+\frac{1}{2}})+(f^{\prime}_{1}(\tilde{U}^{n-\frac{1}{2}})-f^{\prime}_{1}(\tilde{U}^{n+\frac{1}{2}}))\nabla\tilde{U}^{n+\frac{1}{2}}&\\\
&=\frac{\tau}{2}f^{\prime}_{1}(\tilde{U}^{n-\frac{1}{2}})\nabla(\hat{V}^{n-\frac{3}{2}}-3\hat{V}^{n-\frac{1}{2}})+\frac{\tau}{2}f^{\prime\prime}_{1}(\theta\tilde{U}^{n+\frac{1}{2}}+(1-\theta)\tilde{U}^{n-\frac{1}{2}})(\hat{V}^{n-\frac{3}{2}}-3\hat{V}^{n-\frac{1}{2}})\nabla\tilde{U}^{n+\frac{1}{2}}\end{split}$

for some $\theta\in(0,1)$, where in the last step, we have used the first
equation of (3.1), which yields

$\nabla\tilde{U}^{n-\frac{1}{2}}-\nabla\tilde{U}^{n+\frac{1}{2}}={\frac{1}{2}\nabla(4U^{n-1}-3U^{n}-U^{n-2})=\frac{\tau}{2}\nabla(\hat{V}^{n-\frac{3}{2}}-3\hat{V}^{n-\frac{1}{2}}).}$

By (3.12) and the fact that

$\|V^{n}\|+\|\nabla U^{n}\|\lesssim 1,$

we have

$\displaystyle\|\nabla
f_{1}\\!(\tilde{U}^{n-\frac{1}{2}}\\!)\\!-\\!\nabla\\!f_{1}\\!(\tilde{U}^{n+\frac{1}{2}}\\!)\\!\|^{2}\\!\leq\\!\tau^{2}(\epsilon\|\nabla\triangle\tilde{U}^{n\\!-\\!\frac{1}{2}}\|^{2}\\!+\\!\epsilon\|\nabla\triangle\tilde{U}^{n+\frac{1}{2}}\|^{2}\\!+\\!\|\nabla(\hat{V}^{n-\frac{3}{2}}\\!-3\\!\hat{V}^{n-\frac{1}{2}})\|^{2}\\!+\\!C^{2}),$

and thus

$\displaystyle|I_{1}|\leq C\tau\sum_{n=1}^{m}(\|\nabla\triangle
U^{n}\|^{2}+\|{\nabla}V^{n}\|^{2})+C\tau.$

Similarly, there holds

$\displaystyle\|\nabla f_{1}(\tilde{U}^{m+\frac{1}{2}})-\nabla
f_{1}(U^{m+1})\|^{2}\leq\tau^{2}\epsilon\sum_{n=m-1}^{m+1}(\|\nabla
V^{n}\|^{2})+\|\nabla\triangle U^{n}\|^{2})+C_{\epsilon}\tau^{2}.$

Then

$|I|\leq\tau^{2}\sum_{n=m-1}^{m+1}(\|\nabla
V^{n}\|^{2})+\epsilon\|\nabla\triangle{U}^{n}\|^{2})+C_{\epsilon}\tau^{2}+\frac{1}{4}\|\nabla\triangle
U^{m+1}\|^{2}.$

Substituting (3.13), the estimates of $I_{1}$ and $I$ into (3.11) yields

$\displaystyle\frac{1}{2}\|\nabla\triangle U^{m+1}\|^{2}+\|\triangle
V^{m+1}\|^{2}+\lambda\|\triangle U^{m+1}\|^{2}$ $\displaystyle\leq$
$\displaystyle\|\triangle V^{0}\|^{2}+\|\nabla\triangle
U^{0}\|^{2}+\lambda\|\triangle U^{0}\|^{2}$ $\displaystyle+$ $\displaystyle
C\tau\sum_{n=1}^{m+1}(\|\nabla\triangle U^{n}\|^{2}+\|\triangle
V^{n}\|^{2})+C.$

By Gronwall inequality, we have

$\|\nabla\triangle U^{m+1}\|^{2}+\|\triangle V^{m+1}\|^{2}+\lambda\|\triangle
U^{m+1}\|^{2}\leq\|\triangle V^{0}\|^{2}+\|\nabla\triangle
U^{0}\|^{2}+\lambda\|\triangle U^{0}\|^{2}\leq C.$

In case $\lambda=0$, we note that

$\triangle U^{m+1}=\triangle U^{m}+\frac{\tau}{2}(\triangle
V^{{m}+1}+\triangle V^{{m}}).$

Then

$\|\triangle U^{m+1}\|\leq\|\triangle U^{m}\|+\frac{\tau}{2}\|\triangle
V^{m+1}+\triangle V^{m}\|,$

which yields

$\|\triangle U^{m+1}\|\leq\|\triangle U^{0}\|+C\leq C.$

This finishes the proof of (3.10). The proof is complete.

As a direct consequence of (3.10) and the embedding inequality, we have

$\|U^{n}\|_{\infty}\leq C\|\triangle U^{n}\|\leq C,\ \ \forall n.$ (3.14)

###### Remark 1

Following the same argument as that in [28], we conclude that: Assume that
$u_{0}\in H^{3}$, when $\tau$ tends to zero, we have $U^{n}\rightarrow u^{n}$
strongly in $L^{\infty}(0,T;H^{3-\epsilon}),\ \forall\epsilon>0$, weak-star in
$L^{\infty}(0,T;H^{3})$, $V^{n}\rightarrow v^{n}$ weak-star in
$L^{\infty}(0,T;H^{2})$, and $R^{n}\rightarrow r^{n}$ weak-star in
$L^{\infty}(0,T)$;

###### Theorem 1

Suppose $u$ is the solution of (1.1), satisfying

$\|u_{0}\|_{H^{2}}+\|u\|_{L^{\infty}((0,T),H^{2})}+\|u_{t}\|_{L^{2}((0,T),H^{2})}+\|u_{tt}\|_{L^{2}((0,T),H^{2})}\lesssim
1.$

Then (3.1) admits a unique solution $(U^{n},V^{n},R^{n})$ such that

$\|\ u^{n}-U^{n}\|_{H^{1}}+\|v^{n}-V^{n}\|+|r^{n}-R^{n}|\lesssim\tau^{2}.$

Proof. First we denote

$E^{n}_{u}=u^{n}-U^{n},\ \ E^{n}_{v}=v^{n}-V^{n},\ \ E^{n}_{r}=r^{n}-R^{n},\ \
H(u)=\frac{f(u)}{\sqrt{E(u)}}.$

By taking $t=t_{n+\frac{1}{2}}$ in (2.1) and using (3.1), we get

$\displaystyle\begin{split}&D_{\tau}E_{u}^{n+1}=\hat{E}_{v}^{n+\frac{1}{2}}+T_{1},&&\\\
&D_{\tau}E_{v}^{n+1}=\triangle\hat{E}_{u}^{n+\frac{1}{2}}{-}\lambda\hat{E}_{u}^{n+\frac{1}{2}}-\hat{r}^{n+\frac{1}{2}}H(u^{n+\frac{1}{2}})+\hat{R}^{n+\frac{1}{2}}H(\tilde{U}^{n+\frac{1}{2}})+T_{2},&&\\\
&{E_{r}^{n+1}-E_{r}^{n}}=\frac{1}{2}\int_{\Omega}H(u^{n+\frac{1}{2}})({u^{n+1}-u^{n}})-H(\tilde{U}^{n+\frac{1}{2}})({U^{n+1}-U^{n}})d{\bf
x}+T_{3},&&\end{split}$ (3.15)

where $T_{i}$ denote the truncation errors, i.e.,

$\displaystyle
T_{1}=D_{\tau}u^{n+1}-u^{n+\frac{1}{2}}_{t}+v^{n+\frac{1}{2}}-\hat{v}^{n+\frac{1}{2}},$
$\displaystyle
T_{2}=\triangle(u^{n+\frac{1}{2}}-\hat{u}^{n+\frac{1}{2}})+D_{\tau}v^{n+1}-v^{n+\frac{1}{2}}_{t}{-}\lambda
u^{n+\frac{1}{2}}{+}\lambda\hat{u}^{n+\frac{1}{2}}+(\hat{r}^{n+\frac{1}{2}}-r^{n+\frac{1}{2}})H(u^{n+\frac{1}{2}}),$
$\displaystyle
T_{3}=\tau(D_{\tau}r^{n+1}-r^{n+\frac{1}{2}}_{t}){-}\frac{\tau}{2}\int_{\Omega}H(u^{n+\frac{1}{2}})\big{(}D_{\tau}u^{n+1}-u^{n+\frac{1}{2}}_{t}\big{)}d{\bf
x}.$

Multiplying the first equation with $E_{v}^{n+1}-E_{v}^{n}$, the second
equation with $E_{u}^{n+1}-E_{u}^{n}$, and the third equation with
$E_{r}^{n+1}+E^{n}_{r}$ in (3.15), and then summing up three equalities, we
get

$\displaystyle\begin{split}&\frac{1}{2}\big{(}\|E_{v}^{n+1}\|^{2}-\|E_{v}^{n}\|^{2}+\|\nabla
E_{u}^{n+1}\|^{2}-\|\nabla
E_{u}^{n}\|^{2}+\lambda\|E_{u}^{n+1}\|^{2}-\lambda\|E_{u}^{n}\|^{2}\big{)}+|E_{r}^{n+1}|^{2}-|E_{r}^{n}|^{2}\\\
&=(I_{2},DE_{u}^{n+1})+(T_{2},DE_{u}^{n+1})-(T_{1},DE_{v}^{n+1})+2(T_{3}+I_{3})\hat{E}_{r}^{n+\frac{1}{2}},\end{split}$
(3.16)

where

$\displaystyle
I_{2}={\hat{r}^{n+\frac{1}{2}}}(H(\tilde{U}^{n+\frac{1}{2}})-H(u^{n+\frac{1}{2}})),$
(3.17) $\displaystyle
I_{3}=\frac{1}{2}\int_{\Omega}\big{(}H(u^{n+\frac{1}{2}})-H(\tilde{U}^{n+\frac{1}{2}})\big{)}({u^{n+1}-u^{n}})d{\bf
x}.$ (3.18)

We next estimate the terms $T_{i},I_{i},i\leq 3$ respectively. By Taylor
expansion, there holds

$\displaystyle\|T_{1}\|\lesssim\tau^{2},\ \ \|T_{2}\|\lesssim\tau^{2},\ \
|T_{3}|\lesssim\tau^{3}.$

By (3.14) and the fact that $f\in C^{2}(\mathbb{R})$, we have

$|H(U^{n})|+|H^{\prime}(U^{n})|+|f^{\prime}(U^{n})|+|f^{\prime\prime}(U^{n})|\lesssim
1,\ \ \forall n.$

Then there exists some $\theta\in(0,1)$ such that

$\displaystyle\|H(\tilde{U}^{n+\frac{1}{2}})-H(u^{n+\frac{1}{2}})\|$
$\displaystyle=$
$\displaystyle\|H^{\prime}(\theta\tilde{U}^{n+\frac{1}{2}}+(1-\theta)u^{n+\frac{1}{2}})(\tilde{E}_{u}^{n+\frac{1}{2}}+\tilde{u}^{n+\frac{1}{2}}-u^{n+\frac{1}{2}})\|$
$\displaystyle\lesssim$
$\displaystyle\tau^{2}+\|\tilde{E}_{u}^{n+\frac{1}{2}}\|,$

and thus

$\displaystyle\|I_{2}\|\lesssim\tau^{2}+\|\tilde{E}_{u}^{n+\frac{1}{2}}\|,\ \
|I_{3}|\lesssim\tau(\tau^{2}+\|\tilde{E}_{u}^{n+\frac{1}{2}}\|).$

Consequently,

$|(T_{3}+I_{3})\hat{E}_{r}^{n+\frac{1}{2}}|\lesssim\tau|\hat{E}_{r}^{n+\frac{1}{2}}|^{2}+\tau^{-1}(|I_{3}|^{2}+|T_{3}|^{2})\lesssim\tau^{5}+\tau(|\hat{E}_{r}^{n+\frac{1}{2}}|^{2}+\|\tilde{E}_{u}^{n+\frac{1}{2}}\|^{2}).$
(3.19)

Note that

$\|DE_{u}^{n+1}\|=\tau\|\hat{E}_{v}^{n+\frac{1}{2}}+T_{1}\|\lesssim\tau^{3}+\tau\|\hat{E}_{v}^{n+\frac{1}{2}}\|.$

Then

$|(I_{2}+T_{2},DE_{u}^{n+1})|\lesssim\tau\|I_{2}\|^{2}+\tau\|T_{2}\|^{2}+\tau^{-1}\|DE_{u}^{n+1}\|^{2}\lesssim\tau^{5}+\tau(\|\hat{E}_{v}^{n+\frac{1}{2}}\|^{2}+\|\tilde{E}_{u}^{n+\frac{1}{2}}\|{{}^{2}}).$
(3.20)

On the other hand, in light of the second equation of (3.15), we have

$\displaystyle(T_{1},DE_{v}^{n+1})$ $\displaystyle=$
$\displaystyle\tau(T_{1},\triangle\hat{E}_{u}^{n+\frac{1}{2}}{-}\lambda\hat{E}_{u}^{n+\frac{1}{2}}-\hat{r}^{n+\frac{1}{2}}H(u^{n+\frac{1}{2}})+\hat{R}^{n+\frac{1}{2}}H(\tilde{U}^{n+\frac{1}{2}})+T_{2})$
$\displaystyle=$ $\displaystyle\\!-\\!\tau(\nabla
T_{1},\\!\nabla\hat{E}_{u}^{n+\frac{1}{2}})\\!+\\!\tau(T_{1},{-}\lambda\hat{E}_{u}^{n\\!+\\!\frac{1}{2}}\\!-\\!\hat{r}^{n\\!+\\!\frac{1}{2}}H(u^{n\\!+\\!\frac{1}{2}})\\!+\\!\hat{R}^{n\\!+\\!\frac{1}{2}}H(\tilde{U}^{n\\!+\\!\frac{1}{2}})\\!+\\!T_{2}).$

Noticing that

$\displaystyle\|\hat{r}^{n+\frac{1}{2}}H(u^{n+\frac{1}{2}}){-}\hat{R}^{n+\frac{1}{2}}H(\tilde{U}^{n+\frac{1}{2}})\|$
$\displaystyle=$
$\displaystyle\|\hat{E}_{r}^{n+\frac{1}{2}}H(u^{n+\frac{1}{2}})+\hat{R}^{n+\frac{1}{2}}(H(u^{n+\frac{1}{2}})-H(\tilde{U}^{n+\frac{1}{2}}))\|$
$\displaystyle\lesssim$
$\displaystyle|E_{r}^{{n}+\frac{1}{2}}|+\tau^{2}+\|\tilde{E}_{u}^{n+\frac{1}{2}}\|,$

we have

$|(T_{1},DE_{v}^{n+1})|\lesssim\tau^{3}(\|\nabla\hat{E}_{u}^{n+\frac{1}{2}}\|+\|\hat{E}_{u}^{n+\frac{1}{2}}\|^{2}+\tau^{2}+|E_{r}^{{n}+\frac{1}{2}}|+\|\tilde{E}_{u}^{n+\frac{1}{2}}\|).$
(3.21)

Plugging (3.19)-(3.21) into (3.16) yields

$\displaystyle\frac{1}{2}\big{(}\|E_{v}^{n\\!+\\!1}\|^{2}\\!-\\!\|E_{v}^{n}\|^{2}\\!+\\!\|\nabla
E_{u}^{n\\!+\\!1}\|^{2}\\!-\\!\|\nabla
E_{u}^{n}\|^{2}\\!+\\!\lambda\|E_{u}^{n+1}\|^{2}\\!-\\!\lambda\|E_{u}^{n}\|^{2}\big{)}\\!+\\!(E_{r}^{n+1})^{2}\\!-\\!(E_{r}^{n})^{2}$
$\displaystyle\lesssim\tau^{5}+\tau(\|\tilde{E}_{u}^{n+\frac{1}{2}}\|^{2}+\|\hat{E}_{v}^{n+\frac{1}{2}}\|^{2}+\|\nabla\hat{E}_{u}^{n+\frac{1}{2}}\|^{2}+\|\hat{E}_{u}^{n+\frac{1}{2}}\|^{2}+|\hat{E}_{r}^{n+\frac{1}{2}}|^{2}\big{)}.$

Summing up all $n$ from $0$ to $m$ and using the initial values, we get

$\displaystyle\|E_{v}^{m+1}\|^{2}+\|\nabla
E_{u}^{m+1}\|^{2}+\lambda\|E_{u}^{m+1}\|^{2}+|E_{r}^{m+1}|^{2}$
$\displaystyle\leq C\tau^{4}+C\tau\sum_{n=0}^{m+1}(\|E_{v}^{n}\|^{2}+\|\nabla
E_{u}^{n}\|^{2}+\lambda\|E_{u}^{n}\|^{2}+|E_{r}^{n}|^{2}).$

Then the desired result follows from the conclusion in Lemma 1.

##  4 Error estimates for the fully discrete solution

In this section, we establish the error estimates for the solution of (2.2).
By the error decomposition, we have

$w^{n}-w_{h}^{n}=w^{n}-W^{n}+W^{n}-w_{h}^{n}=w^{n}-W^{n}+e_{w}^{n},\ \
w=u,v,r.$

In light of the conclusion in Theorem 1, we only need to estimate the term
$e_{w},w=u,v,r$. To this end, we first define the Ritz projection operator
$R_{h}:H_{0}^{1}(\Omega)\rightarrow V_{h}$ by

$\displaystyle\left(\nabla(v-R_{h}v),\nabla\omega\right)=0,\quad\forall\omega\in
V_{h}.$

Then $e^{n}_{u}$ (similar for $e_{v}^{n}$) can be decomposed into

$e_{u}=U^{n}-u^{n}_{h}=\xi^{n}_{u}+\eta^{n}_{u},\ \
\xi^{n}_{u}=R_{h}U^{n}-u^{n}_{h},\ \ \eta_{u}=U^{n}-R_{h}U^{n}.$

According to the standard FEM theory [32], it holds that

$\displaystyle\|v-R_{h}v\|_{L^{2}}+h\|\nabla(v-R_{h}v)\|_{L^{2}}\leq
Ch^{s}\|v\|_{H^{s}},\quad\forall v\in H^{s}(\Omega)\cap H_{0}^{1}(\Omega)$
(4.1)

for $1\leq s\leq k+1$.

Note that the exact solutions of (3.1) satisfy

$\displaystyle\begin{split}&\big{(}D_{\tau}U^{n+1},w_{h}\big{)}=\big{(}\hat{V}^{n+\frac{1}{2}},w_{h}\big{)},&&\\\
&\big{(}D_{\tau}V^{n+1},\zeta_{h}\big{)}={-}\big{(}\nabla\hat{U}^{n+\frac{1}{2}},\nabla\zeta_{h}\big{)}{-}\lambda(\hat{U}^{n+\frac{1}{2}},\zeta_{h})-\big{(}{\hat{R}^{n+\frac{1}{2}}}H(\tilde{U}^{n+\frac{1}{2}}),\zeta_{h}\big{)},&&\\\
&{R^{n+1}-R^{n}}=\frac{1}{2}\int_{\Omega}H(\tilde{U}^{n+\frac{1}{2}})({U^{n+1}-U^{n}})d{\bf
x},&&\end{split}$ (4.2)

Subtracting (2.2) from (4.2) gives the following error equation

$\displaystyle\begin{split}&\big{(}D_{\tau}\xi_{u}^{n+1},w_{h}\big{)}=\big{(}\hat{\xi}_{v}^{n+\frac{1}{2}},w_{h}\big{)}+R_{1}(w_{h}),&&\\\
&\big{(}D_{\tau}\xi_{v}^{n+1},\zeta_{h}\big{)}=-\big{(}\nabla\hat{\xi}_{u}^{n+\frac{1}{2}},\nabla\zeta_{h}\big{)}{-}\lambda\big{(}\hat{\xi}_{u}^{n+\frac{1}{2}},\zeta_{h}\big{)}-\big{(}{\hat{e}_{r}^{n+\frac{1}{2}}}H(\tilde{u}_{h}^{n+\frac{1}{2}}),\zeta_{h}\big{)}+R_{2}(\zeta_{h})+I_{2}(\zeta_{h}),&&\\\
&{e_{r}^{n+1}-e_{r}^{n}}=\frac{1}{2}\int_{\Omega}H(\tilde{u}_{h}^{n+\frac{1}{2}})({\xi_{u}^{n+1}-\xi_{u}^{n}})d{\bf
x}+R_{3}+I_{3},&&\end{split}$ (4.3)

where $H(u)=\frac{f(u)}{\sqrt{E(u)}}$, and

$\displaystyle\begin{split}&&R_{1}(w_{h})=\big{(}\hat{\eta}_{v}^{n+\frac{1}{2}},w_{h}\big{)}-\big{(}D_{\tau}\eta_{u}^{n+1},w_{h}\big{)},\\\
&&R_{2}(\zeta_{h})=-\big{(}\nabla\hat{\eta}_{u}^{n+\frac{1}{2}},\nabla\zeta_{h}\big{)}{-}\lambda\big{(}\hat{\eta}_{u}^{n+\frac{1}{2}},\zeta_{h}\big{)}-\big{(}D_{\tau}\eta_{v}^{n+1},\zeta_{h}\big{)},\\\
&&I_{2}(\zeta_{h})=\big{(}{\hat{R}^{n+\frac{1}{2}}}(H(\tilde{u}_{h}^{n+\frac{1}{2}})-H(\tilde{U}^{n+\frac{1}{2}})),\zeta_{h}\big{)},\\\
&&R_{3}=\frac{1}{2}\int_{\Omega}H(\tilde{u}_{h}^{n+\frac{1}{2}})({\eta_{u}^{n+1}-\eta_{u}^{n}})d{\bf
x},\\\
&&I_{3}=\frac{1}{2}\int_{\Omega}\big{(}H(\tilde{U}^{n+\frac{1}{2}})-H(\tilde{u}_{h}^{n+\frac{1}{2}})\big{)}({{U}^{n+1}-U^{n}})d{\bf
x}.\end{split}$ (4.4)

###### Theorem 2

Suppose $u$ is the solution of (1.1) satisfying

$\|u_{0}\|_{H^{k+1}}+\|u\|_{L^{\infty}((0,T),H^{k+1})}+\|u_{t}\|_{L^{2}((0,T),H^{k+1})}+\|u_{tt}\|_{L^{2}((0,T),H^{k+1})}\lesssim
1,$

and $(u_{h}^{n},v_{h}^{n},r_{h}^{n})$ is the solution of (2.2) with
$(u_{h}^{0},v_{h}^{0},r_{h}^{0})=(R_{h}u_{0},R_{h}u_{1},\sqrt{E(u_{0})})$.
Then

$\displaystyle\|R_{h}v^{n}-v_{h}^{n}\|+\|R_{h}u^{n}-u_{h}^{n}\|_{H^{1}}+|r^{n}-r_{h}^{n}|\lesssim
h^{k+1}+\tau^{2}.$ (4.5)

Proof. We first estimate the terms $\xi_{u}^{n},\xi_{v}^{n}$. By taking
$(w_{h},\zeta_{h})=(D\xi_{v}^{n+1},D\xi_{u}^{n+1})$ and multiplying the third
equation of (4.3) by $2\hat{e}_{r}^{n+\frac{1}{2}}={e_{r}^{n+1}+e_{r}^{n}}$,
we derive

$\displaystyle\frac{1}{2}\big{(}\|\xi_{v}^{n+1}\|^{2}-\|\xi_{v}^{n}\|^{2}+\|\nabla\xi_{u}^{n+1}\|^{2}-\|\nabla\xi_{u}^{n}\|^{2}+\lambda\|\xi_{u}^{n+1}\|^{2}-\|\lambda\xi_{u}^{n}\|^{2}\big{)}+(e_{r}^{n+1})^{2}-(e_{r}^{n})^{2}$
$\displaystyle=$ $\displaystyle
R_{2}(D\xi_{u}^{n+1})+I_{2}(D\xi_{u}^{n+1})-R_{1}(D\xi_{v}^{n+1})+2(R_{3}+I_{3})\hat{e}_{r}^{n+\frac{1}{2}},$

where $R_{i},I_{i},i\leq 3$ are given in (4.4). Summing up all $n$ from $0$ to
$m$ and using the initial error $\xi_{u}^{0}=\xi_{v}^{0}=0$, we get

$\displaystyle\begin{split}&\frac{1}{2}(\|\xi_{v}^{{m}+1}\|^{2}+\|\nabla\xi_{u}^{m+1}\|^{2}+\lambda\|\xi_{u}^{m+1}\|^{2})+|e_{r}^{m+1}|^{2}&\\\
&=\sum_{n=0}^{m}\left(R_{2}(D\xi_{u}^{n+1})+I_{2}(D\xi_{u}^{n+1})-R_{1}(D\xi_{v}^{n+1})+2(R_{3}+I_{3})\hat{e}_{r}^{n+\frac{1}{2}}\right).&\end{split}$
(4.6)

To estimate the terms in the right hand side of (4.6), we shall first make the
hypothesis that there exists a positive constant $C_{*}$ such that

$\|u_{h}^{n}\|_{L^{\infty}}\leq C_{*}.$ (4.7)

This hypothesis will be verified later by using the method of mathematical
induction.

Due to (4.7) and the fact that $f\in C^{2}(\mathbb{R})$, we have

$|H(u_{h}^{n})|+|H^{\prime}(u_{h}^{n})|+|f^{\prime}(u_{h}^{n})|+|f^{\prime\prime}(u_{h}^{n})|\lesssim
1,\ \ \forall n.$

Then

$|R_{3}|\lesssim\tau h^{k+1},\ \ |R_{2}(D\xi_{u}^{n+1})|\lesssim
h^{k+1}\|D\xi_{u}^{n+1}\|.$

By Taylor expansion, there exists a $\theta\in(0,1)$ such that

$\displaystyle\|H(\tilde{u}_{h}^{n+\frac{1}{2}})-H(\tilde{U}^{n+\frac{1}{2}})\|$
$\displaystyle=$
$\displaystyle\|H^{\prime}(\theta\tilde{u}_{h}^{n+\frac{1}{2}}+(1-\theta)\tilde{U}^{n+\frac{1}{2}})(\tilde{u}_{h}^{n+\frac{1}{2}}-\tilde{U}^{n+\frac{1}{2}})\|$
(4.8) $\displaystyle\lesssim$ $\displaystyle
h^{k+1}+\|\tilde{\xi}_{u}^{n+\frac{1}{2}}\|.$

Then

$\displaystyle|I_{2}(D\xi_{u}^{n+1})|\lesssim(h^{k+1}+\|\tilde{\xi}_{u}^{n+\frac{1}{2}}\|)\|D\xi_{u}^{n+1}\|,\
\ |I_{3}|\lesssim\tau(h^{k+1}+\|\tilde{\xi}_{u}^{n+\frac{1}{2}}\|),$

and thus

$\big{|}\sum_{n=0}^{m}2(R_{3}+I_{3})\hat{e}_{r}^{n+\frac{1}{2}}\big{|}\leq
Ch^{2(k+1)}+\tau\sum_{n=1}^{m}\big{(}\|\tilde{\xi}_{u}^{n+\frac{1}{2}}\|^{2}+|\hat{e}_{r}^{n+\frac{1}{2}}|^{2}\big{)}.$
(4.9)

On the other hand, we choose $w_{h}=D\xi_{u}^{n+1}$ in (4.3) to get

$\|D\xi_{u}^{n+1}\|\lesssim\tau(\|\hat{\xi}_{v}^{n+\frac{1}{2}}\|+h^{k+1}),$

which yields, together with Cauchy-Schwarz inequality

$|I_{2}(D\xi_{u}^{n+1})|+|R_{2}(D\xi_{u}^{n+1})|\lesssim\tau
h^{2(k+1)}+\tau(\|\hat{\xi}_{v}^{n+\frac{1}{2}}\|^{2}+\|\tilde{\xi}_{u}^{n+\frac{1}{2}}\|^{2}).$

Consequently,

$\sum_{n=0}^{m}|I_{2}(D\xi_{u}^{n+1})|+|R_{2}(D\xi_{u}^{n+1})|\lesssim
h^{2(k+1)}+\tau\sum_{n=0}^{m}(\|\hat{\xi}_{v}^{n+\frac{1}{2}}\|^{2}+\|\tilde{\xi}_{u}^{n+\frac{1}{2}}\|^{2}).$
(4.10)

As for the term $\sum_{n=0}^{m}R_{1}(D\xi_{v}^{n+1})$ in (4.6), we recall the
definition of $R_{1}$ in (4.4) to get

$\displaystyle\big{|}\sum_{n=0}^{m}R_{1}(D\xi_{v}^{n+1})\big{|}$
$\displaystyle=$
$\displaystyle\big{|}\big{(}\hat{\eta}_{v}^{m+\frac{1}{2}}-D_{\tau}\eta_{u}^{m+1},\xi_{v}^{m+1}\big{)}+\sum_{n=1}^{m}\big{(}\hat{\eta}_{v}^{n-\frac{1}{2}}-\hat{\eta}_{v}^{n+\frac{1}{2}}+D_{\tau}\eta_{u}^{n+1}-D_{\tau}\eta_{u}^{n},\xi_{v}^{n}\big{)}\big{|}$
$\displaystyle\leq$ $\displaystyle Ch^{k+1}\|\xi_{v}^{m+1}\|+C\tau
h^{k+1}\sum_{n=1}^{m}\|\xi_{v}^{n}\|$ $\displaystyle\leq$ $\displaystyle
Ch^{2(k+1)}+C\tau\sum_{n=1}^{m}\|\xi_{v}^{n}\|^{2}+\frac{1}{4}\|\xi_{v}^{m+1}\|^{2}.$

Substituting the above inequality, (4.9)-(4.10) into (4.6), we get

$\|\xi_{v}^{m+1}\|^{2}+\|\nabla\xi_{u}^{m+1}\|^{2}+\|\xi_{u}^{m+1}\|^{2}+|e_{r}^{m+1}|^{2}\lesssim
h^{2(k+1)}+\tau\sum_{n=0}^{m+1}\big{(}\|\xi_{v}^{n}\|^{2}+|e_{r}^{n}|^{2}+\|\xi_{u}^{n}\|^{2}\big{)}.$

By the Gronwall inequality given in Lemma 1,

$\|\xi_{v}^{n}\|+\|\xi_{u}^{n}\|_{H^{1}}+|e_{r}^{n}|\lesssim h^{k+1},\ \
\forall n\geq 1.$ (4.11)

Then from the triangle inequality and the conclusion in Theorem 1,

$\displaystyle\|R_{h}v^{n}-v_{h}^{n}\|+\|R_{h}u^{n}-u_{h}^{n}\|_{H^{1}}+|r^{n}-r_{h}^{n}|$
$\displaystyle\leq\|\xi_{v}^{n}\|+\|\xi_{u}^{n}\|_{H^{1}}+|e_{r}^{n}|+\|R_{h}(v^{n}-V^{n})\|+\|R_{h}(u^{n}-U^{n})\|+|r^{n}-R^{n}|\lesssim
h^{k+1}+\tau^{2}.$

The proof is complete.

###### Remark 2

Note that the optimal convergence rate for the $H^{1}$ error approximation is
$O(h^{k})$. The error estimate in (4.5) indicates that the Galerkin SAV
solution $u_{h}^{n}$ is superclose to the Ritz projection of the exact
solution $R_{h}u^{n}$ under the $H^{1}$ norm, which is one order higher than
the counterpart optimal converge rate. As a direct consequence of (4.5), we
have the following optimal error estimates:

$\displaystyle\|v^{n}-v_{h}^{n}\|_{H^{1}}+\|u^{n}-u_{h}^{n}\|_{H^{1}}\lesssim
h^{k}+\tau^{2},$
$\displaystyle\|v^{n}-v_{h}^{n}\|+\|u^{n}-u_{h}^{n}\|+|r^{n}-r_{h}^{n}|\lesssim
h^{k+1}+\tau^{2}.$

To end with this section, we prove the inequality (4.7).

###### Lemma 3

Under the conditions of Theorem 2, it holds that

$\|u_{h}^{n}\|_{L^{\infty}}\leq C_{*},\ \ \forall n\geq 1,$

where the constant $C_{*}$ is independent of $\tau$ and $h$.

Proof. We will show the above inequality by induction. To this end, we first
denote by ${\cal I}_{h}U\in V_{h}$ the interpolation function of $U$. By the
approximation theory, we have

$\|{\cal I}_{h}U-U\|+\|{\cal I}_{h}U-R_{h}U\|\lesssim h^{2}\|U\|_{H^{2}},\ \
\|{\cal I}_{h}U\|_{L^{\infty}}\lesssim\|U\|_{L^{\infty}}.$

Note that

$\|u_{h}^{0}\|_{L^{\infty}}=\|R_{h}u^{0}\|_{L^{\infty}}\leq C.$

By (4.11) and the inverse inequality $\|v_{h}\|_{L^{\infty}}\lesssim
h^{-\frac{d}{2}}\|v_{h}\|$ for all $v_{h}\in V_{h}$, we get

$\displaystyle\|u_{h}^{1}\|_{L^{\infty}}$ $\displaystyle\leq$
$\displaystyle\|u_{h}^{1}-R_{h}U^{1}\|_{L^{\infty}}+\|R_{h}U^{1}-{\cal
I}_{h}U^{1}\|_{L^{\infty}}+\|{\cal I}_{h}U^{1}\|_{L^{\infty}}$
$\displaystyle\leq$ $\displaystyle
C(h^{2-\frac{d}{2}}+\|U^{1}\|_{L^{\infty}})\leq
C_{1}(h^{2-\frac{d}{2}}+\|U^{1}\|_{H^{2}}).$

Now we choose a positive constant $h_{1}$ which is small enough to satisfy

$C_{1}h_{1}^{\frac{1}{2}}\leq C.$

Then for $h\in(0,h_{1}]$, we derive that

$\|u_{h}^{1}\|_{L^{\infty}}\leq C+\|U^{1}\|_{H^{2}}\leq C_{3}.$

Therefore, we can choose the positive constant $C_{*}$ independent of $h$ and
$\tau$ such that

$C_{*}\geq\max\\{2\|U^{n}\|_{H^{2}},\|u_{h}^{1}\|_{L^{\infty}}\\}.$

Then (4.7) is valid for $n=1$. Next, suppose (4.7) holds for all $l\leq n-1$.
We will show that it is also valid for $n$. Thanks to (4.11), we have

$\|u_{h}^{n}-R_{h}U^{n}\|_{H^{1}}\leq Ch^{k+1}.$

Then

$\displaystyle\|u_{h}^{n}\|_{L^{\infty}}$ $\displaystyle\leq$
$\displaystyle\|u_{h}^{n}-R_{h}U^{n}\|_{L^{\infty}}+\|R_{h}U^{n}-{\cal
I}_{h}U^{n}\|_{L^{\infty}}+\|{\cal I}_{h}U^{n}\|_{L^{\infty}}$
$\displaystyle\leq$ $\displaystyle
Ch^{-\frac{d}{2}}(\|u_{h}^{n}-R_{h}U^{n}\|+\|R_{h}U^{n}-{\cal
I}_{h}U^{n}\|)+\|U^{n}\|_{L^{\infty}}\leq
C_{1}h^{\frac{1}{2}}+\frac{C_{*}}{2}.$

Let $h_{1}$ be small enough to satisfy

$C_{1}h_{1}^{\frac{1}{2}}\leq\frac{C_{*}}{2}.$

Then for $h\in(0,h_{1}]$, we derive that

$\|u_{h}^{n}\|_{L^{\infty}}\leq C_{1}h^{\frac{1}{2}}+\frac{C_{*}}{2}\leq
C_{*}.$

This completes the induction.

##  5 Numerical simulations

We present several numerical results to confirm our theoretical findings in
this section.

Example 1 Consider the following Klein-Gordon equation

$\displaystyle
u_{tt}=u_{xx}+u_{yy}+u-u^{3}+g_{1}(x,y,t).~{}~{}~{}~{}~{}~{}(x,y,t)\in[0,1]^{2}\times[0,T],$
(5.1)

where $u(x,y,0)$, $u_{t}(x,y,0)$ and $g_{1}(x,y,t)$ are given by the exact
solution

$\displaystyle u(x,t)=\exp(-t)x^{2}(1-x)^{2}y^{2}(1-y)^{2}.$ (5.2)

We test convergence orders of the fully discrete scheme using uniform
triangulation with $M+1$ nodes in each spatial direction, and take $N=M$ and
$N=M^{\frac{3}{2}}$ for the linear finite element method (L-FEM) and quadratic
finite element method (Q-FEM), respectively. We list the errors at time $T=1$
as well as the convergence rates in Table 1. Here and below, we denote

$\|e\|_{0}=\|u^{N}-u_{h}^{N}\|,~{}~{}~{}\|e\|_{1}=\|R_{h}u^{N}-u_{h}^{N}\|_{H^{1}}.$

These results indicate that the fully discrete scheme is convergent and has
the order $\mathcal{O}(\tau^{2}+h^{r+1})$. We also test the unconditional
convergence of the fully discrete scheme with different spatial step-sizes for
every fixed $\tau$. The $l_{2}$-errors at time $T=1$ are shown in Figure 1.
When the temporal stepsize is fixed, the $L_{2}$ errors tend to a constant.
They imply that the error estimates hold without any temporal mesh sizes
restrictions dependent on the spatial mesh sizes.

Figure 1: $L^{2}$-errors of linear and quadratic finite element approximation Table 1:  Errors and convergent orders for 2D problems. |  | L-FEM |  |  | Q-FEM |
---|---|---|---|---|---|---
M | $\|e\|_{0}$ | order | $\|e\|_{1}$ | order | $\|e\|_{0}$ | order | $\|e\|_{1}$ | order
8 | 4.54E-4 | – | 5.56e-4 | – | 2.47E-6 | – | 4.60E-5 | –
16 | 1.11E-4 | 2.03 | 1.61E-4 | 1.83 | 2.71E-7 | 3.19 | 6.85E-6 | 2.74
24 | 4.90e-5 | 2.02 | 7.40E-5 | 1.96 | 7.66E-8 | 3.12 | 2.41E-6 | 2.57
32 | 2.74e-5 | 2.02 | 4.24E-5 | 1.96 | 3.18E-8 | 3.05 | 9.81E-7 | 3.10
40 | 1.75E-5 | 2.01 | 2.73E-5 | 1.97 | 1.62E-8 | 3.02 | 4.88E-7 | 3.12
Figure 2: The evolutions of the discrete energies

Then, we set $g_{1}(x,y,t)=0$, $T=100$, $N=10$, $M=10$, and solve the problem
using L-FEM. The time discretization is achieved by the linearized Crank-
Nicolson (LCN) method, and by the proposed SAV LCN method, respectively. The
evolutions of the discrete energies are shown in Figure 2. Clearly, the
energies obtained by the LCN finite element method increase as time goes on,
while the one obtained by our method remains the same. It implies that
numerical solutions by the SAV approach conserve the energy.

Example 2 Consider the following Sine-Gordon equation

$\displaystyle
u_{tt}=u_{xx}+u_{yy}+u_{zz}+\sin(u)+g_{2}(x,y,z,t),~{}~{}~{}~{}~{}~{}(x,y,z,t)\in[0,1]^{3}\times[0,1],$
(5.3)

where the initial conditions and $g_{2}(x,y,z,t)$ are produced from the exact
solution

$\displaystyle u(x,t)=(1+t^{3})\sin(2\pi x)\sin(2\pi y)\sin(2\pi z).$ (5.4)

Table 2:  Errors and convergent orders for 3D problems. |  | L-FEM |  |  |  | Q-FEM |
---|---|---|---|---|---|---|---
M | $\|e\|_{0}$ | order | $\|e\|_{1}$ | order | M | $\|e\|_{0}$ | order | $\|e\|_{1}$ | order
16 | 6.36E-2 | – | 4.43E-1 | – | 10 | 5.69E-3 | – | 1.29E-1 | –
20 | 4.17E-2 | 1.89 | 2.94E-1 | 1.83 | 12 | 3.24E-3 | 3.09 | 7.41E-2 | 2.96
24 | 2.94E-2 | 1.92 | 2.08E-1 | 1.91 | 14 | 1.99E-3 | 3.16 | 4.91E-2 | 2.67
28 | 2.18E-2 | 1.94 | 1.55E-1 | 1.90 | 16 | 1.31E-3 | 3.13 | 3.40E-2 | 2.75
Figure 3: The evolutions of the discrete energies

We still take $N=M$ and $N=M^{\frac{3}{2}}$ for the linear and quadratic
finite element approximations, respectively. The numerical errors at time
$T=1$ as well as the convergence rates are presented in Table 2. The given
results indicate that the fully discrete scheme has the order
$\mathcal{O}(\tau^{2}+h^{r+1})$.

Next, we set $g_{2}(x,y,z,t)=0$, $T=20$, $N=10$, $M=10$ and solve the problem
by the linear finite element method. The evolutions of the discrete energies
for the 3D problems are displayed in Fig. 3. Clearly, the discrete energies by
the SAV approach remain unchanged, while the ones obtained by the LCN finite
element method increase as time goes on. They further confirm the findings in
this study.

##  6 Conclusion

In this study, we present a linearly implicit numerical schemes for solving
the nonlinear wave equation (1.1). The scheme is developed by combining the
SAV approach with finite element methods, classical Crank-Nicolson methods,
and extrapolation approximation. The fully discrete scheme is proved to be
unconditionally convergent and energy-conserving. Numerical illustrations are
presented to confirm the theoretical findings.

## References

  * [1] G. Akrivis, B. Li, and D. Li,  Energy-Decaying Extrapolated RK–SAV Methods for the Allen–Cahn and Cahn–Hilliard Equations,  SIAM. J. Sci. Comput., 41 (2019) A3703-A3727.
  * [2] R. Bellman,  The stability of solutions of linear differential equations,  Duke Math. J., 10 (1943), 643-647.
  * [3] A. Biswas,  Soliton perturbation theory for phi-four model and nonlinear Klein-Gordon equations.,  Commun. Nonlinear Sci. Numer. Simul., 14 (2009) 3239-3249.
  * [4] L. Brugnano, G. F. Caccia, and F. Iavernaro,  Energy conservation issues in the numerical solution of the semilinear wave equation,  Appl. Math. Comput., 270 (2015) 842-870.
  * [5] W. Cai, C. Jiang, Y. Wang, and Y. Song,  Structure-preserving algorithms for the two-dimensional sine-Gordon equation with Neumann boundary conditions,  J. Comput. Phys., 395 (2019) 166-185.
  * [6] J. Cai and J. Shen,  Two classes of linearly implicit local energy-preserving approach for general multi-symplectic Hamiltonian PDEs,  J. Comput. Phys., (2019) https://doi.org/10.1016/j.jcp.2019.108975.
  * [7] J. R. Cannon and Y. Lin,  Nonclassical H1 projection and Galerkin methods for nonlinear parabolic integro-differential equations,  SIAM. J. Numer. Anal., 25 (1988) 187–201.
  * [8] B. Cano and M. J. Moreta,  Multistep cosine methods for second–order partial differential systems,  IMA J. Numer. Anal., 30 (2010) 431–461.
  * [9] W. Cao, D. Li, and Z. Zhang,  Optimal superconvergence of energy conserving local discontinuous Galerkin methods for wave equations,  Commu. Comput. Phys., 21 (2017) 211–236.
  * [10] C. Carstensen and G. Dolzmann,  Time-Space Discretization of the Nonlinear Hyperbolic System $u_{tt}=div(\sigma(\mathbf{D}U)+\mathbf{D}U_{t})$,  SIAM. J. Numer. Anal., 42 (2004) 75–89.
  * [11] J. E. Dendy,  Galerkin’s Method for Some Highly Nonlinear Problems,  SIAM. J. Numer. Anal., 14 (1977) 327-347.
  * [12] R. K. Dodd, I. C. Eilbeck, J. D. Gibbon, and H. C. Morris,  Solitons and Nonlinear Wave Equations.,  Academic, 1982\.
  * [13] P. J. Drazin and R. S. Johnson,  Solitons: An Introduction,  Phy. Tod., 43 (1990) 70-71.
  * [14] I. Farago,  Finite element method for solving nonlinear parabolic equations,  Computers Math. Applic., 21 (1991) 59–69.
  * [15] S. Garcia,  Improved error estimates for mixed finite-element approximations for nonlinear parabolic equations: the discrete-time case,  Numer. Methods Partial Differ. Equ., 10 (1994) 149–169.
  * [16] L. Gauckler,  Error analysis of trigonometric integrators for semilinear wave equations,  SIAM J. Numer. Anal., 53 (2015) 1082–1106.
  * [17] V. Grimm,  A note on the Gautschi-type method for oscillatory second-order differential equations,  Numer. Math., 102 (2005) 61–66.
  * [18] J. G. Heywood and R. Rannacher,  Finite element approximation of the nonstationary Navier-Stokes problem IV: Error analysis for second-order time discretization,  SIAM. J. Numer. Anal., 27 (1990) 353-384.
  * [19] B. Li, H. Gao, and W. Sun,  Unconditional optimal error estimates of a Crank-Nicolson Galerkin method the nonlinear thermistor equations,  SIAM J. Numer. Anal.,  52 (2014) 933–954.
  * [20] B. Li and W. Sun,  Error analysis of linearized semi-implicit Galerkin finite element methods for nonlinear parabolic equations,  Int. J. Numer. Anal. Model., 10 (2013) 622–633.
  * [21] B. Li and W. Sun,  Unconditional convergence and optimal error estimates of a Galerkin-mixed FEM for incompressible miscible flow in porous media,  SIAM J. Numer. Anal., 51 (2013) 1959–1977.
  * [22] D. Li and W. Sun,  Linearly Implicit and High-Order Energy-Conserving Schemes for Nonlinear Wave Equations,  J. Sci. Comput., 83 (2020) 65.
  * [23] D. Li, J. Wang, and J. Zhang,  Unconditionally convergent $L1$-Galerkin FEMs for nonlinear time-fractional Schrödinger equations,  SIAM. J. Sci. Comput., 39 (2017) A3067–A3088.
  * [24] D. Li, C. Wu, and Z. Zhang,  Linearized Galerkin FEMs for Nonlinear Time Fractional Parabolic Problems with Non$-$smooth Solutions in Time Direction,  J. Sci. Comput., 80 (2019) 403–419.
  * [25] X. Li, J. Wen, D. Li, Mass- and energy-conserving difference schemes for nonlinear fractional Schrödinger equations. Appl Math Lett, 2021, 111: 106686.
  * [26] M. Luskin,  A Galerkin method for nonlinear parabolic equations with nonlinear boundary conditions,  SIAM J. Numer. Anal. 16 (1979) 284–299.
  * [27] H. H. Rachford Jr.,  Two-level discrete-time Galerkin approximations for second order nonlinear parabolic partial differential equations,  SIAM J. Numer. Anal., 10 (1973) 1010–1026.
  * [28] J. Shen and J. Xu,  Convergence and Error Analysis for the Scalar Auxiliary Variable (SAV) Schemes to Gradient Flows,  SIAM J. Numer. Anal., 56 (2018) 2895-2912.
  * [29] J. Shen, J. Xu, and J. Yang,  The scalar auxiliary variable (SAV) approach for gradient flows,  J. Comput. Phys., 353 (2018) 407-416.
  * [30] J. Shen, J. Xu, and J. Yang,  A New Class of Efficient and Robust Energy Stable Schemes for Gradient Flows,  SIAM Rev., 61 (2019) 474-506.
  * [31] Roger Temam,  Infinite-dimensional dynamical systems in mechanics and physics, 2nd edition,  Springer, 1997\.
  * [32] V. Thomée,  Galerkin finite element methods for parabolic problems,  Springer-Verlag, 1997\.
  * [33] B. Wang and X. Wu,  The formulation and analysis of energy-preserving schemes for solving high-dimensional nonlinear Klein-Gordon equations,  IMA J. Numer. Anal., 39 (2019) 2016-2044.
  * [34] A. M. Wazwaz,  New travelling wave solutions to the Boussinesq and the Klein-Gordon equations.,  Commun. Nonlinear Sci. Numer. Simul., 13 (2008) 889-901.
  * [35] X. Wu, B. Wang, and W. Shi.  Efficient energy preserving integrators for oscillatory Hamiltonian systems,  J. Comput. Phys., 235 (2013) 587-605.


# Scalable Quantum Ground State Preparation of the Heisenberg Model: A
Variational Quantum Eigensolver Approach

Jinao Wang
Tsinghua International School
Beijing, China
<EMAIL_ADDRESS>
&Rimika Jaiswal
Department of Physics
University of California, Santa Barbara
Santa Barbara

###### Abstract

Quantum systems have historically been formidable to simulate using classical
computational methods, particularly as the system size grows. The Heisenberg
Model, pivotal in understanding magnetic materials, is a quintessential
example where classical simulations face scalability issues. The Variational
Quantum Eigensolver (VQE) algorithm is a system composed of a quantum circuit
as well as a classical optimizer that can efficiently prepare the Heisenberg
ground state by iteratively optimizing the variational parameters. We assess
the efficacy and scalability of VQE by preparing the ground states of
isotropic and anisotropic Heisenberg models. This paper also aims to provide
insights into the precision and time consumption involved in classical and
optimized sampling approaches in the calculation of expectation values. In
preparing the ground state for the Heisenberg models, this paper paves the way
for more efficient quantum algorithms and contributes to the broader field of
condensed matter physics.

_K_ eywords Many-Body Physics $\cdot$ Quantum Computing $\cdot$ State
Preparation $\cdot$ Lagrangian Physics $\cdot$ Statistical Mechanics

##  1 Introduction

Recent advancements in quantum computing have enabled researchers to apply
quantum circuits in analyzing the quantum states of complex many-body systems,
which have been previously formidable to classical computational approaches.
When it comes to quantum state preparations, classical computation methods
face significant scalability issues as the computational cost grows
exponentially in the Hilbert Space [1]. In addition, as stated by John
Preskill, the current level of quantum technology is best characterized by
Noisy Intermediate Scale Quantum (NISQ) computing that has high error rates
and limited system size [2]

In order to address the limits posed by current quantum hardware, a class of
algorithms called Variational Quantum Algorithms was proposed. VQE, a subclass
of VQAs, will be closely examined in this paper. First proposed by Alberto
Peruzzo et.al [3], VQE is a hybrid quantum-classical algorithm that has
emerged as a promising approach for solving problems in quantum chemistry,
material science, and condensed matter physics using near-term quantum
computers. VQE is an application of the variational principle, where a quantum
computer is trained to prepare the ground state of a given system, such as a
molecule or a spin lattice [4].

The main goal of this paper is to use VQE to prepare the ground state of the
Heisenberg model variationally. The Heisenberg model is a fundamental
theoretical framework in the study of magnetic systems and plays a crucial
role in understanding the behavior of quantum spin systems. Developed by
Werner Heisenberg, this model describes the interactions between magnetic
moments or spins in a lattice, taking into account their quantum mechanical
nature [5]. By determining the ground state of the Heisenberg model, which
corresponds to the lowest energy configuration of the system, physicists can
better understand its low-temperature properties and phase transitions, as
well as providing insights into the nature of quantum correlations and
entanglement in the system [6].

Current attempts to use the Variational Quantum Eigensolver (VQE) to calculate
the Heisenberg model’s ground state have demonstrated promise and versatility,
yet encounter significant challenges. Early research by Cervera-Lierta et al.
and Kokail et al. validated VQE’s potential, delivering accurate results for
small systems but highlighting scalability issues due to limited qubits and
quantum device noise [7, 8].

Moreover, in the study of the anti-ferromagnetic Heisenberg model on the
kagome lattice, VQE results aligned with exact diagonalization for an 18-site
system, except near phase transition points where it erroneously converged to
excited states [9]. Moll et al. also applied VQE within quantum chemistry,
again underscoring the difficulty of larger system optimization due to its
complex landscape [10].

In tackling the scalability issues in ground state preparations using VQE,
this paper focuses not only on obtaining the target ground state but also
analyzing the computation efficiency of VQE and studying the effect of
sampling optimization on the VQE runtime. In addition to the widely explored
1D isotropic Heisenberg model for small lattices, this paper also studies
larger systems with more than 20 qubits. The less explored XXZ anisotropic
Heisenberg model is also examined in this paper, demonstrating VQE’s
capability in a variety of models.

##  2 Methods

###  2.1 The Variatinal Quantum Eigensolver

Suppose we want to prepare a target ground state $|\psi_{target}\rangle$. The
VQE takes in a trivial input state $|\psi_{0}\rangle$ and finds the
appropriate parametrized unitary circuit that can act on the input state to
generate the target state. The optimization process is implemented on a hybrid
system of classical optimizer and quantum circuit, where a set of parameters
are given to a quantum circuit, and the expectation values are calculated,
then sent back to the classical optimizer where the optimizer updates a new
set of parameters to the quantum circuit.

Figure 1: Schematic overview of the structure of the Variational Quantum
Eigensolver (VQE).

###  2.2 1D Isotropic Heisenberg Model

The Hamiltonian for the 1D Heisenberg model is defined as

$H_{T}=J\sum_{i=1}^{L-1}\boldsymbol{S}_{i}\cdot\boldsymbol{S}_{i+1}$ (1)

where S represents the spin operator of a particle, and J is the coupling
constant that will be defined as 1 for the isotropic model. As for the ansatz
design, this paper follows the quantum computing simulation protocol devised
by Wen Wei Ho and Timothy H. Hsieh [6]. The final target state is written as a
linear combination of and , which will be defined as
$H_{1}=\sum_{i=1}^{L/2-1}\boldsymbol{S}_{2i}\cdot\boldsymbol{S}_{2i+1}$ , and
$H_{2}=\sum_{i=1}^{L/2}\boldsymbol{S}_{2i-1}\cdot\boldsymbol{S}_{2i}$ . Then
the structure of the ansatz can be represented as alternating layers of
$H_{1}$ and $H_{2}$ with $p\in\mathbb{Z}$ parameters for $\beta$ and $\gamma$
respectively:

$|\psi(\boldsymbol{\gamma},\boldsymbol{\beta})\rangle_{p}=e^{-i\beta_{p}H_{1}}e^{-i\gamma_{p}H_{2}}\cdots
e^{-i\beta_{1}H_{1}}e^{-i\gamma_{1}H_{2}}\left|\psi_{1}\right\rangle$ (2)

Figure2 shows the detailed layout of the ansatz design for the 1D Heisenberg
model. The circuit is first initialized into entangled bell state pairs
$\bigotimes_{i}\frac{1}{\sqrt{2}}(|\uparrow\downarrow\rangle-|\downarrow\uparrow\rangle)_{2i-1,2i}$
. Then, parameterized $H_{1}$ and $H_{2}$ layers alternate to form $p$ levels,
with each level containing two layers of operation gates. This design of
alternating H1 and H2 layers is adopted from the Quantum Approximate
Optimization Algorithm (QAOA), which is originally developed to solve
combinatorial optimization problems [11]. This ansatz design is effective for
certain problems that have a feasible subspace that is smaller than the full
Hilbert space [5].

The expectation value with respect to the Hamiltonian (also the cost
function)is efficiently defined as:

$F_{p}(\gamma,\beta)={}_{p}\left\langle\psi(\gamma,\beta)\left|H_{T}\right|\psi(\gamma,\beta)\right\rangle_{p}$
(3)

Figure 2: Schematic overview of the structure of the ansatz.

###  2.3 1D Anisotropic XXZ Heisenberg Model

The XXZ model is a generalization of the 1D Heisenberg model given in the form
below

$H=\sum_{i=1}^{N}\left(J_{x}S_{i}^{x}S_{i+1}^{x}+J_{y}S_{i}^{y}S_{i+1}^{y}+J_{z}S_{i}^{z}S_{i+1}^{z}\right)$
(4)

that satisfies the condition $J_{x}=J_{u}\neq J_{z}$. In this paper, $J_{x}=1$
sets the energy scale and $J_{z}=\Delta$ is the anisotropy parameter. Although
the XXZ model is already exactly solvable with the Bethe ansatz [12, 13], its
ground state preparation has not been explored using the VQE.

The ansatz of the XXZ model is similar to the one used in the calculation for
the 1D Heisenberg model (Figure 2). But here, the anisotropy parameter
$\Delta$needs to be passed into both the Hamiltonian and the variational
circuit. When $\Delta=0$, the model becomes the XY model, neglecting the spin
interaction in the Z-axis [19]. When $\Delta=1$, the model recovers isotropy
with an $SU(2)$ symmetry [14]. For $\Delta\in[-1,1]$, the model is in the
critical phase, which will be selected as the value range in the subsequent
experiments [15].

###  2.4 Sampling Optimization

In addition to the VQE, sampling is also applied in the calculation of the
expectation value to further reduce the computational cost as $N$ increases
[16]. Sampling utilizes the probabilistic nature of quantum mechanics by
randomly taking a subset of all the possible outcomes, instead of calculating
all the possible states. In doing so, the task of calculating the expectation
value is reduced from exponential complexity down to polynomial complexity,
with an error of $\frac{1}{\sqrt{M}}$ given by the central limit theorem [17],
where $M$ denotes the number of samples [18, 19].

##  3 Results and Discussions

The ground state of the aforementioned models is prepared for different system
sizes. To demonstrate the scalability of VQE with sampling optimization,
ground states computations are also conducted on larger lattices and analyzed
by recording the computation time.

###  3.1 1D Isotropic Heisenberg Model

Figure 3: Optimization record of 1D Heisenberg Model for N =
4,6,8,10,12,14,16, with sampling optimization applied.

In Figure 3, the optimization was performed on 1D isotropic Heisenberg model
with varying system sizes. In most cases, the calculated ground state
converges to the correct ground state, but in some cases, such as N = 12 and N
= 6, the optimization energy does not converge exactly to the target energy.
This is likely due to the fact that the use of sampling for the calculation of
the ground state led to more randomized results, while in FIG 9 of Appendix
1A, the prepared ground state through exact calculation of the expectation
value demonstrates higher accuracy.

###  3.2 1D Anisotropic Heisenberg Model

Figure 4: Ground states of the 1D Anisotropic XXZ Heisenberg model with
$\Delta\in[-1,1]$ sampled at an interval length of 0.1, sampling optimization
applied.

For the anisotropic XXZ model, the ground state is prepared for a set of delta
values belonging to the critical phase of the XXZ model exhibiting the
characteristics discussed in section 2.2. The sampling results are then
compared with the calculated values obtained through exact diagonalization
shown in figure 4. The $R^{2}$ of both runs calculated with respect to the
computed energy is above 0.8, indicating a strong correlation between the
optimized value and the ground truth, thereby demonstrating the accuracy of
VQE in preparing the ground state of the XXZ anisotropic models.

###  3.3 Running on Larger Systems and Runtime analysis

Figure 5: Optimization results of 1D Isotropic Heisenberg Model for N =
20,24,28, with sampling optimization applied.

For larger systems, the VQE still converges to the expected ground state for
the three runs conducted, demonstrating its stability and robustness even when
the system size increases.

Figure 6: VQE optimization time comparison between exact computation and
sampling optimization.

When comparing the computation time with and without sampling optimization for
the calculation of the expectation value, classical computation is faster on
smaller-scale systems. There is an overlap between the two approaches at N=14.
But as N increases beyond 14, the runtime of the classically computed cost
function increases significantly while the runtime for sampling increases much
slowly (only polynomially)This shows the exponential growth in the computation
complexity of expectation values classically and also proves the efficiency of
sampling optimization that only has polynomial time complexity.

###  3.4 Result Verification

To verify the result of the calculated ground state, the subsystem entropy, as
well as the correlation function of one set of optimized parameters, is
calculated. If the calculated results demonstrate the expected physical
properties, then it can be confirmed that the results obtained through this
research are valid. Due to the length constraint, only the result for
$\Delta=1$ is presented. For the full result, please refer to Appendix 2a and
2b.

Figure 7: Result verification of the prepared ground state for $N=10$ and
$\Delta=1$. 7(a) the correlation of the prepared state with respect to
r(distance between qubits). 7(b) the subsystem entropy of the 1D Heisenberg.

Graph 7(b) shows an oscillating pattern, confirming our system’s
initialization into a dimer state with entangled Bell pairs. High entropy
occurs when subsystem cuts intersect entangled qubits, while low entropy
arises from unentangled pair intersections, confirming the dimerized
structure. Graph 7(a) exhibits an oscillatory correlation function that starts
strong at short distances, indicative of singlets, then rapidly drops to zero,
typical of a strongly dimerized state for the Heisenberg model. Thus, by
analyzing the graph’s characteristics, it can be concluded that the prepared
state is valid.

##  4 Conclusion

In this paper, the optimized parameters are obtained for the variational
quantum circuit at different system sizes and both the isotropic and
anisotropic models. This is valuable for future researchers who wish to obtain
the ground state and its properties directly through variational quantum
circuits. This can be directly useful in the simulations of magnetic systems.

In addition, this paper also demonstrates the capabilities of the VQE for the
XXZ model, which has not been implemented using VQE in the past. Based on the
obtained optimization results, it can be concluded that VQE can reliably be
used to prepare the ground state of the XXZ anisotropic Heisenberg model.

Sampling optimization is verified to improve the time runtime from exponential
complexity down to polynomial complexity while retaining most of the accuracy.

Lastly, VQE also demonstrates robustness as it can still reliably converge to
the desired ground state even on large systems that are out of reach for
classical computations.

## Acknowledgments

I would like to express my heartfelt gratitude to my research partner, Lana
Wong, for her support and encouragement throughout this research. A big thank
you to my TA Mr. Vamshi Chowdary Madala for giving me exhaustive and detailed
advice on my research project, and answering my questions and concerns with
patience.

Additionally, my sincere appreciation goes to the Program Director of the
Research Mentorship Program (RMP) at UCSB, Dr. Lina Kim, for her wonderful
lectures and support in crafting and presenting all crucial components of my
research project.

Finally, I would like to thank all the UCSB faculty and those who helped
organize RMP 2023. This paper would not have been possible without your
collective support and mentorship.

## References

  * [1] G. Carleo, Y. Nomura, and M. Imada, “Constructing exact representations of quantum many-body systems with deep neural networks,” Nature Communications, vol. 9, no. 1, 2018.
  * [2] J. Preskill, “Quantum computing in the nisq era and beyond,” Quantum, vol. 2, p. 79, 2018.
  * [3] A. Peruzzo, J. Mcclean, P. Shadbolt, M.-H. Yung, X.-Q. Zhou, P. J. Love, A. Aspuru-Guzik, and J. L. O’Brien, “A variational eigenvalue solver on a photonic quantum processor,” Nature Communications, vol. 5, no. 1, 2014\.
  * [4] J. Tilly, H. Chen, S. Cao, D. Picozzi, K. Setia, Y. Li, E. Grant, L. Wossnig, I. Rungger, G. H. Booth, and J. Tennyson, “The variational quantum eigensolver: A review of methods and best practices,” Physics Reports, vol. 986, pp. 1–128, nov 2022.
  * [5] M. Cerezo, A. Arrasmith, R. Babbush, S. C. Benjamin, S. Endo, K. Fujii, J. R. McClean, K. Mitarai, X. Yuan, L. Cincio, and et al., “Variational quantum algorithms,” Nature Reviews Physics, vol. 3, no. 9, p. 625–644, 2021\.
  * [6] A. Hunt, “A Look at VQE: Ground States Help Us Simulate — austinjhunt.medium.com.” https://austinjhunt.medium.com/ground-states-help-to-simulanalyzing-and-implementing-variational-quantum-eigensolver-vqe-with-62ec0331b741.  [Accessed 22-08-2023].
  * [7] A. N. Ciavarella, S. Caspar, M. Illa, and M. J. Savage, “State preparation in the heisenberg model through adiabatic spiraling,” Quantum, vol. 7, p. 970, 2023.
  * [8] C. Kokail, C. Maier, R. Van Bijnen, T. Brydges, M. K. Joshi, P. Jurcevic, C. A. Muschik, P. Silvi, R. Blatt, C. F. Roos, and et al., “Self-verifying variational quantum simulation of lattice models,” Nature, vol. 569, no. 7756, p. 355–360, 2019.
  * [9] J. Kattemölle and J. Van Wezel, “Variational quantum eigensolver for the Heisenberg antiferromagnet on the Kagome lattice,” Physical Review B, vol. 106, no. 21, 2022.
  * [10] N. Moll, P. Barkoutsos, L. S. Bishop, J. M. Chow, A. Cross, D. J. Egger, S. Filipp, A. Fuhrer, J. M. Gambetta, M. Ganzhorn, and et al., “Quantum optimization using variational algorithms on near-term quantum devices,” Quantum Science and Technology, vol. 3, no. 3, p. 030503, 2018.
  * [11] E. Farhi, J. Goldstone, and S. Gutmann, “A quantum approximate optimization algorithm,” arXiv pre-print server, 2014.
  * [12] Y. Miao, J. Lamers, and V. Pasquier, “On the q operator and the spectrum of the xxz model at root of unity,” SciPost Physics, vol. 11, no. 3, 2021\.
  * [13] L. Šamaj and Z. Bajnok, XXZ Heisenberg chain: Bethe ansatz and the ground state, p. 158–174.  Cambridge University Press, 2013.
  * [14] J. Yang, X. Wang, and J. Wu, “Magnetic excitations in the one-dimensional heisenberg–ising model with external fields and their experimental realizations,” Journal of Physics A: Mathematical and Theoretical, vol. 56, no. 1, p. 013001, 2023.
  * [15] A. Cinnirella, Dynamics of the Heisenberg XXZ model across the ferromagnetic transition.  PhD thesis.
  * [16] “2.5: Expectation values — chem.libretexts.org.” https://chem.libretexts.org/Courses/Lock_Haven_University_of_Pennsylvania/CHEM321%3A_Quantum_Chemistry_and_Spectroscopy/Quantum_Chemistry_and_Spectroscopy/2%3A_The_Postulates_of_Quantum_Mechanics/2.5%3A_Expectation_values.  Accessed 22-08-2023.
  * [17] C. D. Cushen and R. L. Hudson, “A quantum-mechanical central limit theorem,” Journal of Applied Probability, vol. 8, no. 3, pp. 454–469, 1971.
  * [18] N. C. Rubin, R. Babbush, and J. Mcclean, “Application of fermionic marginal constraints to hybrid quantum algorithms,” New Journal of Physics, vol. 20, no. 5, p. 053020, 2018.
  * [19] A. Arrasmith, L. Cincio, Rolando, and Patrick, “Operator sampling for shot-frugal optimization in variational algorithms,” arXiv pre-print server, 2020.

## Appendix

###  4.1 Exact Computation Optimization Records

The ground state optimization for the isotropic Heisenberg model is obtained
via the exact computation of the expectation values. When compared with the
results obtained from sampling, the results are more accurate but also need
more iterations to converge.

Figure 8: Optimization record of 1D Heisenberg Model for N =
4,6,8,10,12,14,16, with classically computed expectation value Figure 9:
Ground states of the 1D Anisotropic XXZ Heisenberg model with in the range of
[-1,1] sampled at an interval length of 0.1, classically computed expectation
value

###  4.2 Prepared State Verification

####  4.2.1 Correlation Function

The XXZ model is in the critical phase when the anisotropy parameter satisfies
$-1\leq\Delta\leq+1$. This is characterized by long-range correlations
(polynomially decaying correlations). For $\Delta\geq 1$, the model is in a
gapped phase characterized by exponentially decaying correlations. These
trends can be seen from the correlation function plots below.

Figure 10: the correlation of the prepared state for $N=10$ and $\Delta=0.3$
with respect to r(distance between qubits) Figure 11: the correlation of the
prepared state for $N=10$ and $\Delta=1.0$ with respect to r Figure 12: the
correlation of the prepared state for $N=10$ and $\Delta=2.5$ with respect to
r

####  4.2.2 Subsystem Entropy

For $\Delta=1$ , the subsystem entropy displays an alternating profile
characteristic of the dimer-like nature of the ground state of the Heisenberg
model. This vanishes as one goes away from the isotropic limit.

Figure 13: the subsystem entropy of the prepared states for $N=10$ and
$\Delta=0.3$ Figure 14: the subsystem entropy of the prepared states for
$N=10$ and $\Delta=1.0$ Figure 15: the subsystem entropy of the prepared
states for $N=10$ and $\Delta=2.5$