# Error estimates for the scalar auxiliary variable (SAV) scheme to the Cahn-
Hilliard equation

Shu Ma Department of Mathematics, City University of Hong Kong, 83 Tat Chee
Avenue, Kowloon, Hong Kong, P.R. China. Email address<EMAIL_ADDRESS>Weifeng Qiu Department of Mathematics, City University of Hong Kong, 83 Tat
Chee Avenue, Kowloon, Hong Kong, P.R. China. Email address:
<EMAIL_ADDRESS>
###### Abstract

The optimal error estimate that depending only on the polynomial degree of
$\varepsilon^{-1}$ is established for the temporal semi-discrete scheme of the
Cahn-Hilliard equation, which is based on the scalar auxiliary variable (SAV)
formulation. The key to our analysis is to convert the structure of the SAV
time-stepping scheme back to a form compatible with the original format of the
Cahn-Hilliard equation, which makes it feasible to use spectral estimates to
handle the nonlinear term. Based on the transformation of the SAV numerical
scheme, the optimal error estimate for the temporal semi-discrete scheme which
depends only on the low polynomial order of $\varepsilon^{-1}$ instead of the
exponential order, is derived by using mathematical induction, spectral
arguments, and the superconvergence properties of some nonlinear terms.
Numerical examples are provided to illustrate the discrete energy decay
property and validate our theoretical convergence analysis.

Key words:  Cahn-Hilliard equation, SAV formulation, energy decay, spectral
estimates, polynomial order, error estimates.

##  1 Introduction

In this paper we consider the initial boundary value problem for the Cahn-
Hilliard (CH) phase field equation

$\displaystyle\partial_{t}u$ $\displaystyle=\Delta\big{(}-\varepsilon\Delta
u+\frac{1}{\varepsilon}f(u)\big{)}$
$\displaystyle\qquad\mbox{in}\,\,\,{\varOmega}\times(0,T],$ (1.1a)
$\displaystyle\partial_{\bm{n}}u$
$\displaystyle=\partial_{\bm{n}}\big{(}-\varepsilon\Delta
u+\frac{1}{\varepsilon}f(u)\big{)}=0$
$\displaystyle\qquad\mbox{on}\,\,\,\partial{\varOmega}\times(0,T],$ (1.1b)
$\displaystyle u(\cdot,0)$ $\displaystyle=u_{0},$
$\displaystyle\qquad\mbox{in}\,\,\,{\varOmega},$ (1.1c)

where ${\varOmega}\subset R^{d}$, $d=2,3$ is a bounded domain, ${\bm{n}}$ is
the outward normal, $\varepsilon$ is a small parameter, and $f$ is the
derivative of a non-negative potential function $F$ with two local minima,
i.e., $f=F^{\prime}$. For instance, the Ginzburg-Landau energy function

$F(v)=\frac{1}{4}(v^{2}-1)^{2}\quad\mbox{and}\quad f(v)=v^{3}-v.$

In view of its wide application as a phase field model [10, 38], many
numerical methods and analyses have been developed for approximating the Cahn-
Hilliard equation (1.1). On the one hand, most of these works have been
performed for the Cahn-Hilliard equation with fixed $\varepsilon>0$. As
pointed out in [4, 15, 17, 18], error estimates using the direct Gronwall
inequality argument yield a constant factor $e^{T/\varepsilon}$, and as a
result the error grows exponentially as $\varepsilon\to 0$. Such an estimate
is clearly not useful for very small $\varepsilon$, especially in solving the
problem of whether the computed numerical interface converges to the original
sharp interface of the Hele-Shaw problem when $\varepsilon\to 0$, see [14, 27]
for details. To overcome this difficulty, Feng and Prohl [24] first
established a priori error estimates with polynomial dependence on
$\varepsilon^{-1}$ for the time-discrete format of the Cahn-Hilliard equation.
Then, Feng and Wu obtained a posteriori error estimates with polynomial-order
of $\varepsilon^{-1}$ in [25] for the same time-discrete methods. The main
idea of the polynomial-order error estimates of $\varepsilon^{-1}$ is to use
the spectral estimates given by Alikakos and Fusco [3] and Chen [13] for the
linearized Cahn-Hilliard operator to handle the nonlinear term in the error
analysis. After this, spectral estimates were frequently used to eliminate the
exponential dependence on $\varepsilon^{-1}$ in the error analyses of other
numerical methods for the Cahn-Hilliard equations (see [6, 19, 21, 28, 33] and
references therein) and the related phase field equations, including the
Allen-Cahn equations in [1, 6, 7, 8, 20, 22, 23, 26], the Ginzburg-Landau
equations in [5], and the phase field models with nonlinear constitutive laws
in [16].

On the other hand, it is well-know that the Cahn-Hilliard equation (1.1) is
the $H^{-1}$-gradient flow of the energy functional

$\displaystyle E[u]:=\int_{{\varOmega}}\Big{(}\frac{\varepsilon}{2}|\nabla
u|^{2}+\frac{1}{\varepsilon}F(u)\Big{)}\,{\mathrm{d}}x.$ (1.2)

As a result, the solution of the Cahn-Hilliard equation has decaying energy.
Indeed, testing (1.1) by $-\varepsilon\Delta u+\frac{1}{\varepsilon}f(u)$
yields

$\displaystyle
E[u(\cdot,t_{2})]-E[u(\cdot,t_{1})]=-\int_{t_{1}}^{t_{2}}\|\partial_{t}u\|_{H^{-1}}^{2}\,{\mathrm{d}}t\leq
0,$ (1.3)

where $0\leq t_{1}\leq t_{2}\leq T$ and $\|\partial_{t}u\|_{H^{-1}}=\|\nabla
w\|$ with $w=-\varepsilon\Delta u+\frac{1}{\varepsilon}f(u)$.

Accordingly, great efforts have been devoted to the construction of efficient
and accurate numerical methods that preserve the energy decay properties at
the discrete level. In particular, for those widely used linear time-stepping
methods with energy decay, including stabilized semi-implicit schemes [11, 12,
32], invariant energy quadratization (IEQ) methods [34, 35, 36], and the
scalar auxiliary variable (SAV) approach in [29, 30] and [2]. Optimal error
estimates of these time-stepping schemes for the Cahn-Hilliard equations with
fixed $\varepsilon^{-1}$ are now well developed. As far as the error analysis
is concerned, the main difficulty that remains with these methods is how to
establish error bounds that depend only on the polynomial order of
$\varepsilon^{-1}$ rather than the exponential order for small $\varepsilon\to
0$. This is due to the fact that compared to Feng’s previous work [19, 21, 24,
25, 28], the numerical methods based on the IEQ/SAV formulation break the
standard structure of the nonlinear term of the Cahn-Hilliard equation, which
is crucial to the utilization of the spectral arguments. This makes the
spectral estimates ineffective in estimating errors for the IEQ/SAV approach.
To the best of our knowledge, the optimal error estimates depending only on
the polynomial order of $\varepsilon^{-1}$ for the IEQ/SAV methods to the
Cahn-Hilliard equation remains open.

The objective of this paper is to establish error bounds which depend on
$\varepsilon^{-1}$ only in low polynomial order for a semi-discrete methods
based on the SAV formulation of the Cahn-Hilliard equation. The SAV
reformulation of the Cahn-Hilliard equation was introduced in [30, 31] as an
enhanced version of the invariant energy quadratization (IEQ) approach [34,
35, 36, 37], for developing energy-decay methods at the discrete level. By
reconstructing the system based on the SAV reformulation, we obtain a time
semi-discrete scheme, which is linear and easy-to-implement. The SAV
formulation introduces new difficulties to the error analysis for the Cahn-
Hilliard equation due to the presence of a new scalar $r$ in the nonlinear
part (see equation (2.8b)), which alters the structure of the original Cahn-
Hilliard equation and makes the spectral argument not directly applicable. To
improve the current error analysis, Zhang and Yang [39] have recently made a
breakthrough in the estimates of the IEQ method for the Allen-Cahn equation.
They established optimal error bounds on the polynomial dependence of
$\varepsilon^{-1}$ for the IEQ-based numerical schemes. Inspired by [39], we
rewrite the format of the SAV scheme using the new scalar variable into a form
compatible with the original format of the Cahn-Hilliard equation, which makes
it feasible to use spectral estimates in the error analysis. However, this
structural transformation will accordingly introduces a strong perturbation
term that needs to be delicately controlled. Furthermore, unlike the Allen-
Cahn equation which is a gradient flow in $L^{2}({\varOmega})$, the Cahn-
Hilliard equation is a gradient flow in $H^{-1}({\varOmega})$, which makes the
analysis for the Cahn-Hilliard equation in this paper more delicate and
complicated than that for the Allen-Cahn equation. In our analysis, these
difficulties are overcome by combining the following techniques:

  * (1)

To use the spectral estimates of the linearized Cahn-Hilliard operator to deal
with the nonlinear potential term in the error analysis, we reconvert the
structure of the SAV scheme into a form compatible with the original Cahn-
Hilliard equation (1.1).

  * (2)

An inductive argument is used to deal with the difficulties caused by the
strong perturbation term that appear in the structural transformation (see
equation (4.49)). In particular, error bounds of $\|\nabla e^{i-1}-\nabla
e^{i-2}\|$ and $\|e^{i-1}-e^{i-2}\|_{H^{-1}}$ for $i\leq n$ need to be
established. By integrating them into the estimates of the perturbation terms,
the super-convergence characteristics of some of their resulting nonlinear
terms will complete the mathematical induction method.

  * (3)

Given that the solution of the Cahn-Hilliard problem (1.1) preserves the total
mass property (i.e.
$\frac{{\mathrm{d}}}{{\mathrm{d}}t}\int_{{\varOmega}}u(x,t)\,{\mathrm{d}}x=0$),
which is not possessed by the corresponding Allen-Cahn problem, we can
establish its optimal error bounds in the
$L^{\infty}(0,T;H^{-1}({\varOmega}))$-norm with polynomial dependence on
$\varepsilon^{-1}$, with the help of an error estimate in the
$L^{2}(0,T;H^{1}({\varOmega}))$-norm. This is different from the error
estimates in the $L^{\infty}(0,T;L^{2}({\varOmega}))$-norm of the Allen-Cahn
equation given in [39].

To the best of our knowledge, this is the first error estimate of the
polynomial dependence on $\varepsilon^{-1}$ for the IEQ/SAV-type schemes of
the Cahn-Hilliard equation (1.1).

The rest of this paper is organized as follows. In Section 2, we present the
SAV reformulation of the Cahn-Hilliard equation and introduce an equivalent
transformation of the SAV time-stepping scheme. In Section 3, we show the
properties of energy decay and derive the consistency estimates for the
proposed method. In Section 4, we present an error estimate of the semi-
discrete SAV scheme to derive a convergence rate that does not depend on
$\varepsilon^{-1}$ exponentially. The spectrum estimate plays a crucial role
in the proof. Finally, in Section 5, we present a few numerical experiments to
validate the theoretical results.

##  2 Formulation of the Semi-discrete SAV scheme

In this section, we construct a backward Euler implicit-explicit type temporal
semi-discrete numerical scheme based on the SAV reformulation of the CH
equation (1.1), and also present an equivalent formulation of the SAV scheme.

###  2.1 Function spaces

Let $W^{s,p}({\varOmega})$ denote the usual Sobolev spaces, and
$H^{s}({\varOmega})$ denote the Hilbert spaces $W^{s,2}({\varOmega})$ with
norm $\|\cdot\|_{H^{s}}$. Let $\|\cdot\|$ and $(\cdot,\cdot)$ represent the
$L^{2}$ norm and $L^{2}$ inner product, respectively. In addition, define for
$p\geq 0$

$\displaystyle H^{-p}({\varOmega}):=(H^{p}({\varOmega}))^{*},\quad
H_{0}^{-p}({\varOmega}):=\\{u\in H^{-p}({\varOmega})\,|\,\langle
u,1\rangle_{p}=0\\},$ (2.4)

where $\langle\cdot,\cdot\rangle_{p}$ stands for the dual product between
$H^{p}({\varOmega})$ and $H^{-p}({\varOmega})$. We denote
$L_{0}^{2}({\varOmega}):=H_{0}^{0}({\varOmega})$. For $v\in
L_{0}^{2}({\varOmega})$, let $-\Delta^{-1}v:=v_{1}\in H^{1}({\varOmega})\cap
L_{0}^{2}({\varOmega})$, where $v_{1}$ is the solution to

$\displaystyle-\Delta v_{1}=v\,\,\mbox{in}\,\,{\varOmega},\qquad\frac{\partial
v_{1}}{\partial n}=0\,\,\mbox{on}\,\,\partial{\varOmega},$ (2.5)

and $\|v\|_{-1}:=\sqrt{(v,-\Delta^{-1}v)}$.

For $v\in L_{0}^{2}({\varOmega})\cap H^{1}({\varOmega})$, we have the
following inequality

$\displaystyle\|v\|^{2}=(\nabla v,\nabla(-\Delta)^{-1}v)\leq\|\nabla
v\|\|v\|_{-1}.$ (2.6)

We denote by $C$ generic constant and $C_{i}$, $\tilde{C}_{i}$, $\tilde{C}$,
$\kappa_{i}^{\prime}$ and $\kappa_{i}$ specific constants, which are
independent of $\tau$, $h$ and $\varepsilon$, but may possibly depend on the
domain ${\varOmega}$, $T$ and the constants of Sobolev inequalities. We use
notation $\lesssim$ in the sense that $f\lesssim g$ means that $f\leq Cg$ with
positive constant $C$ independent of $\tau$, $h$ and $\varepsilon$.

###  2.2 The SAV reformulation

The SAV formulation of the CH equation (cf. [29, 30]) introduces a scalar
auxiliary variable

$\displaystyle
r=\sqrt{\mbox{$\int_{\varOmega}$}F(u){\mathrm{d}}x+c_{0}}\quad\mbox{with}\quad
g(u)=\frac{f(u)}{\sqrt{\mbox{$\int_{\varOmega}$}F(u){\mathrm{d}}x+c_{0}}},$
(2.7)

with a positive $c_{0}$ (which guarantees that the function $r$ has a positive
lower bound), and reformulate (1.1) as

$\displaystyle\partial_{t}u$ $\displaystyle=\Delta w$
$\displaystyle\qquad\mbox{in}\,\,\,{\varOmega}\times(0,T],$ (2.8a)
$\displaystyle w$ $\displaystyle=-\varepsilon\Delta
u+\frac{1}{\varepsilon}rg(u)$
$\displaystyle\qquad\mbox{in}\,\,\,{\varOmega}\times(0,T],$ (2.8b)
$\displaystyle\frac{{\mathrm{d}}r}{{\mathrm{d}}t}$
$\displaystyle=\frac{1}{2}\big{(}g(u),\partial_{t}u\big{)}$
$\displaystyle\qquad\mbox{in}\,\,\,{\varOmega}\times(0,T],$ (2.8c)
$\displaystyle\partial_{\bm{n}}u$ $\displaystyle=\partial_{\bm{n}}w=0$
$\displaystyle\qquad\mbox{on}\,\,\,\partial{\varOmega}\times(0,T],$ (2.8d)
$\displaystyle u(\cdot,0)$ $\displaystyle=u_{0}$
$\displaystyle\qquad\mbox{in}\,\,\,{\varOmega},$ (2.8e) $\displaystyle r(0)$
$\displaystyle=\sqrt{\mbox{$\int_{\varOmega}$}F(u_{0}){\mathrm{d}}x+c_{0}}.$
(2.8f)

We define an energy functional with respect to $u$ and $r$:

$\displaystyle E[u,r]=\frac{\varepsilon}{2}\|\nabla
u\|^{2}+\frac{1}{\varepsilon}r^{2},$ (2.9)

and taking the $L^{2}$ inner product of the first equation (2.8a) with $w$, of
the second equation (2.8b) with $\partial_{t}u$, and of the third equation
(2.8c) with $\frac{1}{\varepsilon}2r$, performing integration by parts and
summing up the two obtained equations, we find that

$\displaystyle\frac{{\mathrm{d}}}{{\mathrm{d}}t}E[u,r]=-\|\nabla
w\|^{2}=-\|\partial_{t}u\|_{-1}^{2}\leq 0.$ (2.10)

###  2.3 The equivalent formulation of the SAV scheme

Let $\\{t_{n}\\}_{n=0}^{N}$ be a uniform partition of $[0,T]$ with the time
step size $\tau$, where $N$ is a positive integer and hence
$\tau=\frac{T}{N}$. We consider the following temporal semi-discrete SAV
scheme for solving the system (2.8):

$\left\\{\begin{aligned} \frac{u^{n+1}-u^{n}}{\tau}&=\Delta w^{n+1},\\\
w^{n+1}&=-\varepsilon\Delta u^{n+1}+\frac{1}{\varepsilon}r^{n+1}g(u^{n}),\\\
r^{n+1}-r^{n}&=\frac{1}{2}\big{(}g(u^{n}),u^{n+1}-u^{n}\big{)},\\\
\partial_{\bm{n}}u^{n+1}|_{\partial{\varOmega}}&=\partial_{\bm{n}}w^{n+1}|_{\partial{\varOmega}}=0,\end{aligned}\right.$
(2.11)

with $u^{0}=u_{0}$ and $r^{0}=r(0)$ for $n=0,1,\dots,N$.

By taking the $L^{2}$ inner product of the first equation in (2.11) with
$v=1$, we have the following volume conservation property, which is important
to the error estimates.

###### Lemma 2.1

The numerical solution of (2.11) satisfies

$\displaystyle\frac{1}{|{\varOmega}|}\int_{{\varOmega}}u^{n}(x)\,{\mathrm{d}}x=\frac{1}{|{\varOmega}|}\int_{{\varOmega}}u^{0}(x)\,{\mathrm{d}}x,\qquad
n=1,\dots,N,$ (2.12)

and the error function $e^{n}:=u(t_{n})-u^{n}$ satisfies

$\displaystyle\int_{{\varOmega}}e^{n}\,{\mathrm{d}}x=0,\qquad\qquad
n=1,\dots,N.$ (2.13)

Because of (2.12), $u^{n+1}-u^{n}$ and $e^{n}$ belong to
$L_{0}^{2}({\varOmega})$ such that we can define their $\|\cdot\|_{-1}$ norm.
Also, the semi-discrete SAV scheme (2.11) can be written as

$\left\\{\begin{aligned}
\Delta^{-1}\frac{u^{n+1}-u^{n}}{\tau}&=-\varepsilon\Delta
u^{n+1}+\frac{1}{\varepsilon}r^{n+1}g(u^{n}),\\\
r^{n+1}-r^{n}&=\frac{1}{2}\big{(}g(u^{n}),u^{n+1}-u^{n}\big{)},\qquad
n=0,1,\dots,N.\end{aligned}\right.$ (2.14)

In order to avoid the exponentially dependence of the error bound on
$\frac{1}{\varepsilon}$ induced by using the Gronwall inequality, we need to
use a spectral estimate of the linearized Cahn-Hilliard operator, which is
given in [3, 13, 24] and will be described in Section 4.

However, compared with the previous work [24], the SAV method (2.11) alters
the structure of the CH equation such that the spectral argument can not be
applied directly. To achieve the ideal error bound, we need to transform the
structure into a form compatible with the CH equation (1.1) so that the
spectral estimate of the linearized Cahn-Hilliard operator can be used. To
this end, we define $A(v)=\sqrt{\int_{{\varOmega}}F(v)\,{\mathrm{d}}x+c_{0}}$,
the Gateaux derivatives of $A(v)$ can be defined as follows:

$\displaystyle DA(v,w)$
$\displaystyle:=\frac{1}{2}\Big{(}\frac{f(v)}{\sqrt{\int_{{\varOmega}}F(v)\,{\mathrm{d}}x+c_{0}}},w\Big{)}=\frac{1}{2}\big{(}g(v),w\big{)},$
(2.15) $\displaystyle D^{2}A(v,w)$
$\displaystyle:=\frac{1}{2}\frac{\int_{{\varOmega}}f^{\prime}(v)w^{2}\,{\mathrm{d}}x}{\sqrt{\int_{{\varOmega}}F(v)\,{\mathrm{d}}x+c_{0}}}-\frac{1}{4}\frac{\big{(}\int_{{\varOmega}}f(v)w\,{\mathrm{d}}x\big{)}^{2}}{\Big{(}\int_{{\varOmega}}F(v)\,{\mathrm{d}}x+c_{0}\Big{)}^{\frac{3}{2}}}.$
(2.16)

By Taylor expansion, we derive

$\displaystyle
A(u^{i})=A(u^{i-1})+\frac{1}{2}\big{(}g(u^{i-1}),u^{i}-u^{i-1}\big{)}+\frac{1}{2}D^{2}A(\xi_{i}\,;u^{i}-u^{i-1}),$
(2.17)

where $\xi_{i}=\theta u^{i}+(1-\theta)u^{i-1}$ with $\theta\in(0,1)$. Thus we
get

$\displaystyle\frac{1}{2}\big{(}g(u^{i-1}),u^{i}-u^{i-1}\big{)}=A(u^{i})-A(u^{i-1})-\frac{1}{2}D^{2}A(\xi_{i}\,;u^{i}-u^{i-1}),$

which together with the second equation in (2.14) implies

$\displaystyle
r^{i}-r^{i-1}=A(u^{i})-A(u^{i-1})-\frac{1}{2}D^{2}A(\xi_{i}\,;u^{i}-u^{i-1}).$

After summing up the above equation from $i=1$ to $n$, we derive

$\displaystyle
r^{n}-r^{0}=A(u^{n})-A(u^{0})-\frac{1}{2}\sum_{i=1}^{n}D^{2}A(\xi_{i}\,;u^{i}-u^{i-1}).$

Since $r^{0}=A(u^{0})$, we obtain

$\displaystyle
r^{n}=A(u^{n})-\frac{1}{2}\sum_{i=1}^{n}D^{2}A(\xi_{i}\,;u^{i}-u^{i-1}).$
(2.18)

Then the SAV scheme (2.14) can be written as

$\displaystyle\Delta^{-1}\frac{u^{n+1}-u^{n}}{\tau}=\,$
$\displaystyle-\varepsilon\Delta
u^{n+1}+\frac{1}{\varepsilon}r^{n}g(u^{n})+\frac{1}{\varepsilon}g(u^{n})(r^{n+1}-r^{n})$
(2.19) $\displaystyle=\,$ $\displaystyle-\varepsilon\Delta
u^{n+1}+\frac{1}{\varepsilon}r^{n}g(u^{n})+\frac{1}{2\varepsilon}g(u^{n})\big{(}g(u^{n}),u^{n+1}-u^{n}\big{)},$

which together with (2.18) gives

$\displaystyle\Delta^{-1}\frac{u^{n+1}-u^{n}}{\tau}=\,$
$\displaystyle-\varepsilon\Delta
u^{n+1}+\frac{1}{\varepsilon}f(u^{n})-\frac{1}{2\varepsilon}g(u^{n})\sum_{i=1}^{n}D^{2}A(\xi_{i}\,;u^{i}-u^{i-1})$
(2.20)
$\displaystyle+\frac{1}{2\varepsilon}g(u^{n})\big{(}g(u^{n}),u^{n+1}-u^{n}\big{)}.$

The above equation (2.20) provides an equivalent formulation of the semi-
discrete SAV scheme (2.11), which will be frequently used in the subsequent
error analysis.

##  3 Energy decay and consistency analysis

In this subsection, we present several inequalities related to the proposed
numerical method.

###  3.1 Assumption and regularity

The quartic growth of the Ginzburg-Landau energy function
$F(u)=\frac{1}{4}(u^{2}-1)^{2}$ at infinity poses various technical
difficulties for the analysis and approximation of CH equations. Although the
CH equation does not satisfy the maximum principle, if the maximum norm of the
initial condition $u^{0}$ is bounded by $M$, it has been shown in [9] that the
maximum norm of the solution of the CH equation for the truncation potential
$F(u)$ with quadratic growth rate at infinity is bounded. Therefore, it has
been a common practice (cf. [32]) to consider the CH equations with a
truncated $F(u)$. Before presenting the detailed numerical analysis, we first
make some assumptions.

###### Assumption 3.1

We assume that the potential function $F(u)$ whose derivative
$f(u)=F^{\prime}(u)$ satisfies the following condition:

  * (i)

$F\in C^{4}(\mathbb{R})$, $F(\pm 1)=0$, and $F>0$ elsewhere.

  * (ii)

$f=F^{\prime}$, there exists a non-negative constant $L$ such that

$\displaystyle\max_{v\in\mathbb{R}}|f(v)|\leq
L,\quad\max_{v\in\mathbb{R}}|f^{\prime}(v)|\leq
L\quad\mbox{and}\quad\max_{v\in\mathbb{R}}|f^{\prime\prime}(v)|\leq L.$ (3.21)

In order to trace the dependence of the solution on the small parameter
$\varepsilon>0$, we assume that the solution of (1.1) satisfies the following
conditions:

###### Assumption 3.2

Suppose there exist positive $\varepsilon$-independent constants $m_{0}$ and
$\rho_{j}$ for $j=1,2,3$ such that the solution of (1.1) satisfies

$\displaystyle\frac{1}{|{\varOmega}|}\int_{{\varOmega}}u\,{\mathrm{d}}x=m_{0}\in(-1,1),$
(3.22a)
$\displaystyle\operatorname*{ess\,sup}_{t\in[0,\infty]}\Big{\\{}\frac{\varepsilon}{2}\|\nabla
u\|^{2}+\frac{1}{\varepsilon}\|F(u)\|_{L^{1}}\Big{\\}}+\int_{0}^{\infty}\|u_{t}\|_{H^{-1}}^{2}\,{\mathrm{d}}s\lesssim\varepsilon^{-\rho_{1}},$
(3.22b)
$\displaystyle\int_{0}^{\infty}\|u_{t}\|^{2}\,{\mathrm{d}}s\lesssim\varepsilon^{-\rho_{2}},$
(3.22c)
$\displaystyle\int_{0}^{\infty}\|\Delta^{-1}u_{tt}\|_{H^{-1}}^{2}\,{\mathrm{d}}t\lesssim\varepsilon^{-\rho_{3}}.$
(3.22d)

###### Remark 3.1

(a) Note that the conditions (i) and (ii) are satisfied by restricting the
growth of $F(u)$ for $|u|\geq M$. More precisely, for a given $M$, we can
replace $F(v)=\frac{1}{4}(v^{2}-1)^{2}$ by a cut-off function $\hat{F}(v)\in
C^{4}(\mathbb{R})$ as follows:

$\hat{F}(v)=\left\\{\begin{aligned}
&((2M)^{2}-1)2M(v-2M)+\frac{1}{4}((2M)^{2}-1)^{2}&&\mbox{for}\,\,v>2M,\\\
&\Phi_{+}(v)&&\mbox{for}\,\,v\in(M,2M],\\\
&\frac{1}{4}(v^{2}-1)^{2}&&\mbox{for}\,\,v\in[-M,M],\\\
&\Phi_{-}(v)&&\mbox{for}\,\,v\in[-2M,-M),\\\
&-((2M)^{2}-1)2M(v+2M)+\frac{1}{4}((2M)^{2}-1)^{2}&&\mbox{for}\,\,v<-2M.\end{aligned}\right.$
(3.23)

and $\Phi_{+}(v)$ and $\Phi_{-}(v)$ are the unique ninth order polynomials
satisfying the required conditions at $|v|=M$ and $|v|=2M$, respectively. Then
we replace $f(v)=(v^{2}-1)v$ by $\hat{F}^{\prime}(v)$ which is

$\hat{f}(v)=\hat{F}^{\prime}(v)=\left\\{\begin{aligned}
&((2M)^{2}-1)2M&&\mbox{for}\,\,v\geq 2M,\\\
&\Phi^{\prime}_{+}(v)&&\mbox{for}\,\,v\in(M,2M],\\\
&(v^{2}-1)v&&\mbox{for}\,\,v\in[-M,M],\\\
&\Phi^{\prime}_{-}(v)&&\mbox{for}\,\,v\in[-2M,-M),\\\
&-((2M)^{2}-1)2M&&\mbox{for}\,\,v\leq-2M.\end{aligned}\right.$ (3.24)

In simplicity, we still denote the modified function $\hat{F}$ by $F$. It is
then obvious that there exists $L$ such that (3.21) are satisfied with $f$
replaced by $\hat{f}$.

(b) The transformed SAV scheme (2.20) introduced a complicated term. With the
condition (ii) in the Assumption 3.1, we can get

$\displaystyle
D^{2}A(v;w)\leq\big{(}\|f^{\prime}(v)\|_{L^{\infty}}+\|f(v)\|^{2}\big{)}\|w\|^{2}\lesssim\|w\|^{2},$
(3.25)

which will be frequently used to control the difficult term in the error
analysis.

(c) Assumption 3.2 can be achieved in many cases. For example, suppose that
$f$ satisfies Assumption 3.1, $\partial{\varOmega}$ is of class $C^{2,1}$,
$u^{0}\in H^{3}({\varOmega})$, and there exist positive
$\varepsilon$-independent constants $\sigma_{j}$ for $j=1,2,3$ such that

$\displaystyle E(u^{0})=\frac{\varepsilon}{2}\|\nabla
u^{0}\|^{2}+\frac{1}{\varepsilon}\|F(u^{0})\|_{L^{1}}\lesssim\varepsilon^{-2\sigma_{1}},$
(3.26a) $\displaystyle\|-\varepsilon\Delta
u^{0}+\frac{1}{\varepsilon}f(u^{0})\|\lesssim\varepsilon^{-2\sigma_{2}},$
(3.26b) $\displaystyle\|-\varepsilon\Delta
u^{0}+\frac{1}{\varepsilon}f(u^{0})\|_{H^{1}}\lesssim\varepsilon^{-2\sigma_{3}}.$
(3.26c)

Then the estimates (3.22a)–(3.22d) can be derived by standard test function
techniques and satisfy:

$\displaystyle\rho_{1}=2\sigma_{1},\quad\rho_{2}=\max\\{2\sigma_{1}+2,2\sigma_{2}-1\\}\quad\mbox{and}\quad\rho_{3}=\max\\{2\sigma_{1}+6,2\sigma_{2}+3,2\sigma_{3}-1\\}.$

We refer to [24] for their detailed proof.

###  3.2 Energy decay structure

In this subsection, we prove the following energy decay property of the
numerical solution, which comprise of the first theorem of this paper.

###### Theorem 3.1

(energy decay) The scheme (2.14) is unconditionally energy stable in the sense
that

$\displaystyle
E(u^{n+1},r^{n+1})-E(u^{n},r^{n})\leq-\frac{1}{\tau}\|u^{n+1}-u^{n}\|_{-1}^{2}\leq
0\quad\mbox{for }n\geq 1.$ (3.27)

Proof.  Taking the inner product of the first equation in (2.14) with
$u^{n+1}-u^{n}\in L_{0}^{2}({\varOmega})$ and multiplying the second equation
in (2.14) by $\frac{2}{\varepsilon}r^{n+1}$, we derive that

$\displaystyle\frac{1}{\tau}\|u^{n+1}-u^{n}\|_{-1}^{2}+\frac{\varepsilon}{2}\Big{(}\|\nabla
u^{n+1}\|^{2}-\|\nabla u^{n}\|^{2}+\|\nabla u^{n+1}-\nabla u^{n}\|^{2}\Big{)}$
(3.28)
$\displaystyle\qquad\qquad+\frac{1}{\varepsilon}r^{n+1}\big{(}g(u^{n}),u^{n+1}-u^{n}\big{)}=0,$
$\displaystyle\frac{1}{\varepsilon}\Big{(}(r^{n+1})^{2}-(r^{n})^{2}+(r^{n+1}-r^{n})^{2}\Big{)}=\frac{1}{\varepsilon}r^{n+1}\big{(}g(u^{n}),u^{n+1}-u^{n}\big{)}.$
(3.29)

Taking the summation of the above obtained equations, we get

$\displaystyle\frac{1}{\tau}\|u^{n+1}-u^{n}\|_{-1}^{2}+$
$\displaystyle\frac{\varepsilon}{2}\Big{(}\|\nabla u^{n+1}\|^{2}-\|\nabla
u^{n}\|^{2}+\|\nabla u^{n+1}-\nabla u^{n}\|^{2}\Big{)}$ (3.30)
$\displaystyle+$
$\displaystyle\frac{1}{\varepsilon}\Big{(}(r^{n+1})^{2}-(r^{n})^{2}+(r^{n+1}-r^{n})^{2}\Big{)}=0.$

which gives (3.27). $\square$

###### Remark 3.2

After summing up (3.30) from $n=0$ to $N$, we get

$\displaystyle\frac{\varepsilon}{2}\|\nabla
u^{N+1}\|^{2}+\frac{1}{\varepsilon}(r^{N+1})^{2}$
$\displaystyle+\frac{1}{\tau}\sum_{n=0}^{N}\|u^{n+1}-u^{n}\|_{-1}^{2}+\frac{\varepsilon}{2}\sum_{n=0}^{N}\|\nabla
u^{n+1}-\nabla u^{n}\|^{2}$ (3.31)
$\displaystyle+\frac{1}{\varepsilon}\sum_{n=0}^{N}(r^{n+1}-r^{n})^{2}=\frac{\varepsilon}{2}\|\nabla
u^{0}\|^{2}+\frac{1}{\varepsilon}(r^{0})^{2}\lesssim\varepsilon^{-\rho_{1}}.$

It follows from $u^{n+1}-u^{n}\in L_{0}^{2}({\varOmega})$ and (2.6) that

$\displaystyle\sum_{n=0}^{N}\|u^{n+1}-u^{n}\|^{2}\leq\,$
$\displaystyle\sum_{n=0}^{N}\|u^{n+1}-u^{n}\|_{-1}\|\nabla u^{n+1}-\nabla
u^{n}\|$ (3.32) $\displaystyle\leq\,$
$\displaystyle\Big{(}\sum_{n=0}^{N}\|u^{n+1}-u^{n}\|_{-1}^{2}\Big{)}^{\frac{1}{2}}\Big{(}\sum_{n=0}^{N}\|\nabla
u^{n+1}-\nabla u^{n}\|^{2}\Big{)}^{\frac{1}{2}}$ $\displaystyle\lesssim\,$
$\displaystyle\varepsilon^{-(\rho_{1}+\frac{1}{2})}\tau^{\frac{1}{2}}.$

###  3.3 Consistency

To derive the error estimates of the equivalent transformation (2.20) of the
semi-discrete SAV scheme (2.11), we reformulate the CH equation (1.1) as the
truncated form

$\displaystyle\Delta^{-1}\frac{u(t_{n+1})-u(t_{n})}{\tau}$
$\displaystyle=-\varepsilon\Delta
u(t_{n+1})+\frac{1}{\varepsilon}f(u(t_{n}))+\mathcal{R}^{n+1},$ (3.33)

where the truncation error $\mathcal{R}^{n+1}$ is given by

$\displaystyle\mathcal{R}^{n+1}=\Big{[}\Delta^{-1}\frac{u(t_{n+1})-u(t_{n})}{\tau}-\Delta^{-1}\partial_{t}u(t_{n+1})\Big{]}+\frac{1}{\varepsilon}\Big{[}f(u(t_{n+1}))-f(u(t_{n}))\Big{]}.$
(3.34)

###### Lemma 3.1

(consistency estimate) If the solution of (1.1) satisfies (3.22c)–(3.22d),
then the following consistency estimate holds:

$\displaystyle\tau\sum_{n=0}^{N}\|\mathcal{R}^{n+1}\|_{H^{-1}}^{2}\leq
C\varepsilon^{-\max\\{\rho_{2}+2,\,\rho_{3}\\}}\tau^{2}.$ (3.35)

Proof.  By performing standard calculations, we get

$\displaystyle\Big{\|}\Delta^{-1}\frac{u(t_{n+1})-u(t_{n})}{\tau}-\Delta^{-1}\partial_{t}u(t_{n+1})\Big{\|}_{H^{-1}}^{2}=\,$
$\displaystyle\Big{\|}\int_{t_{n}}^{t_{n+1}}(s-t_{n})\Delta^{-1}u_{tt}\,{\mathrm{d}}s\Big{\|}_{H^{-1}}^{2}$
(3.36) $\displaystyle\lesssim\,$
$\displaystyle\tau\int_{t_{n}}^{t_{n+1}}\|\Delta^{-1}u_{tt}\|_{H^{-1}}^{2}\,{\mathrm{d}}s,$

and

$\displaystyle\Big{\|}\frac{1}{\varepsilon}f(u(t_{n+1}))-\frac{1}{\varepsilon}f(u(t_{n}))\Big{\|}_{H^{-1}}^{2}=\,$
$\displaystyle\varepsilon^{-2}\sup_{v\in
H^{1}({\varOmega})}\frac{\big{(}f(u(t_{n+1}))-f(u(t_{n})),v\big{)}^{2}}{\|v\|_{H^{1}}^{2}}$
(3.37) $\displaystyle\leq\,$ $\displaystyle\varepsilon^{-2}\sup_{v\in
H^{1}({\varOmega})}\frac{\|f^{\prime}(\xi^{n})\|_{L^{3}}^{2}\|u(t_{n+1})-u(t_{n})\|^{2}\|v\|_{L^{6}}^{2}}{\|v\|_{H^{1}}^{2}}$
$\displaystyle\lesssim\,$
$\displaystyle\varepsilon^{-2}\|u(t_{n+1})-u(t_{n})\|^{2}$
$\displaystyle\lesssim\,$
$\displaystyle\varepsilon^{-2}\tau\int_{t_{n}}^{t_{n+1}}\|u_{t}\|^{2}\,{\mathrm{d}}s,$

where $\xi^{n}$ is between $u(t_{n})$ and $u(t_{n+1})$. Thus, we have

$\displaystyle\tau\sum_{n=0}^{N}\|\mathcal{R}^{n+1}\|_{H^{-1}}^{2}\lesssim\,$
$\displaystyle\tau^{2}\int_{0}^{T}\|\Delta^{-1}u_{tt}\|_{H^{-1}}^{2}\,{\mathrm{d}}s+\varepsilon^{-2}\tau^{2}\int_{0}^{T}\|u_{t}\|^{2}\,{\mathrm{d}}s$
(3.38) $\displaystyle\lesssim\,$
$\displaystyle\varepsilon^{-\rho_{3}}\tau^{2}+\varepsilon^{-(\rho_{2}+2)}\tau^{2}.$

The proof is completed. $\square$

##  4 Error estimates

In this section, we will derive the error bound of the semi-discrete scheme
(2.11), in which the focus is to obtain the polynomial type dependence of the
error bound on $\varepsilon^{-1}$. If we use the usual error estimate of the
SAV numerical scheme (2.11), the error growth depends on $\varepsilon^{-1}$
exponentially. To avoid the exponential dependence on $\varepsilon^{-1}$
induced by using the Gronwall inequality, we need to use a spectral estimate
of the linearized Cahn-Hilliard operator, which is given in [3, 13, 24].

###### Lemma 4.1

(spectral estimate) Suppose that Assumption 3.1 holds. Then there exist
$0<\varepsilon_{0}<<1$ and a positive $\lambda_{0}$ such that the principle
eigenvalue of the linearized Cahn-Hilliard operator

$\displaystyle\mathcal{L}_{CH}:=\Delta(\varepsilon\Delta-\frac{1}{\varepsilon}f^{\prime}(u)I)$

satisfies for all $t\in[0,T]$

$\displaystyle\lambda_{CH}=\inf_{\begin{subarray}{c}0\neq v\in
H^{1}({\varOmega})\\\ \Delta w=v\end{subarray}}\frac{\varepsilon\|\nabla
v\|^{2}+\frac{1}{\varepsilon}\big{(}f^{\prime}(u(\cdot,t)v,v)\big{)}}{\|\nabla
w\|^{2}}\geq-\lambda_{0},$ (4.39)

for $\varepsilon\in(0,\varepsilon_{0})$, where $I$ denotes the identity
operator and $u$ is the solution of the Cahn-Hilliard problem (1.1).

We will now prove the following error estimates for the semi-discrete
numerical scheme, which is the main result of this paper.

###### Theorem 4.1

(error estimate) We assume that assumptions 3.1 and 3.2 hold and that

$\tau\leq\tilde{C}\varepsilon^{\beta_{0}}\quad\mbox{with}\quad\beta_{0}=\left\\{\begin{aligned}
&\alpha_{0}+8\quad&&\mbox{for}\quad d=2,\\\
&\frac{6\alpha_{0}+36}{5}\quad&&\mbox{for}\quad d=3,\end{aligned}\right.$
(4.40)

then the discrete solution given by (2.11) satisfies the following error
estimate for $e^{n}=u(t_{n})-u^{n}$:

$\displaystyle\max_{1\leq n\leq m}\|e^{n}\|_{-1}^{2}$
$\displaystyle+\frac{1}{2}\sum_{n=1}^{m}\|e^{n}-e^{n-1}\|_{-1}^{2}+\frac{1}{2}\varepsilon^{4}\sum_{n=1}^{m}\tau\|\nabla
e^{n}\|^{2}+\max_{1\leq n\leq m}\tau\varepsilon^{4}\|\nabla e^{m}\|^{2}$
(4.41) $\displaystyle+\frac{1}{2}\varepsilon^{4}\sum_{n=1}^{m}\tau\|\nabla
e^{n}-\nabla e^{n-1}\|^{2}\leq\kappa_{0}\varepsilon^{-\alpha_{0}}\tau^{2}.$

where
$\alpha_{0}:=\max\\{\rho_{1}+3,\,2\rho_{2}+4,\,\rho_{2}+6,\,\rho_{3}+4\\}$,
and the constants $\tilde{C}$ and $\kappa_{0}$ are given in the proof.

Proof.  We use the mathematical induction as follows. The proof is split into
four steps. The first step gives the error estimate for the first step
$t=t^{1}$. Steps two and three use the spectral estimate (4.39) to avoid
exponential blow-up in $\varepsilon^{-1}$ of the error constants. In the last
step, an inductive argument is used to conclude the proof.

Step 1: Estimation of $\|e^{1}\|_{-1}^{2}+\tau\varepsilon\|\nabla
e^{1}\|^{2}$. For $n=0$ in (2.14), we have

$\left\\{\begin{aligned}
\Delta^{-1}\frac{u^{1}-u^{0}}{\tau}&=-\varepsilon\Delta
u^{1}+\frac{1}{\varepsilon}r^{1}g(u^{0})\\\
r^{1}&=r^{0}+\frac{1}{2}\big{(}g(u^{0}),u^{1}-u^{0}\big{)}\end{aligned}\right.$
(4.42)

After plugging the second equation into the first equation, we get

$\displaystyle\Delta^{-1}\frac{u^{1}-u^{0}}{\tau}$
$\displaystyle=-\varepsilon\Delta
u^{1}+\frac{1}{\varepsilon}f(u^{0})+\frac{1}{2\varepsilon}g(u^{0})\big{(}g(u^{0}),u^{1}-u^{0}\big{)}.$
(4.43)

Subtracting (4.43) from (3.33), we get the corresponding error equation

$\displaystyle\Delta^{-1}\frac{e^{1}}{\tau}$ $\displaystyle=-\varepsilon\Delta
e^{1}+\mathcal{R}^{1}-\frac{1}{2\varepsilon}g(u^{0})\big{(}g(u^{0}),u^{1}-u^{0}\big{)}.$
(4.44)

Taking the $L^{2}$ inner product of the above equation with $e^{1}\in
L_{0}^{2}({\varOmega})$, we have

$\displaystyle\|e^{1}\|_{-1}^{2}+\tau\varepsilon\|\nabla
e^{1}\|^{2}=-\tau(\mathcal{R}^{1},e^{1})+\frac{\tau}{2\varepsilon}\big{(}g(u^{0}),u^{1}-u^{0}\big{)}\big{(}g(u^{0}),e^{1}\big{)}.$
(4.45)

Using Lemma 3.1, we have

$\displaystyle\tau(\mathcal{R}^{1},e^{1})\leq\frac{\varepsilon\tau}{2}\|\nabla
e^{1}\|^{2}+C\varepsilon^{-1}\tau\|\mathcal{R}^{1}\|_{H^{-1}}^{2}\leq\frac{\varepsilon\tau}{4}\|\nabla
e^{1}\|^{2}+C\varepsilon^{-\max\\{\rho_{2}+3,\,\rho_{3}+1\\}}\tau^{2}.$ (4.46)

From (3.32), we have

$\displaystyle\frac{\tau}{2\varepsilon}\big{(}g(u^{0}),u^{1}-u^{0}\big{)}\big{(}g(u^{0}),e^{1}\big{)}\leq\,$
$\displaystyle C\tau\varepsilon^{-1}\|g(u^{0})\|^{2}\|u^{1}-u^{0}\|\|e^{1}\|$
(4.47) $\displaystyle\leq\,$ $\displaystyle
C\tau\varepsilon^{-1}\|u^{1}-u^{0}\|\|e^{1}\|_{-1}^{\frac{1}{2}}\|\nabla
e^{1}\|^{\frac{1}{2}}$ $\displaystyle\leq\,$
$\displaystyle\frac{1}{2}\tau^{\frac{1}{2}}\varepsilon^{\frac{1}{2}}\|e^{1}\|_{-1}\|\nabla
e^{1}\|+C\tau^{\frac{3}{2}}\varepsilon^{-\frac{5}{2}}\|u^{1}-u^{0}\|^{2}$
$\displaystyle\leq\,$
$\displaystyle\frac{1}{4}\|e^{1}\|_{-1}^{2}+\frac{\varepsilon\tau}{4}\|\nabla
e^{1}\|^{2}+C\varepsilon^{-(\rho_{1}+3)}\tau^{2}$

Thus, there exists a positive constant $\kappa_{0}^{\prime}$ such that

$\displaystyle\|e^{1}\|_{-1}^{2}+\tau\varepsilon\|\nabla
e^{1}\|^{2}\leq\kappa_{0}^{\prime}\varepsilon^{-\max\\{\rho_{1}+3,\,\rho_{2}+3,\,\rho_{3}+1\\}}\tau^{2}.$
(4.48)

We use the mathematical induction as follows. For $m=1$, (4.41) holds from
(4.48). We suppose (4.41) holds for $m=1,2,\dots,N$, and show that (4.41) is
also valid for $m=N+1$.

Subtracting (2.20) from (3.33), by denoting $e^{n}=u(t_{n})-u^{n}$, we get the
error equation

$\displaystyle\Delta^{-1}\frac{e^{n+1}-e^{n}}{\tau}=\,$
$\displaystyle-\varepsilon\Delta
e^{n+1}+\frac{1}{\varepsilon}\Big{[}f(u(t_{n}))-f(u^{n})\Big{]}$ (4.49)
$\displaystyle+\frac{1}{2\varepsilon}g(u^{n})\sum_{i=1}^{n}D^{2}A(\xi_{i}\,;u^{i}-u^{i-1})$
$\displaystyle-\frac{1}{2\varepsilon}g(u^{n})\big{(}g(u^{n}),u^{n+1}-u^{n}\big{)}+\mathcal{R}^{n+1}.$

Step 2: Estimation of x$\displaystyle\|\nabla
e^{N+1}\|^{2}+\sum_{n=1}^{N}\|\nabla e^{n+1}-\nabla e^{n}\|^{2}$. By taking
the $L^{2}$ inner product of (4.49) with $e^{n+1}-e^{n}\in
L_{0}^{2}({\varOmega})$, we have

$\displaystyle\frac{1}{2}\|\nabla e^{n+1}\|^{2}-\,$
$\displaystyle\frac{1}{2}\|\nabla e^{n}\|^{2}+\frac{1}{2}\|\nabla
e^{n+1}-\nabla
e^{n}\|^{2}+\frac{\tau}{\varepsilon}\Big{\|}\frac{e^{n+1}-e^{n}}{\tau}\Big{\|}_{-1}^{2}$
(4.50) $\displaystyle=\,$
$\displaystyle-\frac{1}{\varepsilon^{2}}\Big{(}f(u(t_{n}))-f(u^{n}),e^{n+1}-e^{n}\Big{)}$
$\displaystyle-\frac{1}{2\varepsilon^{2}}\big{(}g(u^{n}),e^{n+1}-e^{n}\big{)}\sum_{i=1}^{n}D^{2}A(\xi_{i}\,;u^{i}-u^{i-1})$
$\displaystyle+\frac{1}{2\varepsilon^{2}}\big{(}g(u^{n}),e^{n+1}-e^{n}\big{)}(g(u^{n}),u^{n+1}-u^{n})-\frac{1}{\varepsilon}\big{(}\mathcal{R}^{n+1},e^{n+1}-e^{n}\big{)}$
$\displaystyle=:\,$ $\displaystyle K_{1}+K_{2}+K_{3}+K_{4}.$

It follows from (2.6) that

$\displaystyle K_{1}=\,$
$\displaystyle-\frac{1}{\varepsilon^{2}}\Big{(}f(u(t_{n}))-f(u^{n}),e^{n+1}-e^{n}\Big{)}$
(4.51) $\displaystyle\leq\,$
$\displaystyle\varepsilon^{-2}\|f^{\prime}(w^{n})\|_{L^{\infty}}\|e^{n}\|\|e^{n+1}-e^{n}\|$
$\displaystyle\leq\,$ $\displaystyle C\varepsilon^{-2}\|\nabla
e^{n}\|^{\frac{1}{2}}\|e^{n}\|_{-1}^{\frac{1}{2}}\|\nabla e^{n+1}-\nabla
e^{n}\|^{\frac{1}{2}}\|e^{n+1}-e^{n}\|_{-1}^{\frac{1}{2}}$
$\displaystyle\leq\,$ $\displaystyle
C\frac{1}{\gamma_{0}}\varepsilon^{-2}\tau^{\frac{1}{2}}\|\nabla
e^{n}\|\|e^{n}\|_{-1}+2\gamma_{0}\varepsilon^{-2}\tau^{-\frac{1}{2}}\|\nabla
e^{n+1}-\nabla e^{n}\|\|e^{n+1}-e^{n}\|_{-1}$ $\displaystyle\leq\,$
$\displaystyle C\tau\|\nabla
e^{n}\|^{2}+C\varepsilon^{-4}\|e^{n}\|_{-1}^{2}+\gamma_{0}\tau^{-1}\varepsilon^{-4}\|e^{n+1}-e^{n}\|_{-1}^{2}+\gamma_{0}\|\nabla
e^{n+1}-\nabla e^{n}\|^{2}$

where $w_{n}=\theta u^{n}+(1-\theta)u(t_{n})$ with $\theta\in(0,1)$, and
$\gamma_{0}$ is a sufficiently small constant such that the two terms
involving $\|e^{n+1}-e^{n}\|_{-1}^{2}$ and $\|\nabla e^{n+1}-\nabla
e^{n}\|^{2}$ can be absorbed by the left hand side in the Step 4.

We divided the term $K_{2}$ into two parts as

$\displaystyle K_{2}=\,$
$\displaystyle-\frac{1}{2\varepsilon^{2}}\big{(}g(u^{n}),e^{n+1}-e^{n}\big{)}\sum_{i=1}^{n}D^{2}A(\xi_{i}\,;u^{i}-u^{i-1})$
(4.52) $\displaystyle=\,$
$\displaystyle\frac{1}{2\varepsilon^{2}}\big{(}g(u^{n}),e^{n+1}-e^{n}\big{)}\sum_{i=1}^{n}D^{2}A(\xi_{i}\,;e^{i}-e^{i-1})$
$\displaystyle-\frac{1}{2\varepsilon^{2}}\big{(}g(u^{n}),e^{n+1}-e^{n}\big{)}\sum_{i=1}^{n}D^{2}A(\xi_{i}\,;u(t_{i})-u(t_{i-1}))$
$\displaystyle=:\,$ $\displaystyle K_{21}+K_{22}.$

From the assumptions of induction and
$D^{2}A(v;w)\leq\big{(}\|f^{\prime}(v)\|_{L^{\infty}}+\|f(v)\|^{2}\big{)}\|w\|^{2}$,
we have

$\displaystyle\big{(}g(u^{n}),e^{n+1}-e^{n}\big{)}$
$\displaystyle=\frac{1}{\sqrt{\int_{{\varOmega}}F(u^{n})\,{\mathrm{d}}x+c_{0}}}\big{(}f(u^{n}),e^{n+1}-e^{n}\big{)}$
(4.53) $\displaystyle\leq C\|f(u^{n})\|\|e^{n+1}-e^{n}\|$ $\displaystyle\leq
C\|\nabla e^{n+1}-\nabla
e^{n}\|^{\frac{1}{2}}\|e^{n+1}-e^{n}\|_{-1}^{\frac{1}{2}},$
$\displaystyle\sum_{i=1}^{n}D^{2}A(\xi_{i}\,;e^{i}-e^{i-1})$
$\displaystyle\leq\sum_{i=1}^{n}\big{(}\|f^{\prime}(\xi_{i})\|_{L^{\infty}}+\|f(\xi_{i})\|^{2}\big{)}\|e^{i}-e^{i-1}\|^{2}$
(4.54) $\displaystyle\leq
C\Big{(}\sum_{i=1}^{n}\|e^{i}-e^{i-1}\|_{-1}^{2}\Big{)}^{\frac{1}{2}}\Big{(}\sum_{i=1}^{n}\|\nabla
e^{i}-\nabla e^{i-1}\|^{2}\Big{)}^{\frac{1}{2}}$ $\displaystyle\leq
C\kappa_{0}\varepsilon^{-(\alpha_{0}+2)}\tau^{\frac{3}{2}},$
$\displaystyle\sum_{i=1}^{n}D^{2}A(\xi_{i}\,;u(t_{i})-u(t_{i-1}))$
$\displaystyle\leq\sum_{i=1}^{n}\big{(}\|f^{\prime}(\xi_{i})\|_{L^{\infty}}+\|f(\xi_{i})\|^{2}\big{)}\|u(t_{i})-u(t_{i-1})\|^{2}$
(4.55) $\displaystyle\leq
C\tau\sum_{i=1}^{n}\int_{t_{i-1}}^{t_{i}}\|u_{t}\|^{2}\,{\mathrm{d}}s$
$\displaystyle\leq C\tau\int_{0}^{t_{n}}\|u_{t}\|^{2}\,{\mathrm{d}}s$
$\displaystyle\leq C\varepsilon^{-\rho_{2}}\tau.$

Using the estimates above, we have

$\displaystyle K_{21}=\,$
$\displaystyle\frac{1}{2\varepsilon^{2}}\big{(}g(u^{n}),e^{n+1}-e^{n}\big{)}\sum_{i=1}^{n}D^{2}A(\xi_{i}\,;e^{i}-e^{i-1})$
(4.56) $\displaystyle\leq\,$ $\displaystyle C\varepsilon^{-2}\|\nabla
e^{n+1}-\nabla
e^{n}\|^{\frac{1}{2}}\|e^{n+1}-e^{n}\|_{-1}^{\frac{1}{2}}\kappa_{0}\varepsilon^{-(\alpha_{0}+2)}\tau^{\frac{3}{2}}$
$\displaystyle\leq\,$ $\displaystyle
2\gamma_{0}\varepsilon^{-2}\tau^{-\frac{1}{2}}\|\nabla e^{n+1}-\nabla
e^{n}\|\|e^{n+1}-e^{n}\|_{-1}+C\frac{1}{\gamma_{0}}\kappa_{0}^{2}\varepsilon^{-(2\alpha_{0}+6)}\tau^{\frac{7}{2}}$
$\displaystyle\leq\,$
$\displaystyle\gamma_{0}\tau^{-1}\varepsilon^{-4}\|e^{n+1}-e^{n}\|_{-1}^{2}+\gamma_{0}\|\nabla
e^{n+1}-\nabla
e^{n}\|^{2}+C\kappa_{0}^{2}\varepsilon^{-(2\alpha_{0}+6)}\tau^{\frac{7}{2}}.$

and

$\displaystyle K_{22}=\,$
$\displaystyle-\frac{1}{2\varepsilon^{2}}\big{(}g(u^{n}),e^{n+1}-e^{n}\big{)}\sum_{i=1}^{n}D^{2}A(\xi_{i}\,;u(t_{i})-u(t_{i-1}))$
(4.57) $\displaystyle\leq\,$ $\displaystyle
C\varepsilon^{-2}\|e^{n+1}-e^{n}\|_{-1}^{\frac{1}{2}}\|\nabla e^{n+1}-\nabla
e^{n}\|^{\frac{1}{2}}\varepsilon^{-\rho_{2}}\tau$ $\displaystyle\leq\,$
$\displaystyle 2\gamma_{0}\varepsilon^{-2}\tau^{-\frac{1}{2}}\|\nabla
e^{n+1}-\nabla
e^{n}\|\|e^{n+1}-e^{n}\|_{-1}+C\frac{1}{\gamma_{0}}\varepsilon^{-(2\rho_{2}+2)}\tau^{\frac{5}{2}}$
$\displaystyle\leq\,$
$\displaystyle\gamma_{0}\tau^{-1}\varepsilon^{-4}\|e^{n+1}-e^{n}\|_{-1}^{2}+\gamma_{0}\|\nabla
e^{n+1}-\nabla e^{n}\|^{2}+C\varepsilon^{-(2\rho_{2}+2)}\tau^{\frac{5}{2}}.$

In addition, we have

$\displaystyle K_{3}=\,$
$\displaystyle\frac{1}{2\varepsilon^{2}}\big{(}g(u^{n}),e^{n+1}-e^{n}\big{)}\big{(}g(u^{n}),u^{n+1}-u^{n}\big{)}$
(4.58) $\displaystyle=\,$
$\displaystyle-\frac{1}{2\varepsilon^{2}}\big{(}g(u^{n}),e^{n+1}-e^{n}\big{)}\big{(}g(u^{n}),e^{n+1}-e^{n}\big{)}$
$\displaystyle+\frac{1}{2\varepsilon^{2}}\big{(}g(u^{n}),e^{n+1}-e^{n}\big{)}\big{(}g(u^{n}),u(t_{n+1})-u(t_{n})\big{)}$
$\displaystyle=:\,$ $\displaystyle K_{31}+K_{32}.$

where

$\displaystyle K_{31}=\,$
$\displaystyle-\frac{1}{2\varepsilon^{2}}\big{(}g(u^{n}),e^{n+1}-e^{n}\big{)}\big{(}g(u^{n}),e^{n+1}-e^{n}\big{)}$
(4.59) $\displaystyle\leq\,$ $\displaystyle
C\varepsilon^{-2}\|g(u^{n})\|^{2}\|e^{n+1}-e^{n}\|_{-1}\|\nabla e^{n+1}-\nabla
e^{n}\|$ $\displaystyle\leq\,$ $\displaystyle
C\frac{1}{\gamma_{0}}\varepsilon^{-4}\|e^{n+1}-e^{n}\|_{-1}^{2}+\gamma_{0}\|\nabla
e^{n+1}-\nabla e^{n}\|^{2},$

and

$\displaystyle K_{32}=\,$
$\displaystyle\frac{1}{2\varepsilon^{2}}\big{(}g(u^{n}),e^{n+1}-e^{n}\big{)}\big{(}g(u^{n}),u(t_{n+1})-u(t_{n})\big{)}$
(4.60) $\displaystyle\leq\,$ $\displaystyle
C\varepsilon^{-2}\|e^{n+1}-e^{n}\|_{-1}^{\frac{1}{2}}\|\nabla e^{n+1}-\nabla
e^{n}\|^{\frac{1}{2}}\|u(t_{n+1})-u(t_{n})\|$ $\displaystyle\leq\,$
$\displaystyle
2\gamma_{0}\varepsilon^{-2}\tau^{-\frac{1}{2}}\|e^{n+1}-e^{n}\|_{-1}\|\nabla
e^{n+1}-\nabla
e^{n}\|+C\gamma_{0}^{-1}\varepsilon^{-2}\tau^{\frac{3}{2}}\int_{t_{n}}^{t_{n+1}}\|u_{t}\|^{2}\,{\mathrm{d}}s$
$\displaystyle\leq\,$
$\displaystyle\gamma_{0}\tau^{-1}\varepsilon^{-4}\|e^{n+1}-e^{n}\|_{-1}^{2}+\gamma_{0}\|\nabla
e^{n+1}-\nabla
e^{n}\|^{2}+C\varepsilon^{-2}\tau^{\frac{3}{2}}\int_{t_{n}}^{t_{n+1}}\|u_{t}\|^{2}\,{\mathrm{d}}s,$

The last term on the right hand side of (4.50) can be bounded by

$\displaystyle
K_{4}=\frac{1}{\varepsilon}\big{(}\mathcal{R}^{n+1},e^{n+1}-e^{n}\big{)}\leq\gamma_{0}\|\nabla
e^{n+1}-\nabla
e^{n}\|^{2}+\frac{1}{\gamma_{0}}\varepsilon^{-2}\|\mathcal{R}^{n+1}\|_{H^{-1}}^{2}.$
(4.61)

Combining these estimates (4.51)–(4.61) together with (4.50), and taking the
summation for $n=0$ to $N$, we have

$\displaystyle\|\nabla e^{N+1}\|^{2}$
$\displaystyle+(1-6\gamma_{0})\sum_{n=0}^{N}\|\nabla e^{n+1}-\nabla
e^{n}\|^{2}\leq C\tau\sum_{n=0}^{N}\|\nabla
e^{n}\|^{2}+C\varepsilon^{-4}\sum_{n=0}^{N}\|e^{n}\|_{-1}^{2}$ (4.62)
$\displaystyle+(C+4\gamma_{0}\tau^{-1})\varepsilon^{-4}\sum_{n=0}^{N}\|e^{n+1}-e^{n}\|_{-1}^{2}+C\varepsilon^{-2}\sum_{n=0}^{N}\|\mathcal{R}^{n+1}\|_{H^{-1}}^{2}$
$\displaystyle+C\kappa_{0}^{2}\varepsilon^{-(2\alpha_{0}+6)}\tau^{\frac{5}{2}}+C\varepsilon^{-(2\rho_{2}+2)}\tau^{\frac{3}{2}}+C\varepsilon^{-2}\tau^{\frac{3}{2}}\int_{0}^{T}\|u_{t}\|^{2}\,{\mathrm{d}}s.$

Step 3: Estimation of
$\displaystyle\|e^{N+1}\|_{-1}^{2}+\sum_{n=1}^{N}\|e^{n}-e^{n-1}\|_{-1}^{2}$.
Taking the $L^{2}$ inner product of (4.49) with $e^{n+1}\in
L_{0}^{2}({\varOmega})$, we get

$\displaystyle\frac{1}{2\tau}\Big{(}\|e^{n+1}\|_{-1}^{2}$
$\displaystyle-\|e^{n}\|_{-1}^{2}+\|e^{n+1}-e^{n}\|_{-1}^{2}\Big{)}+\varepsilon\|\nabla
e^{n+1}\|^{2}+\frac{1}{\varepsilon}\Big{(}f(u(t_{n}))-f(u^{n}),e^{n+1}\Big{)}$
$\displaystyle+\frac{1}{2\varepsilon}\big{(}g(u^{n}),e^{n+1}\big{)}\sum_{i=1}^{n}D^{2}A(\xi_{i}\,;u^{i}-u^{i-1})$
(4.63)
$\displaystyle-\frac{1}{2\varepsilon}\big{(}g(u^{n}),e^{n+1}\big{)}\big{(}g(u^{n}),u^{n+1}-u^{n}\big{)}+\big{(}\mathcal{R}^{n+1},e^{n+1}\big{)}=0.$

We denote $w^{n}=\theta u^{n}+(1-\theta)u(t_{n})$ with $\theta\in(0,1)$, and
using the Taylor expansion to get

$\displaystyle\frac{1}{\varepsilon}\Big{(}f(u(t_{n}))-f(u^{n}),e^{n+1}\Big{)}=\,$
$\displaystyle\frac{1}{\varepsilon}\Big{(}f^{\prime}(u(t_{n}))e^{n}-\frac{1}{2}f^{\prime\prime}(w^{n})(e^{n})^{2},e^{n+1}\Big{)}$
$\displaystyle=\,$
$\displaystyle\frac{1}{\varepsilon}\Big{(}f^{\prime}(u(t_{n}))e^{n+1},e^{n+1}\Big{)}-\frac{1}{\varepsilon}\Big{(}f^{\prime}(u(t_{n}))(e^{n+1}-e^{n}),e^{n+1}\Big{)}$
$\displaystyle-\frac{1}{2\varepsilon}\Big{(}f^{\prime\prime}(w^{n})(e^{n})^{2},e^{n+1}\Big{)}.$

Then, (4) becomes

$\displaystyle\frac{1}{2\tau}\Big{(}\|e^{n+1}\|_{-1}^{2}-\,$
$\displaystyle\|e^{n}\|_{-1}^{2}+\|e^{n+1}-e^{n}\|_{-1}^{2}\Big{)}+\varepsilon\|\nabla
e^{n+1}\|^{2}+\frac{1}{\varepsilon}\Big{(}f^{\prime}(u(t_{n}))e^{n+1},e^{n+1}\Big{)}$
(4.64) $\displaystyle=\,$
$\displaystyle\frac{1}{\varepsilon}\Big{(}f^{\prime}(u(t_{n}))(e^{n+1}-e^{n}),e^{n+1}\Big{)}+\frac{1}{2\varepsilon}\Big{(}f^{\prime\prime}(w^{n})(e^{n})^{2},e^{n+1}\Big{)}$
$\displaystyle-\frac{1}{2\varepsilon}\big{(}g(u^{n}),e^{n+1}\big{)}\sum_{i=1}^{n}D^{2}A(\xi_{i}\,;u^{i}-u^{i-1})$
$\displaystyle+\frac{1}{2\varepsilon}\big{(}g(u^{n}),e^{n+1}\big{)}\big{(}g(u^{n}),u^{n+1}-u^{n}\big{)}-\big{(}\mathcal{R}^{n+1},e^{n+1}\big{)}$
$\displaystyle=:\,$ $\displaystyle T_{1}+T_{2}+T_{3}+T_{4}+T_{5}.$

Using the Lemma 4.1, there holds

$\displaystyle\varepsilon\|\nabla
e^{n+1}\|^{2}+\frac{1}{\varepsilon}\Big{(}f^{\prime}(u(t_{n}))e^{n+1},e^{n+1}\Big{)}\geq-\lambda_{0}\|e^{n+1}\|_{-1}^{2}.$
(4.65)

If the entire $\varepsilon\|\nabla e^{n+1}\|^{2}$ term is used to control the
the term $\varepsilon^{-1}\big{(}f^{\prime}(u(t_{n}))e^{n+1},e^{n+1}\big{)}$,
we will not be able to control the $\|\nabla e^{n+1}\|^{2}$ terms in $T_{j}$,
$j=1,\dots,5$. So we apply (4.65) with a scaling factor $(1-\eta)$ close to
but smaller than 1, to get

$\displaystyle-(1-\eta)\frac{1}{\varepsilon}\Big{(}f^{\prime}(u(t_{n}))e^{n+1},e^{n+1}\Big{)}\leq(1-\eta)\lambda_{0}\|e^{n+1}\|_{-1}^{2}+(1-\eta)\varepsilon\|\nabla
e^{n+1}\|^{2}.$ (4.66)

On the other hand,

$\displaystyle-\frac{\eta}{\varepsilon}\Big{(}f^{\prime}(u(t_{n}))e^{n+1},e^{n+1}\Big{)}\leq\frac{C\eta}{\varepsilon}\|e^{n+1}\|^{2}\leq\frac{C\eta}{\varepsilon^{2}\eta_{1}}\|e^{n+1}\|_{-1}^{2}+\frac{\eta\eta_{1}}{4}\|\nabla
e^{n+1}\|^{2}.$ (4.67)

The first term $T_{1}$ on the right-hand side of (4.64) can be bounded by

$\displaystyle T_{1}=\,$
$\displaystyle\frac{1}{\varepsilon}\Big{(}f^{\prime}(u(t_{n}))(e^{n+1}-e^{n}),e^{n+1}\Big{)}$
(4.68) $\displaystyle\leq\,$
$\displaystyle\varepsilon^{-1}\|f^{\prime}(u(t_{n}))\|_{L^{\infty}}\|e^{n+1}-e^{n}\|\|e^{n+1}\|$
$\displaystyle\leq\,$ $\displaystyle
C\varepsilon^{-1}\|e^{n+1}-e^{n}\|_{-1}^{\frac{1}{2}}\|\nabla e^{n+1}-\nabla
e^{n}\|^{\frac{1}{2}}\|e^{n+1}\|_{-1}^{\frac{1}{2}}\|\nabla
e^{n+1}\|^{\frac{1}{2}}$ $\displaystyle\leq\,$
$\displaystyle\tau^{-\frac{1}{2}}\|e^{n+1}-e^{n}\|_{-1}\|e^{n+1}\|_{-1}+C\varepsilon^{-2}\tau^{\frac{1}{2}}\|\nabla
e^{n+1}-\nabla e^{n}\|\|\nabla e^{n+1}\|$ $\displaystyle\leq\,$
$\displaystyle\frac{1}{2}\gamma_{0}\tau^{-1}\|e^{n+1}-e^{n}\|_{-1}^{2}+C\frac{1}{\gamma_{0}}\|e^{n+1}\|_{-1}^{2}+\frac{1}{16}\varepsilon^{\eta_{2}}\|\nabla
e^{n+1}\|^{2}+T_{1}^{*}.$

where $\displaystyle T_{1}^{*}:=C\varepsilon^{-(\eta_{2}+4)}\tau\|\nabla
e^{n+1}-\nabla e^{n}\|^{2}$ and $\gamma_{0}$ is sufficiently small. To control
the last term $T_{1}^{*}$ on the right-hand side of (4.68), we assume that
$\tau\leq\tilde{C}_{1}\varepsilon^{\eta_{2}+8}$ to get

$\displaystyle\tau\sum_{n=0}^{N}T_{1}^{*}:=\,$ $\displaystyle
C\varepsilon^{-(\eta_{2}+4)}\tau^{2}\sum_{n=0}^{N}\|\nabla e^{n+1}-\nabla
e^{n}\|^{2}$ $\displaystyle=\,$ $\displaystyle
C\varepsilon^{-(\eta_{2}+8)}\tau\Big{(}\varepsilon^{4}\tau\sum_{n=0}^{N}\|\nabla
e^{n+1}-\nabla e^{n}\|^{2}\Big{)}$ $\displaystyle\leq\,$ $\displaystyle
frac{1}{16}\varepsilon^{4}\tau\sum_{n=0}^{N}\|\nabla e^{n+1}-\nabla
e^{n}\|^{2}.$

Using the Sobolev interpolation inequality, we have for $v\in
L_{0}^{2}({\varOmega})\cap H^{1}({\varOmega})$

$\displaystyle\|v\|_{L^{4}}\leq C\|\nabla
v\|^{\frac{d}{4}}\|v\|^{1-\frac{d}{4}}\leq C\|\nabla
v\|^{\frac{d}{4}}\|v\|_{-1}^{\frac{1}{2}-\frac{d}{8}}\|\nabla
v\|^{\frac{1}{2}-\frac{d}{8}}=C\|\nabla
v\|^{\frac{1}{2}+\frac{d}{8}}\|v\|_{-1}^{\frac{1}{2}-\frac{d}{8}},$ (4.69)

which together with the assumptions of induction yields

$\displaystyle T_{2}=\,$
$\displaystyle\frac{1}{2\varepsilon}\Big{(}f^{\prime\prime}(w^{n})(e^{n})^{2},e^{n+1}\Big{)}$
(4.70) $\displaystyle\leq\,$ $\displaystyle
C\varepsilon^{-1}\|e^{n}\|_{L^{4}}^{2}\|e^{n+1}\|$ $\displaystyle\leq\,$
$\displaystyle C\varepsilon^{-1}\|\nabla
e^{n}\|^{1+\frac{d}{4}}\|e^{n}\|_{-1}^{1-\frac{d}{4}}\|e^{n+1}\|_{-1}^{\frac{1}{2}}\|\nabla
e^{n+1}\|^{\frac{1}{2}}$ $\displaystyle\leq\,$ $\displaystyle
C\kappa_{0}^{2}\varepsilon^{-(2\alpha_{0}+5+d)}\tau^{3-\frac{d}{4}}\|e^{n+1}\|_{-1}^{\frac{1}{2}}\|\nabla
e^{n+1}\|^{\frac{1}{2}}$ $\displaystyle\leq\,$ $\displaystyle
C\varepsilon^{\frac{\eta_{2}}{2}}\|\nabla
e^{n+1}\|\|e^{n+1}\|_{-1}+C\kappa_{0}^{4}\varepsilon^{-(4\alpha_{0}+10+2d+\frac{\eta_{2}}{2})}\tau^{6-\frac{d}{2}}$
$\displaystyle\leq\,$ $\displaystyle
C\|e^{n+1}\|_{-1}^{2}+\frac{1}{16}\varepsilon^{\eta_{2}}\|\nabla
e^{n+1}\|^{2}+C\kappa_{0}^{4}\varepsilon^{-(4\alpha_{0}+10+2d+\frac{\eta_{2}}{2})}\tau^{6-\frac{d}{2}}.$

It follows from (4.54)–(4.55) that

$\displaystyle T_{3}=\,$
$\displaystyle-\frac{1}{2\varepsilon}\big{(}g(u^{n}),e^{n+1}\big{)}\sum_{i=1}^{n}D^{2}A(\xi_{i}\,;u^{i}-u^{i-1})$
(4.71) $\displaystyle\leq\,$
$\displaystyle\frac{1}{\varepsilon}\big{(}g(u^{n}),e^{n+1}\big{)}\sum_{i=1}^{n}D^{2}A(\xi_{i}\,;e^{i}-e^{i-1})+\frac{1}{\varepsilon}\big{(}g(u^{n}),e^{n+1}\big{)}\sum_{i=1}^{n}D^{2}A(\xi_{i}\,;u(t_{i})-u(t_{i-1}))$
$\displaystyle\leq\,$ $\displaystyle
C\varepsilon^{-1}\|g(u^{n})\|\|e^{n+1}\|\Big{(}\kappa_{0}\varepsilon^{-(\alpha_{0}+2)}\tau^{\frac{3}{2}}+\varepsilon^{-\rho_{2}}\tau\Big{)}$
$\displaystyle\leq\,$ $\displaystyle
C\varepsilon^{-1}\|e^{n+1}\|_{-1}^{\frac{1}{2}}\|\nabla
e^{n+1}\|^{\frac{1}{2}}\Big{(}\kappa_{0}\varepsilon^{-(\alpha_{0}+2)}\tau^{\frac{3}{2}}+\varepsilon^{-\rho_{2}}\tau\Big{)}$
$\displaystyle\leq\,$
$\displaystyle\varepsilon^{\frac{\eta_{2}}{2}}\|e^{n+1}\|_{-1}\|\nabla
e^{n+1}\|+C\kappa_{0}^{2}\varepsilon^{-(2\alpha_{0}+6+\frac{\eta_{2}}{2})}\tau^{3}+C\varepsilon^{-(2\rho_{2}+2+\frac{\eta_{2}}{2})}\tau^{2}$
$\displaystyle\leq\,$ $\displaystyle
C\|e^{n+1}\|_{-1}^{2}+\frac{1}{16}\varepsilon^{\eta_{2}}\|\nabla
e^{n+1}\|^{2}+C\kappa_{0}^{2}\varepsilon^{-(2\alpha_{0}+6+\frac{\eta_{2}}{2})}\tau^{3}+C\varepsilon^{-(2\rho_{2}+2+\frac{\eta_{2}}{2})}\tau^{2}.$

In order to estimate the term $T_{4}$, we divided it into two parts as

$\displaystyle T_{4}=\,$
$\displaystyle\frac{1}{2\varepsilon}\big{(}g(u^{n}),e^{n+1}\big{)}\big{(}g(u^{n}),u^{n+1}-u^{n}\big{)}$
(4.72) $\displaystyle=\,$
$\displaystyle-\frac{1}{2\varepsilon}\big{(}g(u^{n}),e^{n+1}\big{)}\big{(}g(u^{n}),e^{n+1}-e^{n}\big{)}$
$\displaystyle+\frac{1}{2\varepsilon}\big{(}g(u^{n}),e^{n+1}\big{)}\big{(}g(u^{n}),u(t_{n+1})-u(t_{n})\big{)}$
$\displaystyle=:\,$ $\displaystyle T_{41}+T_{42}.$

Similarly to the estimate of $T_{1}$, the term $T_{41}$ can be bounded by

$\displaystyle T_{41}=\,$
$\displaystyle-\frac{1}{2\varepsilon}\big{(}g(u^{n}),e^{n+1}\big{)}\big{(}g(u^{n}),e^{n+1}-e^{n}\big{)}$
(4.73) $\displaystyle\leq\,$ $\displaystyle
C\varepsilon^{-1}\|g(u^{n})\|^{2}\|e^{n+1}\|\|e^{n+1}-e^{n}\|$
$\displaystyle\leq\,$ $\displaystyle
C\varepsilon^{-1}\|e^{n+1}-e^{n}\|_{-1}^{\frac{1}{2}}\|\nabla e^{n+1}-\nabla
e^{n}\|^{\frac{1}{2}}\|e^{n+1}\|_{-1}^{\frac{1}{2}}\|\nabla
e^{n+1}\|^{\frac{1}{2}}$ $\displaystyle\leq\,$
$\displaystyle\tau^{-\frac{1}{2}}\|e^{n+1}-e^{n}\|_{-1}\|e^{n+1}\|_{-1}+C\varepsilon^{-2}\tau^{\frac{1}{2}}\|\nabla
e^{n+1}-\nabla e^{n}\|\|\nabla e^{n+1}\|$ $\displaystyle\leq\,$
$\displaystyle\frac{1}{2}\gamma_{0}\tau^{-1}\|e^{n+1}-e^{n}\|_{-1}^{2}+C\frac{1}{\gamma_{0}}\|e^{n+1}\|_{-1}^{2}+\frac{1}{16}\varepsilon^{\eta_{2}}\|\nabla
e^{n+1}\|^{2}+T_{1}^{*},$

which together with the condition
$\tau\leq\tilde{C}_{1}\varepsilon^{\eta_{2}+8}$ yields

$\displaystyle\tau\sum_{n=0}^{N}T_{1}^{*}=C\varepsilon^{-(\eta_{2}+4)}\tau^{2}\sum_{n=0}^{N}\|\nabla
e^{n+1}-\nabla
e^{n}\|^{2}\leq\frac{1}{16}\varepsilon^{4}\tau\sum_{n=0}^{N}\|\nabla
e^{n+1}-\nabla e^{n}\|^{2}.$

In addition, we obtain

$\displaystyle T_{42}=\,$
$\displaystyle\frac{1}{2\varepsilon}\big{(}g(u^{n}),e^{n+1}\big{)}\big{(}g(u^{n}),u(t_{n+1})-u(t_{n})\big{)}$
(4.74) $\displaystyle\leq\,$ $\displaystyle C\varepsilon^{-1}\|\nabla
e^{n+1}\|^{\frac{1}{2}}\|e^{n+1}\|_{-1}^{\frac{1}{2}}\|u(t_{n+1})-u(t_{n})\|$
$\displaystyle\leq\,$ $\displaystyle\varepsilon^{\frac{\eta_{2}}{2}}\|\nabla
e^{n+1}\|\|e^{n+1}\|_{-1}+C\varepsilon^{-(2+\frac{\eta_{2}}{2})}\tau\int_{t_{n}}^{t_{n+1}}\|u_{t}\|^{2}\,{\mathrm{d}}s$
$\displaystyle\leq\,$ $\displaystyle
C\|e^{n+1}\|_{-1}^{2}+\frac{1}{16}\varepsilon^{\eta_{2}}\|\nabla
e^{n+1}\|^{2}+C\varepsilon^{-(2+\frac{\eta_{2}}{2})}\tau\int_{t_{n}}^{t_{n+1}}\|u_{t}\|^{2}\,{\mathrm{d}}s.$

For $T_{5}$, using Cauchy-Schwartz inequality and Poincaré’s inequality for
$e^{n+1}\in L_{0}^{2}({\varOmega})$, we have

$\displaystyle
T_{5}=-\big{(}\mathcal{R}^{n+1},e^{n+1}\big{)}\leq\|\mathcal{R}^{n+1}\|_{H^{-1}}\|e^{n+1}\|_{H^{1}}\leq
8\varepsilon^{-\eta_{2}}\|\mathcal{R}^{n+1}\|_{H^{-1}}^{2}+\frac{1}{16}\varepsilon^{\eta_{2}}\|\nabla
e^{n+1}\|^{2}.$ (4.75)

Combining the estimates (4.68)–(4.75) together with (4.64), taking the
summation for $n=0$ to $N$, and assuming that
$\tau\leq\tilde{C}_{1}\varepsilon^{\eta_{2}+8}$, we have

$\displaystyle\frac{1}{2}\|e^{N+1}\|_{-1}^{2}$
$\displaystyle+\frac{1}{2}(1-2\gamma_{0})\sum_{n=0}^{N}\|e^{n+1}-e^{n}\|_{-1}^{2}+\varepsilon\tau\sum_{n=0}^{N}\|\nabla
e^{n+1}\|^{2}$ (4.76) $\displaystyle\leq\,$
$\displaystyle\Big{[}C+(1-\eta)\lambda_{0}+\frac{C\eta}{\varepsilon^{2}\eta_{1}}\Big{]}\tau\sum_{n=0}^{N}\|e^{n+1}\|_{-1}^{2}+8\varepsilon^{-\eta_{2}}\tau\sum_{n=0}^{N}\|\mathcal{R}^{n+1}\|_{H^{-1}}^{2}$
$\displaystyle+\Big{[}(1-\eta)\varepsilon+\frac{\eta\eta_{1}}{4}+\frac{\varepsilon^{\eta_{2}}}{2}\Big{]}\tau\sum_{n=0}^{N}\|\nabla
e^{n+1}\|^{2}+\frac{1}{8}\varepsilon^{4}\tau\sum_{n=0}^{N}\|\nabla
e^{n+1}-\nabla e^{n}\|^{2}$
$\displaystyle+C\kappa_{0}^{4}\varepsilon^{-(4\alpha_{0}+10+2d+\frac{\eta_{2}}{2})}\tau^{6-\frac{d}{2}}+C\kappa_{0}^{2}\varepsilon^{-(2\alpha_{0}+6+\frac{\eta_{2}}{2})}\tau^{3}$
$\displaystyle+C\varepsilon^{-(2\rho_{2}+2+\frac{\eta_{2}}{2})}\tau^{2}+C\varepsilon^{-(2+\frac{\eta_{2}}{2})}\tau^{2}\int_{0}^{T}\|u_{t}\|^{2}\,{\mathrm{d}}s.$

By taking $\eta=\varepsilon^{3}$, $\eta_{1}=\varepsilon$ and $\eta_{2}=4$, we
have

$\displaystyle(1-\eta)\varepsilon+\frac{\eta\eta_{1}}{4}+\frac{\varepsilon^{\eta_{2}}}{2}=\varepsilon-\frac{1}{4}\varepsilon^{4},$
(4.77)

which together with (3.22c) gives

$\displaystyle\|e^{N+1}\|_{-1}^{2}$
$\displaystyle+(1-2\gamma_{0})\sum_{n=0}^{N}\|e^{n+1}-e^{n}\|_{-1}^{2}+\frac{1}{2}\varepsilon^{4}\tau\sum_{n=0}^{N}\|\nabla
e^{n+1}\|^{2}$ (4.78) $\displaystyle\leq\,$ $\displaystyle
C\tau\sum_{n=0}^{N}\|e^{n+1}\|_{-1}^{2}+\frac{1}{4}\tau\varepsilon^{4}\sum_{n=0}^{N}\|\nabla
e^{n+1}-\nabla
e^{n}\|^{2}+C\varepsilon^{-4}\tau\sum_{n=0}^{N}\|\mathcal{R}^{n+1}\|_{H^{-1}}^{2}$
$\displaystyle+C\kappa_{0}^{4}\varepsilon^{-(4\alpha_{0}+12+2d)}\tau^{6-\frac{d}{2}}+C\kappa_{0}^{2}\varepsilon^{-(2\alpha_{0}+8)}\tau^{3}$
$\displaystyle+C\varepsilon^{-(2\rho_{2}+4)}\tau^{2}+C\varepsilon^{-(\rho_{2}+4)}\tau^{2}.$

Step 4: Completion of the proof. We now conclude the proof by the following
induction argument which is based on the results from Step 1 to Step 3. By
multiplying $\tau\varepsilon^{4}$ on both sides of (4.62), combining the
estimate (4.78), and together with Lemma 3.1, we obtain

$\displaystyle\|e^{N+1}\|_{-1}^{2}$ $\displaystyle+\tau\varepsilon^{4}\|\nabla
e^{N+1}\|^{2}+(1-C\tau-6\gamma_{0})\sum_{n=0}^{N}\|e^{n+1}-e^{n}\|_{-1}^{2}$
(4.79)
$\displaystyle+(\frac{3}{4}-6\gamma_{0})\tau\varepsilon^{4}\sum_{n=0}^{N}\|\nabla
e^{n+1}-\nabla
e^{n}\|^{2}+\frac{1}{2}\varepsilon^{4}\tau\sum_{n=0}^{N}\|\nabla
e^{n+1}\|^{2}$ $\displaystyle\leq\,$ $\displaystyle
C_{0}\tau\sum_{n=0}^{N}\Big{(}\|e^{n}\|_{-1}^{2}+\tau\varepsilon^{4}\|\nabla
e^{n}\|^{2}\Big{)}+C_{1}\kappa_{0}^{4}\varepsilon^{-(4\alpha_{0}+12+2d)}\tau^{6-\frac{d}{2}}$
$\displaystyle+C_{2}\kappa_{0}^{2}\varepsilon^{-(2\alpha_{0}+8)}\tau^{3}+C_{3}\kappa_{0}^{2}\varepsilon^{-2\alpha_{0}}\tau^{\frac{7}{2}}+C_{4}\varepsilon^{-\max\\{\rho_{1}+3,\,2\rho_{2}+4,\,\rho_{2}+6,\,\rho_{3}+4\\}}\tau^{2}.$

in which the term
$\kappa_{0}^{\prime}\varepsilon^{-\max\\{\rho_{1}+3,\,\rho_{2}+3,\,\rho_{3}+1\\}}\tau^{2}$
is absorbed in
$C_{4}\varepsilon^{-\max\\{\rho_{1}+3,\,2\rho_{2}+4,\,\rho_{2}+6,\,\rho_{3}+4\\}}\tau^{2}$.

Suppose that for sufficiently small constant $\gamma_{0}$ satisfying
$\frac{3}{4}-6\gamma_{0}\geq\frac{1}{2}$ and sufficiently small $\tau$
satisfying

$\displaystyle\tau\leq\Big{(}\frac{C_{4}}{C_{1}}\kappa_{0}^{-4}\Big{)}^{\frac{2}{8-d}}\varepsilon^{\frac{6\alpha_{0}+24+4d}{8-d}},\quad\tau\leq\Big{(}\frac{C_{4}}{C_{2}}\kappa_{0}^{-2}\Big{)}\varepsilon^{\alpha_{0}+8},\quad\tau\leq\Big{(}\frac{C_{4}}{C_{3}}\kappa_{0}^{-2}\Big{)}^{\frac{2}{3}}\varepsilon^{\frac{(2\alpha_{0}+4)}{3}},$
(4.80)

then, by denoting
$\alpha_{0}:=\max\\{\rho_{1}+3,\,2\rho_{2}+4,\,\rho_{2}+6,\,\rho_{3}+4\\}$, we
derive

$\displaystyle\|e^{N+1}\|_{-1}^{2}$ $\displaystyle+\tau\varepsilon^{4}\|\nabla
e^{N+1}\|^{2}+\frac{1}{2}\sum_{n=0}^{N}\|e^{n+1}-e^{n}\|_{-1}^{2}+\frac{1}{2}\varepsilon^{4}\sum_{n=0}^{N}\tau\|\nabla
e^{n+1}\|^{2}$ (4.81)
$\displaystyle+\frac{1}{2}\tau\varepsilon^{4}\sum_{n=0}^{N}\|\nabla
e^{n+1}-\nabla
e^{n}\|^{2}\leq\,C_{0}\tau\sum_{n=0}^{N}\Big{(}\|e^{n}\|_{-1}^{2}+\tau\varepsilon^{4}\|\nabla
e^{n}\|^{2}\Big{)}+4C_{4}\varepsilon^{-\alpha_{0}}\tau^{2}.$

We denote $\kappa_{0}:=4C_{4}e^{(C_{0}T)}$ and use the Gronwall’s inequality
to get

$\displaystyle\|e^{N+1}\|_{-1}^{2}$
$\displaystyle+\frac{1}{2}\sum_{n=0}^{N}\|e^{n+1}-e^{n}\|_{-1}^{2}+\frac{1}{2}\varepsilon^{4}\sum_{n=0}^{N}\tau\|\nabla
e^{n+1}\|^{2}+\tau\varepsilon^{4}\|\nabla e^{N+1}\|^{2}$ (4.82)
$\displaystyle+\frac{1}{2}\tau\varepsilon^{4}\sum_{n=0}^{m}\|\nabla
e^{n+1}-\nabla
e^{n}\|^{2}\leq\,4C_{4}e^{(C_{0}T)}\varepsilon^{-\alpha_{0}}\tau^{2}=\kappa_{0}\varepsilon^{-\alpha_{0}}\tau^{2}.$

The induction is completed.

In the above proof, we have used these conditions:

$\displaystyle\tau\leq\tilde{C}_{1}\varepsilon^{12},\quad\tau\leq\tilde{C}_{2}\varepsilon^{\frac{6\alpha_{0}+24+4d}{8-d}},\quad\tau\leq\tilde{C}_{3}\varepsilon^{\alpha_{0}+8},\quad\tau\leq\tilde{C}_{4}\varepsilon^{\frac{(2\alpha_{0}+4)}{3}},$
(4.83)

where
$\tilde{C}_{2}=\Big{(}\frac{C_{4}}{C_{1}}\kappa_{0}^{-4}\Big{)}^{\frac{2}{8-d}}$,
$\tilde{C}_{3}=\Big{(}\frac{C_{4}}{C_{2}}\kappa_{0}^{-2}\Big{)}$, and
$\tilde{C}_{4}=\Big{(}\frac{C_{4}}{C_{3}}\kappa_{0}^{-2}\Big{)}^{\frac{2}{3}}$.
By denoting

$\beta_{0}=\left\\{\begin{aligned} &\alpha_{0}+8\quad&&\mbox{for}\quad d=2\\\
&\frac{6\alpha_{0}+36}{5}\quad&&\mbox{for}\quad
d=3\end{aligned}\right.\qquad\mbox{and}\qquad\tilde{C}=\min\\{\tilde{C}_{1},\tilde{C}_{2},\tilde{C}_{3},\tilde{C}_{4}\\},$
(4.84)

we specify the final condition on $\tau$, that is,
$\tau\leq\tilde{C}\varepsilon^{\beta_{0}}$.

$\square$

##  5 Numerical experiments

In this section, we present a two-dimensional numerical test to validate the
theoretical results on the energy decay properties proved in Theorem 3.1, as
well as the convergence rates of the proposed method given in Theorem 4.1. All
the computations are performed using the software package NGSolve
(https://ngsolve.org).

We solve the Cahn-Hilliard equation (1.1) on the two-dimensional square
${\varOmega}=[0,1]\times[0,1]$ under Neumann boundary conditions by using the
proposed scheme (2.11) with the following initial condition

$\displaystyle u_{0}(x,y)=$
$\displaystyle\tanh\big{(}((x-0.65)^{2}+(y-0.5)^{2}-0.1^{2})/\varepsilon\big{)}$
(5.85)
$\displaystyle\times\tanh\big{(}((x-0.35)^{2}+(y-0.5)^{2}-0.125^{2})/\varepsilon\big{)},$

where $\tanh(x):=(e^{x}-e^{-x})/(e^{x}+e^{-x})$. This type of initial
condition is also adopted in [21, 25], where the set of the zero-level of the
initial function $u_{0}(x,y)$ encloses two circles of radius $0.1$ and
$0.125$, respectively.

To obtain a $C^{4}$ potential function $F(v)$ that satisfies the assumption
3.1, we modify the common double-well potential
$F(v)=\frac{1}{4}(v^{2}-1)^{2}$ by setting $M=2$ in (3.23) to get a cut-off
function $\hat{F}(v)\in C^{4}(\mathbb{R})$. Correspondingly, the ninth-order
polynomials $\Phi_{+}(v)$ and $\Phi_{-}(v)$ in (3.23) are determined with the
following conditions

$\left\\{\begin{aligned}
&\Phi_{+}^{(i)}(M)=F^{(i)}(M)\quad\mbox{and}\quad\Phi_{-}^{(i)}(-M)=F^{(i)}(-M)\qquad\mbox{for}\,\,\,i=0,1,2,3,4,\\\
&\Phi_{+}(2M)=\Phi_{-}(-2M)=\frac{1}{4}((2M)^{2}-1)^{2},\\\
&\Phi_{+}^{(1)}(2M)=\Phi_{-}^{(1)}(-2M)=(2M)^{3}-2M,\\\
&\Phi_{+}^{(i)}(2M)=\Phi_{-}^{(i)}(-2M)=0\qquad\mbox{for}\,\,\,i=2,3,4.\end{aligned}\right.$
(5.86)

Note that the truncation point $M=2$ used here are for convenience only. For
simplicity, we still denote the modified function $\hat{F}(u)$ by $F(u)$.

The spatial discretization is done by using the Galerkin finite element
method. Let $S_{h}$ denotes the $P_{s}$ conforming finite element space
defined by

$\displaystyle S_{h}:=\\{v_{h}\in C(\bar{\varOmega});\,\,v_{h}|_{K}\in
P_{s}(K),\,\,\,\forall\,K\in\mathcal{T}_{h}\\},$

where $\mathcal{T}_{h}$ is a quasi-uniform triangulation of ${\varOmega}$. We
introduce space notation $S_{h}^{\circ}:=\\{v_{h}\in
S_{h};\,\,(v_{h},1)=0\\}$, and define the discrete inverse Laplace operator
$-\Delta_{h}^{-1}:L_{0}^{2}({\varOmega})\to S_{h}^{\circ}$ such that

$\displaystyle\big{(}\nabla(-\Delta_{h}^{-1})v,\nabla\eta_{h}\big{)}=(v,\eta_{h})\qquad\forall\,\,\,\eta_{h}\in
S_{h}.$ (5.87)

Since the exact solution of the considered problem is not known, we compute
the orders of convergence by the formula

$\mbox{order of
convergence}=\log\Bigg{(}\frac{\|u_{N}^{(\tau)}-u_{N}^{(\tau/2)}\|_{h,-1}}{\|u_{N}^{(\tau/2)}-u_{N}^{(\tau/4)}\|_{h,-1}}\Bigg{)}/\log(2)$

based on the finest three meshes, where $u_{N}^{(\tau)}$ denotes the numerical
solution at $t_{N}=T$ computed by using a stepsize $\tau$, and
$\|v\|_{h,-1}:=\sqrt{(v,-\Delta_{h}^{-1}v)}$ for $v\in
L_{0}^{2}({\varOmega})$.

The time discretization errors in $\|\cdot\|_{h,-1}$-norm are presented in
Figure 1 (left) for four different $\varepsilon=0.08,0.06,0.05,0.04$ at
$T=0.005$, where we have used finite elements of degree $s=3$ with a
sufficiently spatial mesh $h=1/64$ so that the error from spatial
discretization is negligibly small in observing the temporal convergence
rates. From Figure 1 (left), we see that the error of time discretization is
$O(\tau)$, which is consistent with the theoretical results proved in Theorem
4.1. In addition, Figure 1 (right) shows the evolution in time of the discrete
SAV energy for four different $\varepsilon$, which should be decreasing
according to Theorem 3.1. This graph clearly confirms this decay property.
Therefore, the numerical experiments are in accordance with our theoretical
results.

Figure 1: (left) Time discretization errors; (right) Evolution of the SAV
energy.

Figure 2 shows snapshots of the numerical interface for four different
$\varepsilon=0.08,0.06,0.05,0.04$ at six fixed time points. They clearly
indicate that at each time point, as $\varepsilon$ tends to zero, the
numerical interface converges to the sharp interface of the Hele-haw flow,
which is consistent with the phenomenon stated in [21, 25]. It also shows that
for larger $\varepsilon$, the numerical interface evolves faster in time.

Figure 2: Snapshots of the zero-level sets of the numerical solutions.

Acknowledgements. The work of Shu Ma was partially supported by the Research
Grants Council of the Hong Kong Special Administrative Region, China. (Project
Nos. CityU 11302718, CityU 11300621). The work of Weifeng Qiu was partially
supported by the Research Grants Council of the Hong Kong Special
Administrative Region, China. (Project Nos. CityU 11302718, CityU 11300621).

## References

  * [1] G. Akrivis and B. Li.  Error estimates for fully discrete bdf finite element approximations of the Allen-Cahn equation.  IMA Journal of Numerical Analysis, 42(1):363–391, 2022.
  * [2] G. Akrivis, B. Li, and D. Li.  Energy-decaying extrapolated RK–SAV methods for the Allen-Cahn and Cahn-Hilliard equations.  SIAM Journal on Scientific Computing, 41(6):A3703–A3727, 2019.
  * [3] N. D. Alikakos and G. Fusco.  The spectrum of the Cahn-Hilliard operator for generic interface in higher space dimensions.  Indiana University mathematics journal, 42(2):637–674, 1993.
  * [4] J. W. Barrett and J. F. Blowey.  An error bound for the finite element approximation of the Cahn-Hilliard equation with logarithmic free energy.  Numerische Mathematik, 72(1):1–20, 1995.
  * [5] S. Bartels.  Robust a priori error analysis for the approximation of degree-one Ginzburg-Landau vortices.  ESAIM: Mathematical Modelling and Numerical Analysis, 39(5):863–882, 2005.
  * [6] S. Bartels and R. Müller.  Error control for the approximation of Allen-Cahn and Cahn-Hilliard equations with a logarithmic potential.  Numerische Mathematik, 119(3):409–435, 2011.
  * [7] S. Bartels and R. Müller.  Quasi-optimal and robust a posteriori error estimates in ${L}^{\infty}({L}^{2})$ for the approximation of Allen-Cahn equations past singularities.  Mathematics of computation, 80(274):761–780, 2011.
  * [8] S. Bartels, R. Müller, and C. Ortner.  Robust a priori and a posteriori error analysis for the approximation of Allen-Cahn and Ginzburg-Landau equations past topological changes.  SIAM Journal on Numerical Analysis, 49(1):110–134, 2011.
  * [9] L. A. Caffarelli and N. E. Muler.  An ${L}^{\infty}$ bound for solutions of the Cahn-Hilliard equation.  Archive for Rational Mechanics and Analysis, 133(2):129–144, 1995\.
  * [10] J. W. Cahn and J. E. Hilliard.  Free energy of a nonuniform system. I. interfacial free energy.  The Journal of chemical physics, 28(2):258–267, 1958.
  * [11] Y. Cai, H. Choi, and J. Shen.  Error estimates for time discretizations of Cahn-Hilliard and Allen-Cahn phase-field models for two-phase incompressible flows.  Numerische Mathematik, 137(2):417–449, 2017.
  * [12] Y. Cai and J. Shen.  Error estimates for a fully discretized scheme to a Cahn-Hilliard phase-field model for two-phase incompressible flows.  Mathematics of Computation, 87(313):2057–2090, 2018.
  * [13] X. Chen.  Spectrum for the Allen-Chan, Chan-Hillard, and phase-field equations for generic interfaces.  Communications in Partial Differential Equations, 19(7-8):1371–1395, 1994.
  * [14] X. Chen.  Global asymptotic limit of solutions of the Cahn-Hilliard equation.  Journal of Differential Geometry, 44(2):262–311, 1996.
  * [15] Q. Du and R. A. Nicolaides.  Numerical analysis of a continuum model of phase transition.  SIAM Journal on Numerical Analysis, 28(5):1310–1322, 1991.
  * [16] C. Eck, B. Jadamba, and P. Knabner.  Error estimates for a finite element discretization of a phase field model for mixtures.  SIAM Journal on Numerical Analysis, 47(6):4429–4445, 2010.
  * [17] C. M. Elliott and D. A. French.  A nonconforming finite-element method for the two-dimensional Cahn-Hilliard equation.  SIAM Journal on Numerical Analysis, 26(4):884–903, 1989.
  * [18] C. M. Elliott, D. A. French, and F. Milner.  A second order splitting method for the Cahn-Hilliard equation.  Numerische Mathematik, 54(5):575–590, 1989.
  * [19] X. Feng and O. Karakashian.  Fully discrete dynamic mesh discontinuous Galerkin methods for the Cahn-Hilliard equation of phase transition.  Mathematics of Computation, 76(259):1093–1117, 2007.
  * [20] X. Feng and Y. Li.  Analysis of symmetric interior penalty discontinuous Galerkin methods for the Allen-Cahn equation and the mean curvature flow.  IMA Journal of Numerical Analysis, 35(4):1622–1651, 2015.
  * [21] X. Feng, Y. Li, and Y. Xing.  Analysis of mixed interior penalty discontinuous Galerkin methods for the Cahn-Hilliard equation and the hele–shaw flow.  SIAM Journal on Numerical Analysis, 54(2):825–847, 2016.
  * [22] X. Feng and A. Prohl.  Numerical analysis of the Allen-Cahn equation and approximation for mean curvature flows.  Numerische Mathematik, 94(1):33–65, 2003.
  * [23] X. Feng and A. Prohl.  Analysis of a fully discrete finite element method for the phase field model and approximation of its sharp interface limits.  Mathematics of Computation, 73(246):541–567, 2004.
  * [24] X. Feng and A. Prohl.  Error analysis of a mixed finite element method for the Cahn-Hilliard equation.  Numerische Mathematik, 99(1):47–84, 2004.
  * [25] X. Feng and H. Wu.  A posteriori error estimates for finite element approximations of the Cahn-Hilliard equation and the hele-shaw flow.  Journal of Computational Mathematics, pages 767–796, 2008.
  * [26] X. Feng and H.-j. Wu.  A posteriori error estimates and an adaptive finite element method for the Allen-Cahn equation and the mean curvature flow.  Journal of Scientific Computing, 24(2):121–146, 2005.
  * [27] R. L. Pego.  Front migration in the nonlinear Cahn-Hilliard equation.  Proceedings of the Royal Society of London. A. Mathematical and Physical Sciences, 422(1863):261–278, 1989.
  * [28] A. Prohl and X. H. Feng.  Numerical analysis of the Cahn-Hilliard equation and approximation for the hele-shaw problem.  Interfaces and Free Boundaries, 7(1):1–28, 2005.
  * [29] J. Shen and J. Xu.  Convergence and error analysis for the scalar auxiliary variable (SAV) schemes to gradient flows.  SIAM Journal on Numerical Analysis, 56(5):2895–2912, 2018.
  * [30] J. Shen, J. Xu, and J. Yang.  The scalar auxiliary variable (SAV) approach for gradient flows.  Journal of Computational Physics, 353:407–416, 2018.
  * [31] J. Shen, J. Xu, and J. Yang.  A new class of efficient and robust energy stable schemes for gradient flows.  SIAM Review, 61(3):474–506, 2019.
  * [32] J. Shen and X. Yang.  Numerical approximations of Allen-Cahn and Cahn-Hilliard equations.  Discrete & Continuous Dynamical Systems, 28(4):1669, 2010.
  * [33] L. Wang and H. Yu.  On efficient second order stabilized semi-implicit schemes for the Cahn-Hilliard phase-field equation.  Journal of Scientific Computing, 77(2):1185–1209, 2018.
  * [34] X. Yang.  Linear, first and second-order, unconditionally energy stable numerical schemes for the phase field model of homopolymer blends.  Journal of Computational Physics, 327:294–316, 2016.
  * [35] X. Yang and L. Ju.  Efficient linear schemes with unconditional energy stability for the phase field elastic bending energy model.  Computer Methods in Applied Mechanics and Engineering, 315:691–712, 2017.
  * [36] X. Yang and L. Ju.  Linear and unconditionally energy stable schemes for the binary fluid–surfactant phase field model.  Computer Methods in Applied Mechanics and Engineering, 318:1005–1029, 2017.
  * [37] X. Yang, J. Zhao, Q. Wang, and J. Shen.  Numerical approximations for a three-component Cahn-Hilliard phase-field model based on the invariant energy quadratization method.  Mathematical Models and Methods in Applied Sciences, 27(11):1993–2030, 2017.
  * [38] P. Yue, J. J. Feng, C. Liu, and J. Shen.  A diffuse-interface method for simulating two-phase flows of complex fluids.  Journal of Fluid Mechanics, 515:293–317, 2004.
  * [39] G. Zhang and X. Yang.  Error estimates with only a low-order polynomial dependence for the fully-discrete finite element invariant energy quadratization scheme of the Allen-Cahn equation.  To be appear, 2022.


# Evidence against the Efimov effect in ${}^{12}\mathrm{C}$ from spectroscopy
and astrophysics

J. Bishop<EMAIL_ADDRESS>Cyclotron Institute, Texas A&M University,
College Station, TX 77843, USA  Department of Physics & Astronomy, Texas A&M
University, College Station, TX 77843, USA  G.V. Rogachev  Cyclotron
Institute, Texas A&M University, College Station, TX 77843, USA  Department of
Physics & Astronomy, Texas A&M University, College Station, TX 77843, USA
Nuclear Solutions Institute, Texas A&M University, College Station, TX 77843,
USA  S. Ahn  Cyclotron Institute, Texas A&M University, College Station, TX
77843, USA  Center for Exotic Nuclear Studies, Institute for Basic Science,
Daejeon, 34126, Korea  E. Aboud  Cyclotron Institute, Texas A&M University,
College Station, TX 77843, USA  Department of Physics & Astronomy, Texas A&M
University, College Station, TX 77843, USA  M. Barbui  Cyclotron Institute,
Texas A&M University, College Station, TX 77843, USA  A. Bosh  Cyclotron
Institute, Texas A&M University, College Station, TX 77843, USA  Department of
Physics & Astronomy, Texas A&M University, College Station, TX 77843, USA  J.
Hooker  Department of Physics & Astronomy, University of Tennessee, Knoxville,
TN 37996, USA  C. Hunt  Cyclotron Institute, Texas A&M University, College
Station, TX 77843, USA  Department of Physics & Astronomy, Texas A&M
University, College Station, TX 77843, USA  H. Jayatissa  Physics Division,
Argonne National Laboratory, Argonne, IL 60439, USA  E. Koshchiy  Cyclotron
Institute, Texas A&M University, College Station, TX 77843, USA  R. Malecek
Department of Physics and Astronomy, Louisiana State University, Baton Rouge,
LA 70803, USA  S.T. Marley  Department of Physics and Astronomy, Louisiana
State University, Baton Rouge, LA 70803, USA  M. Munch  Department of Physics
and Astronomy, Aarhus University, DK-8000 Aarhus C, Denmark  E.C. Pollacco
IRFU, CEA, Saclay, Gif-Sur-Ivette, France  C.D. Pruitt  Department of
Chemistry, Washington University, St. Louis, MO 63130, USA  Lawrence Livermore
National Laboratory, Livermore, CA 94550, USA  B.T. Roeder  Cyclotron
Institute, Texas A&M University, College Station, TX 77843, USA  A.
Saastamoinen  Cyclotron Institute, Texas A&M University, College Station, TX
77843, USA  L.G. Sobotka  Department of Chemistry, Washington University, St.
Louis, MO 63130, USA  S. Upadhyayula  TRIUMF, 4004 Westbrook Mall, Vancouver,
British Columbia V6T 2A3, Canada

###### Abstract

Background

The Efimov effect is a universal phenomenon in physics whereby three-body
systems are stabilized via the interaction of an unbound two-body sub-systems.
A hypothetical state in ${}^{12}\mathrm{C}$ at 7.458 MeV excitation energy,
comprising of a loose structure of three $\alpha$-particles in mutual two-body
resonance, has been suggested in the literature to correspond to an Efimov
state in nuclear physics. The existence of such a state has not been
demonstrated experimentally.

Purpose

Using a combination of $\gamma$-spectroscopy, charged-particle spectroscopy
and astrophysical rate calculations, strict limits on the existence of such a
state have been established here.

Method

Using the combined data sets from two recent experiments, one with the TexAT
TPC to measure $\alpha$-decay and the other with Gammasphere to measure
$\gamma$-decay of states in ${}^{12}\mathrm{C}$ populated by
${}^{12}\mathrm{N}$ and ${}^{12}\mathrm{B}$ $\beta$-decay respectively, we
achieve high sensitivity to states in close-proximity to the
$\alpha$-threshold in ${}^{12}\mathrm{C}$.

Results

No evidence of a state at 7.458 MeV is seen in either data set. Using a
likelihood method, the 95% C.L. $\gamma$-decay branching ratio is determined
as a function of the $\beta$-decay feeding strength relative to the Hoyle
state. In parallel, calculations of the triple-alpha reaction rate show the
inclusion of the Efimov corresponds to a large increase in the reaction rate
around $5\times 10^{7}$ K.

Conclusion

From decay spectroscopy - at the 95% C.L., the Efimov state cannot exist at
7.458 MeV with any $\gamma$-decay branching ratio unless the $\beta$-strength
is less than 0.7% of the Hoyle state. This limit is evaluated for a range of
different excitation energies and the results are not favorable for existence
of the hypothetical Efimov state in ${}^{12}\mathrm{C}$. Furthermore, the
triple-alpha reaction rate with the inclusion of a state between 7.43 and 7.53
MeV exceeds the rate required for stars to undergo the red giant phase.

###### pacs:

Valid PACS appear here

††preprint: APS/123-QED

##  I Introduction

The Efimov effect is a universal quantum phenomena present in several areas of
Physics. The effect is observed for three-body systems that are comprised of
subsystems where the sub-unit two-body systems are unbound but have a large
s-wave scattering length. Vitaly Efimov found that under these conditions the
long-range three-body attraction arises and that this attraction can support a
family of three-body states. A detailed review of Efimov physics and
experimental evidence for it can be found in Ref. Naidon and Endo (2017). The
classical interpretation is that the three-body force binds the system via the
‘shuttling’ back and forth of one of the particles. This then creates an
infinite series of states given by a universal scaling law. While predicted 50
years ago, this phenomenon took 35 years to observe experimentally. The first
clear evidence for Efimov effect was reported for the system of ultracold gas
of caesium atoms in an external magnetic field Kraemer  _et al._ (2006). In
nuclear systems, which were the original focus of Efimov’s investigation, the
situation is more complicated. In principle, a $J^{\pi}=0^{+}$ (corresponding
to L=0) 3-$\alpha$ state in ${}^{12}\mathrm{C}$, where the 2-$\alpha$ sub-
systems are unbound but form a long-lived resonant state, can be seen as
Efimov trimer. In his original paper Efimov argued that the 0+ excited state
at 7.65 MeV in ${}^{12}\mathrm{C}$ (Hoyle state) possibly originates due to
this interesting three-body quantum phenomena which we now call the Efimov
effect Efimov (1970). Microscopic three-body continuum calculations that
utilize phenomenological $\alpha$-$\alpha$ potential Suno  _et al._ (2015)
indicate that the centrifugal (three-body) potential is small compared to the
nuclear and Coulomb potentials and as a result the Hoyle state is probably not
related to the Efimov effect. Faddeev calculations for the 3$\alpha$ system
have previously been performed and do not predict an additional near-threshold
$0^{+}$ that would correspond to an Efimov state Shinsho and Hiroyuki (1989).
More recently, existence of an Efimov state in ${}^{12}\mathrm{C}$ at an
excitation energy that corresponds to a mutual ${}^{8}\mathrm{Be}$(g.s.)
resonance for all three $\alpha$-particles was suggested in Zheng  _et al._
(2018a); Zheng and Bonasera (2020). This excitation energy is given, in units
of MeV, by the following simple relation which takes into account that the
narrow resonance of ${}^{8}\mathrm{Be}$(g.s.) is unbound by 91.84 keV with
respect to decay to two $\alpha$-particles:

$\displaystyle E=\frac{2}{3}\sum^{3}_{i}E_{rel_{i}}-Q=2\times
0.09184+7.2747=7.458.$ (1)

The main goal of this paper is to search for any evidence of this hypothetical
state or, if not, place experimental limits on its existence. We utilize the
results of two recent experiments. The first one is the study that is
sensitive to three $\alpha$-particles decay channel of the near
$\alpha$-threshold excited states in ${}^{12}\mathrm{C}$ populated in the
$\beta^{+}$ decay of ${}^{12}\mathrm{N}$ Bishop  _et al._ (2020a, b). The
second is the $\gamma$-spectroscopy study of states in ${}^{12}\mathrm{C}$
populated in $\beta^{-}$ decay of ${}^{12}\mathrm{B}$ Munch  _et al._ (2016).
Combining these two data sets we demonstrate that there is no evidence for a
resonance at or near 7.458 MeV. Stringent experimental limits on its existence
have been established. We also examine the astrophysical implications of such
a state existing below the Hoyle state using a simplistic model to calculate
the triple-alpha reaction rate and demonstrate the incompatibility of this
result with current astrophysical observations.

##  II Experimental limits of an Efimov state

There has not been a large amount of experimental activity to investigate if
there is a resonance in ${}^{12}\mathrm{C}$ below the Hoyle state. One
investigation claiming to observe an Efimov state in a heavy-ion Zn+Zn/Ni+Ni
collisions at 35 MeV/u was recently published Zhang  _et al._ (2019). By
examining events with 3 $\alpha$-particles and measuring their relative
energies, potential Efimov states were found by looking at 3 $\alpha$-particle
triplets where the relative energy between all $\alpha$-particles is
consistent with 92 keV. Due to the reaction mechanism used, there is a
dominant contribution from uncorrelated $\alpha$-particles which are required
to be accounted for via mixing. The remaining spectrum is then accounted for
using an arbitrary fit function on top of Breit Wigner peaks which show up
around 0.1 MeV suggested mutual ${}^{8}\mathrm{Be}$ resonances. Due to low
statistics and dependence on this fit function, it is therefore difficult to
definitely claim evidence of a peak. However, this study represented a
dedicated effort to experimentally observe such a state. While in-medium
effects afforded by a heavy-ion reaction may enhance the production of an
Efimov state, the reaction complexity also introduces numerous sources of
background that are difficult to account for. A cleaner population method is
through the $\beta$-decay of ${}^{12}\mathrm{N}/^{12}\mathrm{B}$ which
populates $0^{+}$, $1^{+}$, and $2^{+}$ states in accordance to the
$\beta$-decay selection rules. Through the recently-developed $\beta$-delayed
charged-particle decay technique using the TexAT TPC Bishop  _et al._ (2020a),
one is afforded a good probe on any resonance existing at a low relative
energy above the 3$\alpha$ threshold. This has not been previously identified
due to the difficulties associated with low-energy measurements using
implantation in silicon detector arrays. Despite the advantages afforded using
this technique, the measurement of a 3-particle final state only 180 keV above
threshold is still challenging.

###  II.1 Efimov state limits from $\beta$-delayed charged-particle decays in
a TPC

To establish a limit on the population of the Efimov state via $\beta$-decay
of ${}^{12}\mathrm{N}$, the TexAT TPC was used to study the $\beta$-delayed
charged-particle decay of ${}^{12}\mathrm{N}$ Bishop  _et al._ (2020a). A beam
of ${}^{12}\mathrm{N}$ with an energy of 24 MeV was stopped inside of the
active-area of the TexAT TPC by 20 Torr of CO2 gas. This was achieved on a
single implantation decay basis which allowed for the matching between the
implantation event and the subset of events where the $\beta$-decay populates
states above the $\alpha$-threshold in ${}^{12}\mathrm{C}$. The trigger
condition for the decay was that more than one Micromegas pad fired. The
trigger threshold is 10 keV and therefore one is sensitive to decay events
down to $E_{x}\sim$ 7.4 MeV. The 3D tracks from the 3 $\alpha$-particles
arising from these events were then reconstructed and the total energy
deposited was used to calculate the excitation energy. This is shown in Fig. 1
where the majority of events seen are from the Gaussian tails of the Hoyle
state at 7.654 MeV. One may also determine the energy of any 3$\alpha$ decay
by looking at the total length of the tracks and converting the range of the
$\alpha$-particles to an energy and summing. Taking events which lay within
the region of interest for both of these plots, one may then manually check
events to deduce the origin of these counts. The peak in Fig. 1 just above the
threshold corresponds to events where the implanting beam is scattered
(primarily off the entrance window) and the ${}^{12}\mathrm{N}$ beam is
implanted either on the cathode or anode. When it subsequently decays, part of
the energy of the event is lost as one or two $\alpha$-particles deposit their
energy into the anode or cathode directly rather than liberating electrons in
the gas volume. To determine the contribution from states between the absolute
threshold at the Hoyle state, a Gaussian tail was fitted for the Hoyle state
contribution using the previously obtained experimental resolution of 55 keV.
To remain conservative, any background above this Gaussian component was
counted as a potential Efimov contribution: 137.1 counts in total. The 95%
confidence limit (C.L.) was therefore taken as 160.5 for $\alpha$-events from
Poisson statistics.

To convert this into a limit, the expected number of 3$\alpha$ Efimov-state
decays is given by:

$\displaystyle
N_{\text{ES}}=N_{\text{Hoyle}}\frac{\beta_{\text{ES}}}{\beta_{\text{Hoyle}}}\frac{\text{BR}_{\alpha}(\text{ES})}{\text{BR}_{\alpha}(\text{Hoyle})},$
(2)

where $\beta_{i}$ describes the $\beta$-decay feeding strength from
${}^{12}\mathrm{N}$,
$\text{BR}_{\alpha}=\frac{\Gamma_{\alpha}}{\Gamma_{\text{tot}}}$ for the ES
and the Hoyle state. Re-arranging the terms, one arrives at:

$\displaystyle\frac{\beta_{\text{ES}}}{\beta_{\text{Hoyle}}}\text{BR}_{\alpha}(\text{ES})=\text{BR}_{\alpha}(\text{Hoyle})\frac{N_{\text{ES}}}{N_{\text{Hoyle}}}.$
(3)

The factor on the left-hand side of Eq. 3 allows us to place a limit on the
$\alpha$ branching ratio multiplied by the $\beta$-feeding strength relative
to the Hoyle state. Inserting the values,
$\text{BR}_{\alpha}(\text{Hoyle})=99.9\%$, $N_{\text{Hoyle}}=$ 23,276 and
$N_{\text{ES}}<160.5$. This 95% confidence limit is represented in Fig. 2
(slanted-right blue hash region) where
$\frac{\beta_{\text{ES}}}{\beta_{\text{Hoyle}}}$ is plotted against
$\text{BR}_{\gamma}\text{(ES)}=1-\text{BR}_{\alpha}(\text{ES})$. The product
of the $\beta$-strengths and the gamma branching-ratio
$(\frac{\beta_{\text{ES}}}{\beta_{\text{Hoyle}}}\text{BR}_{\alpha}(\text{ES}))$
is $<0.69\%$ from Eq. 3 therefore constraining $x(1-y)<0.69\%$ in the
exclusion plot.

Figure 1: Excitation function obtained from 3-$\alpha$ decay energy
deposition. The dotted red line corresponds to a Gaussian fit for the Hoyle
state contribution. The yield above this is conservatively taken as possibly
arising from additional states such as the Efimov state (ES). The peak around
7.3 MeV corresponds to deposition events on the cathode.

Figure 2: 95% C.L. exclusion plot for a 7.458 MeV Efimov state. The gamma and
alpha experimental limits exclude regions of the plot relating to the strength
of the $\beta$-feeding strength relative to the Hoyle state (abscissa) and the
gamma-ray branching ratio (ordinate). With a particular range of the gamma
branching ratio, the existence of a state cannot be conclusively excluded if
the $\beta$-feeding strength ratio is $<7\times 10^{-3}$.

###  II.2 Efimov state limits from $\beta$-delayed gamma-spectroscopy

Aside from observation in the 3$\alpha$ channel, there have been studies of
the $\gamma$-decay spectrum associated with $\beta$-delayed population of
${}^{12}\mathrm{C}$ using ${}^{12}\mathrm{B}$ Munch  _et al._ (2016). By using
Gammasphere to study for coincident gamma rays with the 4.44 MeV first-excited
state in ${}^{12}\mathrm{C}$, one may examine any possible contribution from
${}^{12}\mathrm{C}$(ES) $\rightarrow{{}^{12}\mathrm{C}(2_{1}^{+}})$ which
should correspond to $E_{\gamma}=(7.458-4.444)\,\mathrm{MeV}=3.014$ MeV.
Examining the region of interest from these previous data Munch  _et al._
(2016), no such yield can be seen in Fig. 3. The yield from a potential Efimov
state is $<11.4$ counts at 95% C.L. in comparison to the Hoyle state which has
a yield of $58~{}\pm~{}9$. Following the same prescription above, the number
of $\gamma$ rays from Efimov state decays is given by:

$\displaystyle
N_{\text{ES}}=N_{\text{Hoyle}}\frac{\beta_{\text{ES}}}{\beta_{\text{Hoyle}}}\frac{\text{BR}_{\gamma}(\text{ES})}{\text{BR}_{\gamma}(\text{Hoyle})}.$
(4)

Separating the unknowns onto the left-hand side from those known on the right-
hand side, one gets:

$\displaystyle\frac{\beta_{\text{ES}}}{\beta_{\text{Hoyle}}}\text{BR}_{\gamma}(\text{ES})=\text{BR}_{\gamma}(\text{Hoyle})\frac{N_{\text{ES}}}{N_{\text{Hoyle}}}.$
(5)

The known values are $N_{\text{ES}}<11.4$, $N_{\text{Hoyle}}$ = $58~{}\pm~{}9$
and $\text{BR}_{\gamma}(\text{Hoyle})=0.042\%$. This allows for an additional
exclusion region from the product of the $\beta$-strengths and the gamma
branching-ratio of $<0.008\%$ (red left-slanted hash region in Fig. 2).

Figure 3: Gammasphere energy spectrum gated on 4439 keV transition. The peak
at 3214 keV is from the Hoyle $\rightarrow 2_{1}^{+}$ transition. The yield in
the 3014 keV region (red arrow) shows no sign of a peak from the Efimov
$\rightarrow 2_{1}^{+}$ transition. Data from Munch  _et al._ (2016).

The combination of the gamma and alpha measurements in Fig. 2 demonstrates
that an Efimov state at 7.458 MeV is excluded at the 95% C.L. unless the
$\beta$-feeding strength is $<7\times 10^{-3}$ the strength of the Hoyle state
and has a specific decay $\gamma$-decay branch. Requiring the gamma branching
ratio to be $\approx 100\%$, this 95% C.L. extends down to $<10^{-4}$.

The existence of an Efimov state at an energy other than 7.458 MeV may also be
examined using the same approach as detailed above. The product of the
$\beta$-strengths and the gamma branching-ratio at the 95% C.L. does not
exceed 0.016%. The technique used for the alpha-limit is robust up until the
separation of the Hoyle peak from any potential second peak is roughly
2.35$\sigma$, i.e. up until $E_{x}$= 7.525 MeV where the gamma- and alpha-
branching ratios are comparable. The gamma-limit is highly applicable until
the gamma branching ratio is such that the predicted number of Efimov gamma
decays is less than the 95% C.L. or our detection limit. Given that Hoyle
state is well measured via gamma-decays in the Gammasphere data and the
background is fairly flat around the Hoyle peak, this implies an excitation
energy exceeding the Hoyle state is required whereby the gamma branching ratio
continues to rapidly decrease, and as such, this scenario may be definitely
excluded.

##  III Astrophysical limitations on low-lying states

In the triple-alpha process, there is a strong contribution from low-lying
$0^{+}$ states Ishikawa (2013); Nguyen  _et al._ (2012). It is for this reason
that the Hoyle state so successfully enhances the triple-alpha reaction rate
by 7 orders of magnitude, overcoming the A=5, 8 bottleneck in helium-burning
stars. Any additional low-lying state therefore must also contribute to a
large degree, particularly at lower temperatures. To understand the role that
the Efimov state may have, one can examine the expected astrophysical reaction
rate including such a resonance in addition to the Hoyle state. The NACRE
formulation was used Angulo  _et al._ (1999) and three cases were examined,
the results of which appear in Fig. 4.

The first case (dashed-dotted black line) is the NACRE parameter case Angulo
_et al._ (1999) which shows the current reaction rate. The second case (solid
blue line) uses the most up-to-date parameters for the Hoyle state, the
$2_{2}^{+}$ in ${}^{12}\mathrm{C}$ and the ${}^{8}\mathrm{Be}$ (g.s) Kibédi
_et al._ (2020); Kelley  _et al._ (2017); Tilley  _et al._ (2004). The main
influence corresponds to the modified radiative and alpha-widths of the Hoyle
state and the revised width of the ${}^{8}\mathrm{Be}$ (g.s). The final case
(dashed red line) is whereby, in addition to the states included in the second
case, an additional resonance corresponding to the Efimov state was included.
For such a state, in Table 1, naïve estimates of the widths are calculated via
assuming the same underlying reduced width and radiative width for the Efimov
state as the Hoyle state when scaled appropriately by their energies (for
penetrability and $E_{\gamma}^{5}$ for the $\alpha$ and $\gamma$ widths
accordingly). This is verified by more developed investigations into the
properties of an Efimov state Zheng  _et al._ (2018a). For the Efimov state,
$\Gamma_{\text{rad}}\gg\Gamma_{\alpha_{0}}$ so
$\Gamma_{\text{rad}}\approx\Gamma_{\text{tot}}$. Therefore, the rate is driven
predominantely by $\Gamma_{\alpha_{0}}$.

The contribution of the additional Efimov resonance in the third case clearly
demonstrates a phenomenal increase (a factor of $\sim$ 40,000 at $5\times
10^{7}$ K ) in the reaction rate at lower temperatures which are vitally
important for many stars. Understanding the role that the triple-alpha
reaction rate has on the dynamics of stars has previously been examined Suda
_et al._ (2011) and it was concluded that the reaction rate at $10^{7.8}$ K
must be less than $\sim 10^{-29}$ cm-6 s-1 mol-2 in order for stars to undergo
the red giant phase. This astrophysical limit on the rate is shown in Fig. 4
as an upside-down green triangle and the rate from case three (including the
Efimov state) can clearly be seen to exceed this limit and does so by a factor
of $\sim 24$. As such, the existence of such a state is incompatible with
current astrophysical models. Given the contribution of the Efimov state is
limited by $\Gamma_{\alpha}$, to quench the contribution of this state to be
compatible with the reaction rate limit, this width is required to be
decreased by an order of magnitude. The Efimov state has a well-developed
cluster structure therefore a partial alpha-width a factor of ten lower than
the Hoyle state is clearly not possible. Such a state is additionally excluded
down to an extremely small beta-feeding ratio of $10^{-4}$ from the combined
$\alpha/\gamma$-spectroscopy for such a state.

Figure 4: Triple-alpha reaction rates as a function of the astrophysical
temperature using the NACRE formulism Angulo  _et al._ (1999). The dot-dashed
black line corresponds to the pre-existing NACRE rate. The solid blue line
includes the Hoyle state and the $2_{2}^{+}$ state in ${}^{12}\mathrm{C}$
using the latest accepted values. The third line, dashed red, also includes an
additional state corresponding to the Efimov state using parameters from Table
1. The upper limit calculated from Suda Suda  _et al._ (2011) is shown as a
green triangle whereby the red giant phase requires a reaction rate of less
than $\sim 10^{-29}$ cm-6 s-1 mol-2 at a temperature of $10^{7.8}$ K. The
inclusion of the Efimov state clearly violates this reaction rate limit.

Modifying the energy of this additional state results in a rate of more than
$10^{-29}$ cm-6 s-1 mol-2 for energies of $7.43<E_{x}<7.53$ MeV therefore
resonances in this region may also be explicitly excluded.

##  IV Conclusion

Using combined data from particle and gamma-decay from
${}^{12}\mathrm{N}/^{12}\mathrm{B}$, one may place incredibly restrictive
limits on the beta-feeding strength and therefore existence of an Efimov state
below the Hoyle state in carbon-12. An Efimov state in such close proximity to
the alpha-threshold would decay almost solely via gamma-decay and therefore
the beta-feeding strength relative to the Hoyle state is expected to be
$<10^{-4}$ at the 95% C.L. Furthermore, it can clearly be seen that the Efimov
state has a tremendous impact on stellar evolution. A low-lying state $0^{+}$
(whether Efimov or merely an $\alpha$-clustered state) would have a
significant impact on stellar abundances that the existence of such a
resonance is incompatible with current astrophysical models whereby the
modified triple-alpha reaction rate would eliminate the red-giant phase of
stars.

Table 1: Naïve estimates of the Efimov decay mechanism. The $\alpha$ widths, $\Gamma_{\alpha}$, assume the same reduced width and the radiative width assumes the dominance of the E2 transition to the 4.44 MeV $2^{+}$ state and is scaled as $E_{\gamma}^{5}$. Parameter | Hoyle state | Efimov state
---|---|---
$E_{x}$ | 7.654 MeV | 7.458 MeV
$\Gamma_{\alpha}$ | 9.3 eV Kelley  _et al._ (2017) | 79.8 neV
$\Gamma_{\mathrm{rad}}$ | 5.1 meV Kibédi  _et al._ (2020) | 3.72 meV
$\Gamma_{\alpha}/\Gamma_{\mathrm{tot}}$ | $>99.9\%$ | $0.002\%$

##  V Acknowledgements

J.B. acknowledges and thanks A. Bonasera for discussions on the Efimov state.
This work was supported by the U.S. Department of Energy, Office of Science,
Office of Nuclear Science, under award no. DE-FG02-93ER40773 and by National
Nuclear Security Administration through the Center for Excellence in Nuclear
Training and University Based Research (CENTAUR) under grant number DE-
NA0003841. G.V.R. also acknowledges the support of the Nuclear Solutions
Institute.

## References

  * Naidon and Endo (2017) P. Naidon and S. Endo, Reports on Progress in Physics 80, 056001 (2017).
  * Kraemer  _et al._ (2006) T. Kraemer, M. Mark, P. Waldburger, J. G. Danzl, C. Chin, B. Engeser, A. D. Lange, K. Pilch, A. Jaakkola, H.-C. Nägerl, and R. Grimm, Nature 440, 315 (2006).
  * Efimov (1970) V. Efimov, Phys. Lett. B 33, 563 (1970).
  * Suno  _et al._ (2015) H. Suno, Y. Suzuki, and P. Descouvemont, Phys. Rev. C 91, 014004 (2015).
  * Shinsho and Hiroyuki (1989) O. Shinsho and K. Hiroyuki, Nuclear Physics A 493, 91 (1989).
  * Zheng  _et al._ (2018a) H. Zheng, A. Bonasera, M. Huang, and S. Zhang, Phys. Lett. B 779, 460 (2018a).
  * Zheng and Bonasera (2020) H. Zheng and A. Bonasera, Journal of Physics Communications 4, 085011 (2020).
  * Bishop  _et al._ (2020a) J. Bishop, G. V. Rogachev, S. Ahn, E. Aboud, M. Barbui, P. Baron, A. Bosh, E. Delagnes, J. Hooker, C. Hunt, H. Jayatissa, E. Koshchiy, R. Malecek, S. Marley, R. O’Dwyer, E. Pollacco, C. Pruitt, B. Roeder, A. Saastamoinen, L. Sobotka, and S. Upadhyayula, Nuclear Instruments and Methods in Physics Research Section A: Accelerators, Spectrometers, Detectors and Associated Equipment 964, 163773 (2020a).
  * Bishop  _et al._ (2020b) J. Bishop, G. V. Rogachev, S. Ahn, E. Aboud, M. Barbui, A. Bosh, C. Hunt, H. Jayatissa, E. Koshchiy, R. Malecek, S. T. Marley, E. C. Pollacco, C. D. Pruitt, B. T. Roeder, A. Saastamoinen, L. G. Sobotka, and S. Upadhyayula, Phys. Rev. C 102, 041303(R) (2020b).
  * Munch  _et al._ (2016) M. Munch, M. Alcorta, H. O. U. Fynbo, M. Albers, S. Almaraz-Calderon, M. L. Avila, A. D. Ayangeakaa, B. B. Back, P. F. Bertone, P. F. F. Carnelli, M. P. Carpenter, C. J. Chiara, J. A. Clark, B. DiGiovine, J. P. Greene, J. L. Harker, C. R. Hoffman, N. J. Hubbard, C. L. Jiang, O. S. Kirsebom, T. Lauritsen, K. L. Laursen, S. T. Marley, C. Nair, O. Nusair, D. Santiago-Gonzalez, J. Sethi, D. Seweryniak, R. Talwar, C. Ugalde, and S. Zhu, Phys. Rev. C 93, 065803 (2016).
  * Zhang  _et al._ (2019) S. Zhang, A. Bonasera, M. Huang, H. Zheng, D. X. Wang, J. C. Wang, L. Lu, G. Zhang, Z. Kohley, Y. G. Ma, and S. J. Yennello, Phys. Rev. C 99, 044605 (2019).
  * Ishikawa (2013) S. Ishikawa, Phys. Rev. C 87, 055804 (2013).
  * Nguyen  _et al._ (2012) N. B. Nguyen, F. M. Nunes, I. J. Thompson, and E. F. Brown, Phys. Rev. Lett. 109, 141101 (2012).
  * Angulo  _et al._ (1999) C. Angulo, M. Arnould, M. Rayet, P. Descouvemont, D. Baye, C. Leclercq-Willain, A. Coc, S. Barhoumi, P. Aguer, C. Rolfs, R. Kunz, J. Hammer, A. Mayer, T. Paradellis, S. Kossionides, C. Chronidou, K. Spyrou, S. Degl’Innocenti, G. Fiorentini, B. Ricci, S. Zavatarelli, C. Providencia, H. Wolters, J. Soares, C. Grama, J. Rahighi, A. Shotter, and M. Lamehi Rachti, Nuclear Physics A 656, 3 (1999).
  * Kelley  _et al._ (2017) J. Kelley, J. Purcell, and C. Sheu, Nuclear Physics A 968, 71 (2017).
  * Tilley  _et al._ (2004) D. Tilley, J. Kelley, J. Godwin, D. Millener, J. Purcell, C. Sheu, and H. Weller, Nuclear Physics A 745, 155 (2004).
  * Zheng  _et al._ (2018b) H. Zheng, A. Bonasera, M. Huang, and S. Zhang, Physics Letters B 779, 460 (2018b).
  * Suda  _et al._ (2011) T. Suda, R. Hirschi, and M. Y. Fujimoto, The Astrophysical Journal 741, 61 (2011).
  * Kibédi  _et al._ (2020) T. Kibédi, B. Alshahrani, A. E. Stuchbery, A. C. Larsen, A. Görgen, S. Siem, M. Guttormsen, F. Giacoppo, A. I. Morales, E. Sahin, G. M. Tveten, F. L. B. Garrote, L. C. Campo, T. K. Eriksen, M. Klintefjord, S. Maharramova, H.-T. Nyhus, T. G. Tornyi, T. Renstrøm, and W. Paulsen, Phys. Rev. Lett. 125, 182701 (2020).