Stability and convergence of relaxed scalar auxiliary variable schemes for Cahn–Hilliard systems with bounded mass source

Kei Fong Lam ¹¹1Department of Mathematics, Hong Kong Baptist University, Kowloon Tong, Hong Kong (

\{

akflam,wangru_22

\}

@math.hkbu.edu.hk). Ru Wang ¹¹1Department of Mathematics, Hong Kong Baptist University, Kowloon Tong, Hong Kong (

\{

akflam,wangru_22

\}

@math.hkbu.edu.hk).

Abstract

The scalar auxiliary variable (SAV) approach of Shen et al. (2018), which presents a novel way to discretize a large class of gradient flows, has been extended and improved by many authors for general dissipative systems. In this work we consider a Cahn–Hilliard system with mass source that, for image processing and biological applications, may not admit a dissipative structure involving the Ginzburg–Landau energy. Hence, compared to previous works, the stability of SAV-discrete solutions for such systems is not immediate. We establish, with a bounded mass source, stability and convergence of time discrete solutions for a first-order relaxed SAV scheme in the sense of Jiang et al. (2022), and apply our ideas to Cahn–Hilliard systems appearing in diblock co-polymer phase separation, tumor growth, image inpainting and segmentation.

Keywords: Scalar auxiliary variable; Cahn–Hilliard equation; mass source; energy stability; convergence analysis AMS (MOS) Subject Classification: 35K35, 35K55, 65M12, 65Z05

1 Introduction

We are interested in the stability and convergence of numerical schemes based on the scalar auxiliary variable (SAV) approach [34, 35] for Cahn–Hilliard systems of the form:


$\displaystyle\partial_{t}\varphi$	$\displaystyle=\mathrm{div}\,(m(\varphi)\nabla\mu)+f\quad$	$\displaystyle\text{ in }\Omega\times(0,T),$	(1.1a)
$\displaystyle\mu$	$\displaystyle=-\varepsilon\Delta\varphi+\varepsilon^{-1}F^{\prime}(\varphi)-g\quad$	$\displaystyle\text{ in }\Omega\times(0,T).$	(1.1b)

In (1.1), $\varphi$ and $\mu$ denote the phase field variable and its corresponding chemical potential, $\varepsilon>0$ is a small fixed parameter, $m(\varphi)$ is a phase-dependent mobility function, $F$ is a potential function with derivative $F^{\prime}$ , and $f,g$ are source terms that can be prescribed or may depend on $\varphi$ . We furnish (1.1a)-(1.1b) with homogeneous Neumann conditions $\partial_{\bm{n}}\varphi=\partial_{\bm{n}}\mu=0$ on $\partial\Omega\times(0,T)$ , and initial condition $\varphi(0)=\varphi_{0}$ in $\Omega$ . The above system admits the following energy identity:

\displaystyle\frac{d}{dt}E_{\varepsilon}(\varphi)+\int_{\Omega}m(\varphi)|\nabla\mu|^{2}\,\mathrm{dx}=\int_{\Omega}f\mu+g\partial_{t}\varphi\,\mathrm{dx}

(1.2)

involving the Ginzburg–Landau functional $E_{\varepsilon}(\varphi)=\int_{\Omega}\frac{\varepsilon}{2}|\nabla\varphi|^{2}+\frac{1}{\varepsilon}F(\varphi)\,\mathrm{dx}$ .

Under suitable choices of $f$ and $g$ , such types of Cahn–Hilliard systems have seen recent applications in image inpainting [3], image segmentation [41], phase separation of diblock copolymer [19, 32, 31], evolution of brain metabolites concentration [29], and tumor growth [9, 17, 40]. Thus, efficient, stable and convergent numerical schemes for systems such as (1.1) are necessary to provide meaningful simulations to these areas of applications. We refer to [1, 12] for some previous numerical contributions for (1.1) with $g=0$ .

In the absence of the source terms $f$ and $g$ , and the mobility $m(\varphi)$ set to 1, the resulting Cahn–Hilliard equation can be viewed as a $H^{-1}$ -gradient flow of the Ginzburg–Landau functional $E_{\varepsilon}(\varphi)$ . The only nonlinearity is the derivative term $F^{\prime}(\varphi)$ , where a typical choice for the potential function $F$ is the quartic polynomial $F(s)=\frac{1}{4}(s^{2}-1)^{2}$ that has two minima at $s=\pm 1$ . Various unconditionally stable numerical schemes have been proposed to approximate similar gradient-flow types of phase field equations, among which we mention the Convex-concave splitting (CCS) approach [10, 11]; the Stabilized linearly implicit (SLI) approach [37]; the Invariant energy quadratization (IEQ) approach [6] (which is based on an earlier Lagrange multiplier approach [21]); and the Scalar auxiliary variable approach [34, 35] which is based on the simple idea of introducing a scalar variable

\displaystyle q(t)=\Big{(}\int_{\Omega}\frac{1}{\varepsilon}F(\varphi(t))\,\mathrm{dx}+C_{0}\Big{)}^{1/2}=:Q(\varphi(t))>0,

(1.3)

for potentials $F$ such that $\int_{\Omega}\frac{1}{\varepsilon}F(\varphi)\,\mathrm{dx}$ is bounded strictly from below by $-C_{0}$ with constant $C_{0}>0$ , so that an ordinary differential equation can be obtained:

q_{t}=\frac{1}{2\varepsilon Q(\varphi)}\int_{\Omega}F^{\prime}(\varphi)\partial_{t}\varphi\,\mathrm{dx},\quad q(0)=Q(\varphi_{0}).

The equivalent SAV formulation of (1.1) in strong form reads as


$\displaystyle\partial_{t}\varphi$	$\displaystyle=\mathrm{div}\,(m(\varphi)\nabla\mu)+f\quad$	$\displaystyle\text{ in }\Omega\times(0,T),$	(1.4a)
$\displaystyle\mu$	$\displaystyle=-\varepsilon\Delta\varphi+\frac{q(t)}{\varepsilon Q(\varphi)}F^{\prime}(\varphi)-g\quad$	$\displaystyle\text{ in }\Omega\times(0,T),$	(1.4b)
$\displaystyle q_{t}$	$\displaystyle=\frac{1}{2\varepsilon Q(\varphi)}\int_{\Omega}F^{\prime}(\varphi)\partial_{t}\varphi\,\mathrm{dx}\quad$	$\displaystyle\text{ in }(0,T),$	(1.4c)

which admits an analogous energy identity to (1.2) but with $E_{\varepsilon}$ replaced by a modified energy $G_{\varepsilon}=\frac{\varepsilon}{2}\|\nabla\varphi\|^{2}+|q|^{2}$ . A first order time discretization would take $F^{\prime}(\varphi)/(\varepsilon Q(\varphi))$ explicitly and $q$ implicitly when approximating (1.4b), yielding a numerical scheme that is linear in the next iteration for $(\varphi,\mu,q)$ . In comparison to the CCS approach which results in nonlinear discrete systems, or the SLI approach which requires a quadratic truncation of $F$ and large stabilization parameters, the IEQ and SAV approaches allow for unconditionally stable linear schemes when applied to gradient flows [35, 38] without the need to modify the potential $F$ . Between the latter two, the SAV approach is viewed as an enhancement of the IEQ approach, whose main idea of introducing an auxiliary function $u=\sqrt{F(\varphi)}$ and replacing $(q(t)F^{\prime}(\varphi))/(\varepsilon Q(\varphi))$ with $(F^{\prime}(\varphi)u)/(\varepsilon\sqrt{F(\varphi)})$ in (1.4b), and replacing (1.4c) with a PDE $\partial_{t}u=(F^{\prime}(\varphi)\partial_{t}\varphi)/(2\sqrt{F(\varphi)})$ demands $F$ to be a strictly positive potential function.

Since its introduction, many variants and improvements of the SAV approach have been developed, among which we mention the Lagrange multiplier approach [5] that dissipates the original energy as opposed to a modified energy; the relaxed SAV approach [15] penalizing large differences between $q^{n}$ and $Q(\varphi^{n})$ at the discrete level; SAV variants [23, 27] removing the lower bound assumption on $\int_{\Omega}F(\varphi)\,\mathrm{dx}$ ; generalized SAV (GSAV) [14, 42] and GSAV with relaxation (R-GSAV) [43] approaches applicable to general dissipative systems of the form

\displaystyle\partial_{t}\phi+\mathcal{A}(\phi)+g(\phi)=0\quad\text{ such that }\quad\frac{d}{dt}E(\phi)=-\mathcal{K}(\phi),

(1.5)

with positive operator $\mathcal{A}$ , semi-linear/quasi-linear operator $g$ , free energy functional $E$ and dissipative functional $\mathcal{K}$ such that $\mathcal{K}(\phi)>0$ for all $\phi$ .

Concerning the numerical analysis, assuming sufficient regularity on the solution to the underlying PDE, which is either a gradient flow of some energy functional or a dissipative system of the form (1.5), first error estimates for the standard SAV approach were established in [4, 33] for the Allen–Cahn and Cahn–Hilliard equations; in [26] for a Cahn–Hilliard–Stokes system; in [44] for a Cahn–Hilliard-Hele-Shaw system. For the GSAV and R-GSAV approaches, error analysis with backwards differentiation formula (BDF) time stepping can be found in [13, 43]. On the other hand, convergence of discrete solutions without the assumed smoothness of the continuous solution can be found in [33] for the $L^{2}$ - and $H^{-1}$ -gradient flow of $E_{\varepsilon}$ , and in [28] for the bulk-surface Cahn–Hilliard system of [25].

The motivation of the present work stems from the observations that (i) Cahn–Hilliard systems with mass sources are becoming ubiquitous in various areas of applied sciences; (ii) for a generic source term $f=f(x)$ or $f=f(x,\varphi)$ , the system (1.1) is no longer a gradient flow of $E_{\varepsilon}$ , and thus new ideas are needed to establish the stability of SAV discrete solutions that was automatically guaranteed in previous works. In our analysis we notice that if $f\in L^{2}(\Omega)$ , or $|f(\varphi)|\leq C(1+|\varphi|)$ having linear growth, then stability estimates require the potential $F$ to have at most quadratic growth (see Section 2.3). This restriction has also been observed in [16, 18] for a Cahn–Hilliard tumor model (5.11) with non-constant mobility function $m(\varphi)$ . This is similar to the SLI approach where $F$ needs to be modified so that $\|F^{\prime\prime}\|_{L^{\infty}}<\infty$ . However, as we recall one of the advantages of the SAV approach over the SLI approach is that $F$ needs no ad-hoc modification, the basis of our contribution is to provide stability estimates for SAV discrete solutions to (1.1) whilst keeping $F$ unmodified, which can be attained when $f$ is a bounded function. This boundedness assumption is not restrictive, as seen from several applications of (1.1) outlined in Section 5, and it is possible to extend our stability analysis to Cahn–Hilliard systems with non-bounded sources $f=f(\varphi)$ provided the time evolution of the mean value of $\varphi$ can be controlled.

Our main contributions are as follows: (i) we formulate a linear numerical scheme for (1.1) based on the RSAV approach of [15], and demonstrate first results on the stability of discrete solutions for bounded $f$ and potential $F$ that admits quartic polynomial growth; (ii) we prove first convergence results for the time discrete solutions to a weak solution of (1.1). To the best of our knowledge, SAV schemes for (1.1) has not received much attention in the literature, which we attribute to the lack of an obvious Lyapunov functional, or a dissipative structure akin to (1.5) for the equation. Furthermore, we are not aware of any results concerning the convergence of RSAV schemes of [15] (besides the R-GSAV variant in [43] for dissipative systems). As a consequence, notable modified Cahn–Hilliard models in the literature for diblock co-polymer phase separation, tumor growth, image inpainting and segmentation can now be approach numerically with the SAV framework, thereby extending its applicability to systems that need not exhibit a dissipative structure.

Plan of the paper: We present the RSAV time discretization for (1.1) in Section 2, and establish stability estimates. Section 3 is devoted to the convergence analysis of the time discrete solutions, and we discuss supporting numerical results in Section 4. In Section 5 we outline SAV schemes, their stability and numerical simulations for Cahn–Hilliard system with mass source in biological science and image processing.

Notation: For any $p\in[1,\infty]$ and $k>0$ , the standard Lebesgue and Sobolev spaces over $\Omega$ are denoted by $L^{p}(\Omega)$ and $W^{k,p}(\Omega)$ with the corresponding norms $\|\cdot\|_{L^{p}}$ and $\|\cdot\|_{W^{k,p}}$ . In the special case $p=2$ , these become Hilbert spaces and we employ the notation $H^{k}:=H^{k}(\Omega)=W^{k,2}(\Omega)$ with the corresponding norm $\|\cdot\|_{H^{k}}$ . We denote the topological dual of $H^{1}(\Omega)$ by $(H^{1}(\Omega))^{*}$ and the corresponding duality pairing by $\langle\cdot,\cdot\rangle$ . We use $\|\cdot\|$ and $(\cdot,\cdot)$ for the norm and inner product of $L^{2}(\Omega)$ .

2 RSAV scheme

2.1 First order time discretization

Dividing the time interval $[0,T]$ into a uniform partition of subintervals $[t^{n-1},t^{n}]$ with $\tau=t^{n}-t^{n-1}$ , for $n=1,\dots,N_{\tau}$ , where $\tau N_{\tau}\leq T$ . Denoting by $f^{k},g^{k}$ as approximations of $f,g$ evaluated at $t=t^{k}$ for $k=1,\dots,N_{\tau}$ , then a first order time discretization of (1.4) based on the SAV approach with relaxation proposed in [15] reads as follows: Let $\tau>0$ denote a fixed time step and given $(\varphi^{n-1},q^{n-1},f^{n-1},g^{n-1})$ , find $(\varphi^{n},\mu^{n},r^{n})$ satisfying (in a suitable sense)


$\displaystyle\varphi^{n}-\varphi^{n-1}$	$\displaystyle=\tau\,\mathrm{div}\,(m(\varphi^{n-1})\nabla\mu^{n})+\tau f^{n-1},$	(2.1a)
$\displaystyle\mu^{n}$	$\displaystyle=-\varepsilon\Delta\varphi^{n}+\frac{r^{n}}{\varepsilon Q(\varphi^{n-1})}F^{\prime}(\varphi^{n-1})-g^{n-1},$	(2.1b)
$\displaystyle r^{n}-q^{n-1}$	$\displaystyle=\frac{1}{2\varepsilon Q(\varphi^{n-1})}(F^{\prime}(\varphi^{n-1}),\varphi^{n}-\varphi^{n-1}),$	(2.1c)

initialized with $\varphi^{0}=\varphi_{0}$ and $q^{0}=Q(\varphi_{0})$ . Then, we define

\displaystyle q^{n}:=\zeta_{n}^{\tau}r^{n}+(1-\zeta_{n}^{\tau})Q(\varphi^{n})

(2.2)

for some $\zeta_{n}^{\tau}\in\mathcal{V}_{n}^{\tau}$ , where given $\tau,\mu^{n},r^{n}$ and manually assigned fixed constants $\eta\in(0,1)$ and $M>0$ , the set $\mathcal{V}_{n}^{\tau}$ is defined as the set of constants $\zeta\in[0,1]$ such that

	$\displaystyle 0\geq R_{n}^{\tau}(\zeta)$	$\displaystyle:=\|\zeta r^{n}+(1-\zeta)Q(\varphi^{n})\|^{2}+\frac{1}{2}\|\zeta r^{n}+(1-\zeta)Q(\varphi^{n})-q^{n-1}\|^{2}$		(2.3)
		$\displaystyle\quad-\tau\eta m_{0}\\|\nabla\mu^{n}\\|^{2}-\|r^{n}\|^{2}-\frac{1}{2}\|r^{n}-q^{n-1}\|^{2}-\tau M,$		(2.3)

where $m_{0}$ is the positive lower bound on the mobility function, see ( $\mathrm{A2}$ ) below. Our definition of $\mathcal{V}_{n}^{\tau}$ differs slightly from [15]. The standard SAV method is obtained when $\zeta=1$ . Note that $R_{n}(1)\leq-\tau M<0$ and so $1\in\mathcal{V}_{n}^{\tau}$ , i.e., $\mathcal{V}_{n}^{\tau}$ is non-empty. In fact, $R_{n}^{\tau}(\zeta)=a_{n}^{\tau}\zeta^{2}+b_{n}^{\tau}\zeta+c_{n}^{\tau}$ with coefficients

	$\displaystyle a_{n}^{\tau}$	$\displaystyle=\frac{3}{2}(r^{n}-Q(\varphi^{n}))^{2},\quad b_{n}^{\tau}=(r^{n}-Q(\varphi^{n}))(Q(\varphi^{n})-q^{n-1}),$
	$\displaystyle c_{n}^{\tau}$	$\displaystyle=\|Q(\varphi^{n})\|^{2}-\|r^{n}\|^{2}+\frac{1}{2}[\|Q(\varphi^{n})-q^{n-1}\|^{2}-\|r^{n}-q^{n-1}\|^{2}]-\tau\eta m_{0}\\|\nabla\mu^{n}\\|^{2}-\tau M.$

Since $R_{n}^{\tau}(1)<0$ , by continuity we can find $w_{n}^{\tau}\in[0,1)$ such that $R_{n}^{\tau}(\zeta)\leq 0$ for all $\zeta\in[w_{n}^{\tau},1]$ , and so we can pick any $\zeta_{n}^{\tau}\in(w_{n}^{\tau},1]$ for the update rule (2.2).

Remark 2.1.

It would be ideal to pick $\zeta_{n}^{\tau}=0$ so that $q^{n}=Q(\varphi^{n})$ , which we term as the idealized update approach when discussing convergence in Section 3 below. However, for the Cahn–Hilliard system (1.1) with generic mass source $f$ , we cannot guarantee if $R_{n}^{\tau}(0)\leq 0$ is satisfied for all $n=1,\dots,N_{\tau}$ as there is no available comparison between $|Q(\varphi^{n})|^{2}$ and $|r^{n}|^{2}$ using the discrete energy inequality (2.7) of (2.1). For the numerical simulations performed in [15], the following optimal value for $\zeta_{n}^{\tau}$ is suggested:

\displaystyle\zeta_{n}^{\tau}=\max\Big{(}0,\frac{-b_{n}^{\tau}-\sqrt{(b_{n}^{\tau})^{2}-4a_{n}^{\tau}c_{n}^{\tau}}}{2a_{n}^{\tau}}\Big{)},

(2.4)

where the values of $a_{n}^{\tau}$ , $b_{n}^{\tau}$ and $c_{n}^{\tau}$ at each iteration are defined above. In our numerical tests reported in Section 4, we found that the optimal value $\zeta_{n}^{\tau}=0$ can be achieved with moderate values of $M$ and $\eta$ .

2.2 Stability estimates

We make the following assumptions:

$(\mathrm{A1})$

$\Omega\subset\mathbb{R}^{d}$ , $d\in\{2,3\}$ is a bounded domain with convex or $C^{1,1}$ -boundary $\partial\Omega$ .
$(\mathrm{A2})$

The constant $\varepsilon>0$ is fixed, and the mobility function $m\in C^{0}(\mathbb{R})$ satisfies

$m_{0}\leq m(s)\leq m_{1}\quad\forall s\in\mathbb{R}$

for some constants $m_{1}\geq m_{0}>0$ .
$(\mathrm{A3})$

The initial condition satisfies $\phi_{0}\in H^{1}(\Omega)$ .
$(\mathrm{A4})$

The function $g:\Omega\times(0,T)\to\mathbb{R}$ satisfies $g\in L^{\infty}(0,T;L^{2}(\Omega))$ .
$(\mathrm{A5})$

The potential $F$ is a non-negative function with $F\in C^{1}(\mathbb{R})$ , and there exist constants $c_{0}>0$ , $c_{1}\in\mathbb{R}$ , $c_{2}>0$ such that

$c_{0}|s|^{3}-c_{1}\leq|F^{\prime}(s)|\leq c_{2}(1+|s|^{3})\quad\forall s\in\mathbb{R},$

and the source term $f:\Omega\times(0,T)\to\mathbb{R}$ satisfies $f\in L^{\infty}(0,T;L^{\infty}(\Omega))$ .

In light of ( $\mathrm{A4}$ ) and ( $\mathrm{A5}$ ) for the source functions, we introduce the Steklov averages $f^{\tau}$ and $g^{\tau}$ defined as

f^{\tau}(t,x)=\frac{1}{\tau}\int_{t}^{t+\tau}\tilde{f}(y,x)\,\mathrm{dy},\quad g^{\tau}(t,x)=\frac{1}{\tau}\int_{t}^{t+\tau}\tilde{g}(y,x)\,\mathrm{dy}

for $t\in[0,T-\tau]$ , where $\tilde{f}$ and $\tilde{g}$ are zero extensions of $f$ and $g$ on $t\in\mathbb{R}\setminus[0,T]$ . Then, well-known properties of Steklov averages yield $f^{\tau}\in C^{0}([0,T];L^{\infty}(\Omega))$ , $g^{\tau}\in C^{0}([0,T];L^{2}(\Omega))$ , $\|f^{\tau}\|_{L^{\infty}(0,T;L^{\infty})}\leq\|f\|_{L^{\infty}(0,T;L^{\infty})}$ , $\|g^{\tau}\|_{L^{\infty}(0,T;L^{2})}\leq\|g\|_{L^{\infty}(0,T;L^{2})}$ , with $f^{\tau}\to f$ in $L^{r}(0,T;L^{\infty}(\Omega))$ and $g^{\tau}\to g$ in $L^{r}(0,T;L^{2}(\Omega))$ as $\tau\to 0$ for any $1\leq r<\infty$ . Hence, we can choose $f^{n-1}(x):=f^{\tau}(t^{n-1},x)$ and $g^{n-1}(x):=g^{\tau}(t^{n-1},x)$ in (2.1).

Remark 2.2 (Unique solvability).

The existence of time discrete solutions to (2.1) can be inferred almost analogously as in Section 4.1 by replacing the corresponding matrices with appropriate differential operators, and so we omit the details. Concerning uniqueness of solutions, as (2.1) is linear in $(\varphi^{n},\mu^{n},r^{n})$ , the differences denoted by $(\varphi,\mu,r)$ between any two solution triplets satisfy formally in the strong sense

\varphi=\tau\mathrm{div}\,(m(\varphi^{n-1})\nabla\mu),\quad\mu=-\varepsilon\Delta\varphi+\frac{rF^{\prime}(\varphi^{n-1})}{\varepsilon Q(\varphi^{n-1})},\quad r=\frac{(F^{\prime}(\varphi^{n-1}),\varphi)}{2\varepsilon Q(\varphi^{n-1})}.

Testing with $\mu$ , $\varphi$ and $2r$ , respectively, and upon summing we get

\tau\int_{\Omega}m(\varphi^{n})|\nabla\mu|^{2}\,\mathrm{dx}+\varepsilon\|\nabla\varphi\|^{2}+2|r|^{2}=0.

Hence $r=0$ , while $\varphi$ and $\mu$ are constants. Since $\nabla\mu=\bm{0}$ , we infer from the first equation $\varphi=0$ , and consequently $\mu=0$ . Then, $q^{n}$ is uniquely determined from (2.2).

Remark 2.3.

( $\mathrm{A5}$ ) essentially requires that the potential $F$ has quartic polynomial growth at infinity, which is fulfilled by the classical example $F(s)=\frac{1}{4}(s^{2}-1)^{2}$ . We show in Remark 2.5 that this seems to be sharp for our main stability result.

Remark 2.4.

It is also possible to consider $g\in L^{2}(0,T;H^{1}(\Omega))$ as oppose to ( $\mathrm{A4}$ ), since we can define a new variable $\hat{\mu}:=\mu+g$ and shift $g$ from (1.1b) into the diffusion term of (1.1a), see Section 5.4 for an example.

Our main result is the following stability estimate for the time discrete solutions.

Theorem 2.1 (Stability).

Suppose, for any $\tau\in(0,1)$ , the time discrete system (2.1) is solvable. Then, under ( $\mathrm{A1}$ )-( $\mathrm{A5}$ ), there exists $\tau_{*}\in(0,1)$ depending only on model parameters, such that for all $\tau\in(0,\tau_{*})$ , the following estimate holds for the corresponding discrete solutions $\{\varphi^{k},\mu^{k},r^{k},q^{k}\}_{k=1}^{N_{\tau}}$ , where $\tau N_{\tau}\leq T$ , with a positive constant $C$ independent of $\tau$ and $N_{\tau}$ :

		$\displaystyle\\|\varphi^{k}\\|_{H^{1}}^{2}+\|q^{k}\|^{2}+\|r^{k}\|^{2}+\sum_{n=1}^{k}\Big{(}\\|\varphi^{n}-\varphi^{n-1}\\|_{H^{1}}^{2}+\|r^{n}-q^{n-1}\|^{2}+\|q^{n}-q^{n-1}\|^{2}\Big{)}$		(2.5)
		$\displaystyle\quad+\tau\sum_{n=1}^{k}\Big{(}\\|\mu^{n}\\|_{H^{1}}^{2}+\\|\varphi^{n}\\|_{H^{2}}^{4}\Big{)}+\sum_{n=1}^{k}\frac{1}{\tau}\\|\varphi^{n}-\varphi^{n-1}\\|_{(H^{1})^{*}}^{2}\leq C,$		(2.5)

and for any $h\in\{2,\dots,N_{\tau}\}$ ,

\displaystyle\tau\sum_{n=0}^{N_{\tau}-1}\|\varphi^{n+1}-\varphi^{n}\|^{2}\leq C\tau^{3/2},\quad\tau\sum_{n=0}^{N_{\tau}-h}\|\varphi^{n+h}-\varphi^{n}\|^{2}\leq Ch\tau.

(2.6)

Proof.

Let $\{\varphi^{k},\mu^{k},r^{k},q^{k}\}_{k=1}^{N_{\tau}}$ denote the discrete solutions to (2.1)-(2.2) originating from the initial conditions $\varphi^{0}:=\varphi_{0}$ and $q^{0}:=Q(\varphi_{0})$ . For $k=n$ , we test (2.1a) with $\mu^{n}$ and also with $\varphi^{n}$ , (2.1b) with $\varphi^{n}-\varphi^{n-1}$ and (2.1c) with $2r^{n}$ . Then, employing the identity $(a-b)a=\tfrac{1}{2}(a^{2}-b^{2}+(a-b)^{2})$ we first obtain upon summing

	$\displaystyle\frac{1}{2}\Big{(}\\|\varphi^{n}\\|^{2}+\varepsilon\\|\nabla\varphi^{n}\\|^{2}+2\|r^{n}\|^{2}\Big{)}-\frac{1}{2}\Big{(}\\|\varphi^{n-1}\\|^{2}+\varepsilon\\|\nabla\varphi^{n-1}\\|^{2}+2\|q^{n-1}\|^{2}\Big{)}$
	$\displaystyle\qquad+\frac{1}{2}\Big{(}\\|\varphi^{n}-\varphi^{n-1}\\|^{2}+\varepsilon\\|\nabla(\varphi^{n}-\varphi^{n-1})\\|^{2}\Big{)}+\|r^{n}-q^{n-1}\|^{2}+\int_{\Omega}\tau m(\varphi^{n-1})\|\nabla\mu^{n}\|^{2}\,\mathrm{dx}$
	$\displaystyle\quad=\tau(f^{n-1},\mu^{n})+\tau(f^{n-1},\varphi^{n})+(g^{n-1},\varphi^{n}-\varphi^{n-1})-\tau(m(\varphi^{n-1})\nabla\mu^{n},\nabla\varphi^{n}).$

Combining with the following inequality derived from (2.2) and (2.3)

|q^{n}|^{2}+\frac{1}{2}|q^{n}-q^{n-1}|^{2}\leq\tau M+\tau\eta m_{0}\|\nabla\mu^{n}\|^{2}+|r^{n}|^{2}+\frac{1}{2}|r^{n}-q^{n-1}|^{2},

we obtain the inequality

		$\displaystyle G^{n}-G^{n-1}+\frac{1}{2}\\|\varphi^{n}-\varphi^{n-1}\\|^{2}+\frac{\varepsilon}{2}\\|\nabla(\varphi^{n}-\varphi^{n-1})\\|^{2}+\frac{1}{2}\|q^{n}-q^{n-1}\|^{2}$		(2.7)
		$\displaystyle\qquad+\frac{1}{2}\|r^{n}-q^{n-1}\|^{2}+\int_{\Omega}(1-\eta)\tau m_{0}\|\nabla\mu^{n}\|^{2}\,\mathrm{dx}$
		$\displaystyle\quad\leq\tau(f^{n-1},\mu^{n})+\tau M+\tau(f^{n-1},\varphi^{n})+(g^{n-1},\varphi^{n}-\varphi^{n-1})$
		$\displaystyle\qquad-\tau(m(\varphi^{n-1})\nabla\mu^{n},\nabla\varphi^{n}),$

where

\displaystyle G^{k}:=\frac{1}{2}\|\varphi^{k}\|^{2}+\frac{\varepsilon}{2}\|\nabla\varphi^{k}\|^{2}+|q^{k}|^{2}

(2.8)

plays the role of a discrete energy functional. The latter three terms on the right-hand side of (2.7) can be handled in a standard way, where for the remainder of the proof, the symbol $C$ denotes a generic constant independent of $\tau$ and $n\in\{1,\dots,N_{\tau}\}$ whose value may change line from line and within the same line:

		$\displaystyle\tau(f^{n-1},\varphi^{n})+(g^{n-1},\varphi^{n}-\varphi^{n-1})-\tau(m(\varphi^{n-1})\nabla\mu^{n},\nabla\varphi^{n})$		(2.9)
		$\displaystyle\quad\leq C\tau+C\tau\\|\varphi^{n}\\|_{H^{1}}^{2}+C\\|g^{n-1}\\|^{2}+\frac{1}{4}\\|\varphi^{n}-\varphi^{n-1}\\|^{2}+\frac{\tau(1-\eta)m_{0}}{4}\\|\nabla\mu^{n}\\|^{2}.$		(2.9)

For the first term we use of the boundedness of $f$ and the Poincaré inequality to infer

	$\displaystyle\|(f^{n-1},\mu^{n})\|$	$\displaystyle\leq C\\|\mu^{n}\\|_{L^{1}}\leq C\\|\nabla\mu^{n}\\|_{L^{2}}+C\|(\mu^{n},1)\|$		(2.10)
		$\displaystyle\leq C+\frac{(1-\eta)m_{0}}{4}\\|\nabla\mu^{n}\\|^{2}+C\|(\mu^{n},1)\|.$		(2.10)

To close the estimate, it remains to estimate the mean value of $\mu^{n}$ , where from (2.1b),

\displaystyle C|(\mu^{n},1)|\leq C\frac{|r^{n}|}{\varepsilon Q(\varphi^{n-1})}|(F^{\prime}(\varphi^{n-1}),1)|+C|(g^{n-1},1)|.

(2.11)

By virtue of ( $\mathrm{A5}$ ), we have $|F^{\prime}(s)|\leq CF(s)^{\frac{3}{4}}+C$ and so

\displaystyle|(F^{\prime}(\varphi^{n-1}),1)|\leq C(F(\varphi^{n-1}),1)^{\frac{3}{4}}+C.

(2.12)

Combining with the facts $\varepsilon Q(\varphi^{n-1})^{2}=(F(\varphi^{n-1}),1)+\varepsilon C_{0}$ and $Q(\varphi^{n-1})\geq\sqrt{C_{0}}$ , and the Sobolev embedding $H^{1}(\Omega)\subset L^{4}(\Omega)$ , we deduce the following key estimate:

$\displaystyle C\frac{\|r^{n}\|\|(F^{\prime}(\varphi^{n-1}),1)\|}{\varepsilon Q(\varphi^{n-1})}$	$\displaystyle\leq C\frac{\|r^{n}\|((F(\varphi^{n-1}),1)^{\frac{3}{4}}+1)}{Q(\varphi^{n-1})}\leq C\|r^{n}\|\,(Q(\varphi^{n-1})^{\frac{1}{2}}+1)$	(2.13)
	$\displaystyle\leq\frac{1}{8}\|r^{n}\|^{2}+C\big{(}1+Q(\varphi^{n-1})\big{)}\leq\frac{1}{8}\|r^{n}\|^{2}+C\big{(}1+(F(\varphi^{n-1}),1)^{1/2}\big{)}$
	$\displaystyle\leq\frac{1}{8}\|r^{n}\|^{2}+C\big{(}1+\\|\varphi^{n-1}\\|_{L^{4}}^{2}\big{)}\leq\frac{1}{8}\|r^{n}\|^{2}+C\big{(}1+\\|\varphi^{n-1}\\|_{H^{1}}^{2}\big{)}.$

Substituting into (2.11) yields the following estimate for the mean value

	$\displaystyle C\|(\mu^{n},1)\|$	$\displaystyle\leq C+C\\|g^{n-1}\\|^{2}+\frac{1}{8}\|r^{n}\|^{2}+C\\|\varphi^{n-1}\\|_{H^{1}}^{2}$		(2.14)
		$\displaystyle\leq C+C\\|g^{n-1}\\|^{2}+\frac{1}{4}\|r^{n}-q^{n-1}\|^{2}+C\|q^{n-1}\|^{2}+C\\|\varphi^{n-1}\\|_{H^{1}}^{2},$		(2.14)

and consequently

	$\displaystyle\tau\|(f^{n-1},\mu^{n})\|$	$\displaystyle\leq C\tau+\frac{(1-\eta)\tau m_{0}}{4}\\|\nabla\mu^{n}\\|^{2}+\frac{\tau}{4}\|r^{n}-q^{n-1}\|^{2}$		(2.15)
		$\displaystyle\quad+C\tau\\|g^{n-1}\\|^{2}+C\tau\|q^{n-1}\|^{2}+C\tau\\|\varphi^{n-1}\\|_{H^{1}}^{2}.$		(2.15)

Using (2.9) and (2.15) in (2.7), we arrive at the following:

	$\displaystyle G^{n}+\frac{\min(1,2\varepsilon)}{4}\\|\varphi^{n}-\varphi^{n-1}\\|_{H^{1}}^{2}+\frac{1-\tau}{4}\|r^{n}-q^{n-1}\|^{2}+\frac{1}{2}\|q^{n}-q^{n-1}\|^{2}+\frac{m_{0}(1-\eta)\tau}{2}\\|\nabla\mu^{n}\\|^{2}$
	$\displaystyle\quad\leq G^{n-1}+C\tau(1+\\|g^{n-1}\\|^{2}+G^{n-1}+G^{n})+C\\|g^{n-1}\\|^{2}.$

For $\tau<1$ and arbitrary $k\in\{1,\dots,N_{\tau}\}$ , taking the sum of the above inequality from $n=1$ to $n=k$ yields

	$\displaystyle(1-C\tau)G^{k}+\frac{\min(1,2\varepsilon)}{4}\sum_{n=1}^{k}\\|\varphi^{n}-\varphi^{n-1}\\|_{H^{1}}^{2}+\frac{(1-\tau)}{4}\sum_{n=1}^{k}\|r^{n}-q^{n-1}\|^{2}$
	$\displaystyle\qquad+\frac{1}{2}\sum_{n=1}^{k}\|q^{n}-q^{n-1}\|^{2}+\frac{\tau(1-\eta)m_{0}}{2}\sum_{n=1}^{k}\\|\nabla\mu^{n}\\|^{2}$
	$\displaystyle\quad\leq G^{0}+C\sum_{n=1}^{k}\\|g^{n-1}\\|^{2}+C\tau\sum_{n=1}^{k}(1+\\|g^{n-1}\\|^{2}+G^{n-1}).$

Let $\tau_{*}$ be a fixed constant so that the prefactor $1-C\tau_{*}$ is positive, Then, for all $\tau\in(0,\tau_{*})$ , by the discrete Gronwall inequality it holds that for all $k=1,\dots,N_{\tau}$ ,

		$\displaystyle\\|\varphi^{k}\\|_{H^{1}}^{2}+\|q^{k}\|^{2}+\sum_{n=1}^{k}\\|\varphi^{n}-\varphi^{n-1}\\|_{H^{1}}^{2}+\sum_{n=1}^{k}\|r^{n}-q^{n-1}\|^{2}$		(2.16)
		$\displaystyle\quad+\sum_{n=1}^{k}\|q^{n}-q^{n-1}\|^{2}+\tau\sum_{n=1}^{k}\\|\nabla\mu^{n}\\|^{2}\leq C.$		(2.16)

In turn, we find that for all $k=1,\dots,N_{\tau}$ ,

\displaystyle|r^{k}|^{2}\leq 2|r^{k}-q^{k-1}|^{2}+2|q^{k-1}|^{2}\leq C.

(2.17)

Returning to the mean value estimate (2.14) and by the Poincaré inequality, we find

\displaystyle\tau\sum_{n=1}^{k}|(\mu^{n},1)|\leq C\quad\text{ and }\quad\tau\sum_{n=1}^{k}\|\mu^{n}\|^{2}\leq C.

(2.18)

Next, we test (2.1a) with an arbitrary $\zeta\in H^{1}(\Omega)$ , leading to

\frac{1}{\tau}|(\varphi^{n}-\varphi^{n-1},\zeta)|\leq m_{1}\|\nabla\mu^{n}\|\|\nabla\zeta\|+\|f^{n}\|\|\zeta\|\leq C\Big{(}1+\|\nabla\mu^{n}\|\Big{)}\|\zeta\|_{H^{1}}.

This shows that $\|\tfrac{1}{\tau}(\varphi^{n}-\varphi^{n-1})\|_{(H^{1})^{*}}\leq C(1+\|\nabla\mu^{n}\|)$ and hence

\sum_{n=1}^{k}\frac{1}{\tau}\|\varphi^{n}-\varphi^{n-1}\|_{(H^{1})^{*}}^{2}\leq C.

To show (2.6), we test (2.1a) with $\varphi^{m+h}-\varphi^{m}$ where $m=0,\dots,N_{\tau}-h$ and $h=1,\dots,N_{\tau}$ to obtain

0=(\varphi^{n}-\varphi^{n-1}-\tau f^{n-1},\varphi^{m+h}-\varphi^{m})+\tau(m(\varphi^{n-1})\nabla\mu^{n},\nabla(\varphi^{m+h}-\varphi^{m})).

Summing from $n=m+1,\dots,m+h$ yields

	$\displaystyle\\|\varphi^{m+h}-\varphi^{m}\\|^{2}$	$\displaystyle\leq C\tau\sum_{n=m+1}^{m+h}\Big{(}\\|f^{n-1}\\|+\\|\nabla\mu^{n}\\|\Big{)}\\|\varphi^{m+h}-\varphi^{m}\\|_{H^{1}}$
		$\displaystyle\leq C\tau\sum_{k=1}^{h}\Big{(}1+\\|\nabla\mu^{m+k}\\|\Big{)}\\|\varphi^{m+h}-\varphi^{m}\\|_{H^{1}}.$

Multiplying both sides by $\tau$ , summing from $m=0,\dots,N_{\tau}-h$ and applying Hölder’s inequality leads to

	$\displaystyle\tau\sum_{m=0}^{N_{\tau}-h}\\|\varphi^{m+h}-\varphi^{m}\\|^{2}\leq C\tau^{2}\sum_{k=1}^{h}\sum_{m=0}^{N_{\tau}-h}\Big{(}1+\\|\nabla\mu^{m+k}\\|\Big{)}\\|\varphi^{m+h}-\varphi^{m}\\|_{H^{1}}$
	$\displaystyle\quad\leq C\tau\sum_{k=1}^{h}\Big{(}\tau\sum_{m=0}^{N_{\tau}-h}\\|\varphi^{m+h}-\varphi^{m}\\|_{H^{1}}^{2}\Big{)}^{1/2}\Big{(}1+\Big{(}\tau\sum_{m=0}^{N_{\tau}-h}\\|\nabla\mu^{m+k}\\|^{2}\Big{)}^{1/2}\Big{)}$
	$\displaystyle\quad\leq C\tau\sum_{k=1}^{h}\Big{(}N_{\tau}\tau\max_{n=1,\dots,N_{\tau}}\\|\varphi^{n}\\|_{H^{1}}^{2}\Big{)}^{1/2}\Big{(}1+\Big{(}\tau\sum_{n=0}^{N_{\tau}}\\|\nabla\mu^{n}\\|^{2}\Big{)}^{1/2}\Big{)}\leq Ch\tau.$

When $h=1$ , a higher exponent can be inferred by using (2.16):

\tau\sum_{m=0}^{N_{\tau}-1}\|\varphi^{m+1}-\varphi^{m}\|^{2}\leq C\tau\Big{(}\tau\sum_{m=0}^{N_{\tau}-1}\|\varphi^{m+1}-\varphi^{m}\|_{H^{1}}^{2}\Big{)}^{1/2}\leq C\tau^{3/2}.

Furthermore, we test (2.1b) with $-\Delta\varphi^{n}$ , which yields

	$\displaystyle\frac{\varepsilon}{2}\\|\Delta\varphi^{n}\\|^{2}$	$\displaystyle\leq C\\|g^{n-1}\\|^{2}+C\frac{\|r^{n}\|^{2}}{Q(\varphi^{n-1})^{2}}\\|F^{\prime}(\varphi^{n-1})\\|^{2}$		(2.19)
		$\displaystyle\quad+C\\|\nabla\mu^{n}\\|\\|\nabla\varphi^{n}\\|\leq C+C\\|\nabla\mu^{n}\\|,$		(2.19)

after applying (2.16), (2.17), the lower bound $Q(\varphi^{n-1})\geq\sqrt{C_{0}}$ , and ( $\mathrm{A5}$ ). Invoking elliptic regularity (e.g. [20, Thm. 2.4.2.7]), we find that

\displaystyle\tau\sum_{n=1}^{k}\|\varphi^{n}\|_{H^{2}}^{4}\leq C\tau\sum_{n=1}^{k}\big{(}\|\Delta\varphi^{n}\|^{4}+\|\varphi^{n}\|^{4}\big{)}\leq C.

(2.20)

This completes the proof. ∎

Remark 2.5 (Sharpness of ( $\mathrm{A5}$ )).

If $F$ satisfies the growth assumptions

c|s|^{p}-C\leq F(s)\leq C(1+|s|^{p}),\quad|F^{\prime}(s)|\leq C(1+|s|^{p-1}),

for some $p\in[2,\infty)$ if $d=2$ , or $p\in[2,6)$ if $d=3$ , via a similar calculation as (2.13), we find that

	$\displaystyle\frac{\|r^{n}\|}{\varepsilon Q(\varphi^{n-1})}\|(F^{\prime}(\varphi^{n-1}),1)\|$	$\displaystyle\leq C\frac{\|r^{n}\|}{Q(\varphi^{n-1})}(F(\varphi^{n-1}),1)^{\frac{p-1}{p}}$
		$\displaystyle\quad\leq\frac{1}{8}\|r^{n}\|^{2}+C(F(\varphi^{n-1}),1)^{\frac{p-2}{p}}\leq\frac{1}{8}\|r^{n}\|^{2}+C\\|\varphi^{n-1}\\|_{L^{p}}^{p-2},$

where we have omitted lower order terms. Let $\alpha=\frac{p-2}{p}$ if $d=2$ , or $\alpha=\frac{3(p-2)}{2p}$ if $d=3$ , so that $\alpha\in(0,1)$ , and by invoking the Gagliardo–Nirenberg inequality,

\displaystyle\|\varphi^{n-1}\|_{L^{p}}^{p-2}

\displaystyle\leq C\|\varphi^{n-1}\|_{L^{2}}^{(1-\alpha)(p-2)}\|\varphi^{n-1}\|_{H^{1}}^{\alpha(p-2)}\leq C\|\varphi^{n-1}\|_{L^{2}}^{\frac{s(1-\alpha)(p-2)}{s-1}}+C\|\varphi^{n-1}\|_{H^{1}}^{\alpha(p-2)s},

for some $s>1$ . Hence, to control $\|\varphi^{n-1}\|_{L^{p}}^{p-2}$ with $\|\varphi^{n-1}\|_{H^{1}}^{2}$ , we demand

\alpha(p-2)s=\frac{s(1-\alpha)(p-2)}{s-1}=2,

and a short calculation shows that $p=4$ . Hence, to use the idea of (2.13) to proceed, the quartic growth assumption ( $\mathrm{A5}$ ) on $F$ seems to be sharp.

2.3 Discussion on the polynomial growth of $F$

The boundedness assumption ( $\mathrm{A5}$ ) on $f$ is essential for the above analysis to go through if one chooses $F$ to be a polynomial potential with super-quadratic growth. In particular, if $f$ is assumed to belong to $L^{2}(0,T;L^{2}(\Omega))$ , or treating $f=\hat{f}(\varphi)$ where $|\hat{f}(s)|\leq C(1+|s|)$ has linear growth, then the analysis restricts $F$ to have at most quadratic growth. Indeed, we revisit (2.10) where upon employing the Poincaré inequality:

|(f^{n-1},\mu^{n})|\leq C\|f^{n-1}\|\|\mu^{n}\|\leq C\|f^{n-1}\|\big{(}\|\nabla\mu^{n}\|+|(\mu^{n},1)|\big{)}.

Control of this term requires an a priori estimate on the square of the mean value in terms of $|r^{n}|^{2}$ and $\|\varphi^{n}\|_{H^{1}}^{2}$ . From (2.11) we see

|(\mu^{n},1)|^{2}\leq C\frac{|r^{n}|^{2}}{Q(\varphi^{n-1})^{2}}|(F^{\prime}(\varphi^{n-1}),1)|^{2}+C\|g^{n-1}\|^{2}.

If $|F^{\prime}(s)|\sim|s|^{p}$ for some exponent $p\geq 1$ , then applying a similar argument (omitting the lower order terms), we find that

\displaystyle\frac{|r^{n}|^{2}}{Q(\varphi^{n-1})^{2}}|(F^{\prime}(\varphi^{n-1}),1)|^{2}\leq C\frac{|r^{n}|^{2}}{Q(\varphi^{n-1})^{2}}|(F(\varphi^{n-1}),1)^{\frac{2p}{p+1}}\leq C|r^{n}|^{2}Q(\varphi^{n-1})^{2\frac{p-1}{p+1}}.

In order for this term to be controlled by the left-hand side of (2.7), we are forced to take $p=1$ , which limits $F$ to at most quadratic growth at infinity.

Remark 2.6.

We mention that solutions to models with quadratic potentials (e.g. as limits of a SLI-based numerical scheme) may exhibit behavior different from those corresponding to models with super-quadratic potentials. For instance, as noted in [30, Remark 2.5] for a Cahn–Hilliard tumor model (see (5.11) below), establishing the long-time behavior of weak solutions seem to demand $F$ to have a polynomial growth at infinity faster than cubic.

3 Convergence of the time discrete solutions

In this section we show a weak solution to (1.1) can be obtained in the limit $\tau\to 0$ . As the spatial discretization can be done in a standard way, we omit the details and refer to recent works [18, 28] for the convergence analysis of fully discrete solutions. For fixed $\tau$ , we introduce the affine linear and piecewise constant extensions of time discrete functions $a^{n}$ , $n=1,\dots,N_{\tau}$ :

	$\displaystyle a^{\tau}(t):=\frac{t-t^{n-1}}{\tau}a^{n}+\frac{t^{n}-t}{\tau}a^{n-1}$	$\displaystyle\text{ for }t\in[t^{n-1},t^{n}],$		(3.1)
	$\displaystyle a^{\tau,+}(t):=a^{n},\quad a^{\tau,-}(t):=a^{n-1}$	$\displaystyle\text{ for }t\in(t^{n-1},t^{n}],$		(3.1)

for $a\in\{\varphi,\mu,r,q\}$ , where we set $\varphi^{\tau}(0)=\varphi_{0}$ and $q^{\tau}(0)=Q(\varphi_{0})$ . Upon noting that

	$\displaystyle a^{\tau}(t)-a^{\tau,-}(t)=\frac{t-t^{n-1}}{\tau}(a^{n}-a^{n-1})$	$\displaystyle\text{ for }t\in(t^{n-1},t^{n}],$		(3.2)
	$\displaystyle a^{\tau,+}(t)-a^{\tau}(t)=\frac{t^{n}-t}{\tau}(a^{n}-a^{n-1})$	$\displaystyle\text{ for }t\in(t^{n-1},t^{n}],$		(3.2)

we deduce from Theorem 2.1 the following uniform estimates (where the notation $a^{\tau,(\pm)}$ is a shorthand for $\{a^{\tau},a^{\tau,+},a^{\tau,-}\}$ , while $a^{\tau,\pm}$ is a shorthand for $\{a^{\tau,+},a^{\tau,-}\}$ ):

		$\displaystyle\\|\varphi^{\tau,(\pm)}\\|_{L^{\infty}(0,T;H^{1})}^{2}+\\|\varphi^{\tau,(\pm)}\\|_{L^{4}(0,T;H^{2})}^{4}+\\|q^{\tau,(\pm)}\\|_{L^{\infty}(0,T)}^{2}+\\|r^{\tau,+}\\|_{L^{\infty}(0,T)}^{2}$		(3.3)
		$\displaystyle\qquad+\frac{1}{\tau}\\|r^{\tau,+}-q^{\tau,-}\\|_{L^{2}(0,T)}^{2}+\\|\mu^{\tau,+}\\|_{L^{2}(0,T;H^{1})}^{2}+\\|\partial_{t}\varphi^{\tau}\\|_{L^{2}(0,T;(H^{1})^{*})}^{2}$
		$\displaystyle\qquad+\frac{1}{\tau}\\|q^{\tau}-q^{\tau,\pm}\\|_{L^{2}(0,T)}^{2}+\frac{1}{\tau}\\|\varphi^{\tau}-\varphi^{\tau,\pm}\\|_{L^{2}(0,T;H^{1})}^{2}\leq C,$

with the last two terms coming from taking the $L^{2}$ -norm of (3.2) for $a\in\{\varphi,q\}$ . Furthermore, from (2.6), we infer that

	$\displaystyle\tau^{2}\\|\partial_{t}\varphi^{\tau}\\|_{L^{2}(0,T;L^{2})}^{2}=\\|\varphi^{\tau,+}-\varphi^{\tau,-}\\|_{L^{2}(0,T;L^{2})}^{2}$	$\displaystyle\leq C\tau^{\frac{3}{2}},$		(3.4)
	$\displaystyle\\|\varphi^{\tau,(\pm)}(\cdot+h\tau)-\varphi^{\tau,(\pm)}\\|_{L^{2}(0,T-h\tau;L^{2})}^{2}$	$\displaystyle\leq Ch\tau\quad\text{ for any }h\in\{1,\dots,N_{\tau}\}.$		(3.5)

Then, for arbitrary $\phi\in L^{2}(0,T;H^{1}(\Omega))$ testing (2.1a) and (2.1b) with $\phi$ , summing from $n=1$ to $N_{\tau}$ , and using the above definitions, we arrive at


$\displaystyle 0$	$\displaystyle=\int_{0}^{T}(\partial_{t}\varphi^{\tau}-f^{\tau},\phi)+(m(\varphi^{\tau,-})\nabla\mu^{\tau,+},\nabla\phi)\,\mathrm{dt},$	(3.6a)
$\displaystyle 0$	$\displaystyle=\int_{0}^{T}\varepsilon(\nabla\varphi^{\tau,+},\nabla\phi)+\frac{r^{\tau,+}}{\varepsilon Q(\varphi^{\tau,-})}(F^{\prime}(\varphi^{\tau,-}),\phi)-(\mu^{\tau,+}+g^{\tau},\phi)\,\mathrm{dt}.$	(3.6b)

The main difference compared to the convergence analysis of the standard SAV scheme [33] is that the left-hand side of (2.1c) does not translate well into a time derivative of $r^{\tau}$ or $q^{\tau}$ . However, as we note below, for strict RSAV schemes (i.e., $\zeta_{n}^{\tau}<1$ in (2.2)) it suffices to show that the limits of $r^{\tau,+}$ and $Q(\varphi^{\tau,-})$ coincide, and equation (2.1c) does not play a role.

In the following we explore the convergence analysis for two choices of $\zeta_{n}^{\tau}$ . As the choices are uniform in $n$ we use the notation $\zeta^{\tau}$ . The first considers $\zeta^{\tau}=\zeta$ with a fixed $\zeta\in[0,1]$ , covering both the standard SAV scheme $q^{n}=r^{n}$ with $\zeta=1$ , and the idealized update approach $q^{n}=Q(\varphi^{n})$ with $\zeta=0$ . Note that aside from $\zeta=1$ , it is unknown whether for other values $\zeta\in[0,1)$ the criterion (2.3) is fulfilled for all $n=1,\dots,N_{\tau}$ . This motivates the second choice where we consider $\zeta^{\tau}=1-e_{\tau}\tau$ , with a non-negative sequence $\{e_{\tau}\}_{\tau>0}$ . This results in a method that asymptotically close to the standard SAV scheme with the relaxation effect almost completely negated for small time steps. Recalling the discussion before Remark 2.1, we can pick $e_{\tau}=\max_{n}(1-w_{n}^{\tau})$ , so that $\zeta^{\tau}=1-e_{\tau}\tau\in\mathcal{V}_{n}^{\tau}$ for all $n=1,\dots,N_{\tau}$ .

Theorem 3.1 (Convergence).

For any $\tau\in(0,\tau_{*})$ , let $\{(\varphi^{k},\mu^{k},r^{k},q^{k})\}_{k=1}^{N_{\tau}}$ denote the time discrete solutions to the relaxed SAV scheme (2.1)-(2.2), where we assume that either a fixed $\zeta\in[0,1]$ belongs to $\mathcal{V}_{n}^{\tau}$ for all $n=1,\dots,N_{\tau}$ , or there exists a non-negative convergent sequence $\{e_{\tau}\}_{\tau>0}$ with $\lim_{\tau\to 0}e_{\tau}=e_{*}\in[0,\infty)$ such that $1-e_{\tau}\tau\in\mathcal{V}_{n}^{\tau}$ for all $n=1,\dots,N_{\tau}$ . Then, there exists a non-relabelled subsequence $\tau\to 0$ and limit functions

	$\displaystyle\varphi$	$\displaystyle\in L^{\infty}(0,T;H^{1}(\Omega))\cap H^{1}(0,T;(H^{1}(\Omega))^{*})\cap L^{4}(0,T;H^{2}(\Omega)),$
	$\displaystyle\mu$	$\displaystyle\in L^{2}(0,T;H^{1}(\Omega)),\quad q\in L^{\infty}(0,T),$

such that the interpolation functions $\{\varphi^{\tau,(\pm)},\mu^{\tau,+},r^{\tau,+},q^{\tau,(\pm)}\}$ satisfy the following compactness assertions for any $s<\infty$ in $d=2$ and $s<6$ in $d=3$ :

	$\displaystyle\varphi^{\tau,(\pm)}\to\varphi$	$\displaystyle\text{ weakly* in }L^{\infty}(0,T;H^{1}(\Omega))\cap L^{4}(0,T;H^{2}(\Omega)),$
	$\displaystyle\partial_{t}\varphi^{\tau}\to\partial_{t}\varphi$	$\displaystyle\text{ weakly in }L^{2}(0,T;(H^{1}(\Omega))^{*}),$
	$\displaystyle\varphi^{\tau,(\pm)}\to\varphi$	$\displaystyle\text{ strongly in }L^{2}(0,T;L^{s}(\Omega))\text{ and a.e.~{}in }\Omega\times(0,T),$
	$\displaystyle\mu^{\tau,+}\to\mu$	$\displaystyle\text{ weakly in }L^{2}(0,T;H^{1}(\Omega)),$
	$\displaystyle q^{\tau,(\pm)}\text{ and }r^{\tau,+}\to q$	$\displaystyle\text{ weakly* in }L^{\infty}(0,T),$
	$\displaystyle Q(\varphi^{\tau,(\pm)})\to Q(\varphi)$	$\displaystyle\text{ strongly in }L^{2}(0,T).$

The limit pair $(\varphi,\mu)$ is a weak solution to (1.1) in the sense


$\displaystyle\langle\partial_{t}\varphi,\phi\rangle+(m(\varphi)\nabla\mu,\nabla\phi)-(f,\phi)$	$\displaystyle=0,$	(3.7a)
$\displaystyle\varepsilon(\nabla\varphi,\nabla\phi)+\varepsilon^{-1}(F^{\prime}(\varphi),\phi)-(\mu+g,\phi)$	$\displaystyle=0,$	(3.7b)

for all $\phi\in H^{1}(\Omega)$ and for a.e. $t\in(0,T)$ , while $q(t)=Q(\varphi(t))$ holds for a.e. $t\in(0,T)$ .

Proof.

The weak/weak* convergences of $\varphi^{\tau,(\pm)}$ , $\partial_{t}\varphi^{\tau}$ , $\mu^{\tau,+}$ , $q^{\tau,(\pm)}$ and $r^{\tau,+}$ follow from standard compactness results in Bochner spaces. Then, from (3.3) we infer

\displaystyle\|r^{\tau,+}-q^{\tau,-}\|_{L^{2}(0,T)}^{2}+\|q^{\tau}-q^{\tau,\pm}\|_{L^{2}(0,T)}^{2}+\|\varphi^{\tau}-\varphi^{\tau,\pm}\|_{L^{2}(0,T;H^{1})}^{2}\leq C\tau,

(3.8)

which yields the uniqueness of weak limits for $\{\varphi^{\tau,(\pm)}\}$ and for $\{q^{\tau,(\pm)},r^{\tau,+}\}$ . The strong convergence of $\varphi^{\tau,(\pm)}$ in $L^{2}(0,T;L^{s}(\Omega))$ comes from the application of [39, §8, Thm. 5] with the uniform boundedness of $\varphi^{\tau,(\pm)}$ in $L^{\infty}(0,T;H^{1}(\Omega))$ , the compact embedding $H^{1}(\Omega)\subset\subset L^{s}(\Omega)$ and the uniform time translation estimate (3.5). Let $F^{(0)}(s)=F(s)$ and $F^{(1)}(s)=F^{\prime}(s)$ . Then, a short calculation with the estimate

|F^{(i)}(s_{1})-F^{(i)}(s_{2})|\leq C\big{(}1+|s_{1}|^{3-i}+|s_{2}|^{3-i}\big{)}|s_{1}-s_{2}|\text{ for }i=0,1,

together with the boundedness of $\varphi^{\tau,(\pm)},\varphi$ in $L^{4}(0,T;H^{2}(\Omega))\cap L^{\infty}(0,T;H^{1}(\Omega))$ yields $F^{\prime}(\varphi^{\tau,(\pm)})\in L^{\infty}(0,T;L^{2}(\Omega))$ and

	$\displaystyle\\|F^{\prime}(\varphi^{\tau,(\pm)})-F^{\prime}(\varphi)\\|_{L^{\frac{6}{5}}}$	$\displaystyle\leq C\big{(}1+\\|\varphi^{\tau,(\pm)}\\|_{L^{6}}^{2}+\\|\varphi\\|_{L^{6}}^{2}\big{)}\\|\varphi^{\tau,(\pm)}-\varphi\\|\leq C\\|\varphi^{\tau,(\pm)}-\varphi\\|,$		(3.9)
	$\displaystyle\\|F(\varphi^{\tau,(\pm)})-F(\varphi)\\|_{L^{1}}$	$\displaystyle\leq C\big{(}1+\\|\varphi^{\tau,(\pm)}\\|_{H^{1}}^{3}+\\|\varphi\\|_{H^{1}}^{3}\big{)}\\|\varphi^{\tau,(\pm)}-\varphi\\|\leq C\\|\varphi^{\tau,(\pm)}-\varphi\\|.$		(3.9)

Together with the following estimate

|Q(\varphi^{\tau,(\pm)})-Q(\varphi)|\leq\frac{1}{\varepsilon}\frac{\|F(\varphi^{\tau,(\pm)})-F(\varphi)\|_{L^{1}}}{|Q(\varphi^{\tau,(\pm)})+Q(\varphi)|}\leq C\|\varphi^{\tau,(\pm)}-\varphi\|,

we deduce the strong convergence $Q(\varphi^{\tau,(\pm)})\to Q(\varphi)$ in $L^{2}(0,T)$ .

Then, using the boundeness of $m$ and the a.e. convergence of $\varphi^{\tau,-}$ we can infer that $m(\varphi^{\tau,-})\nabla\zeta\to m(\varphi)\nabla\zeta$ strongly in $L^{2}(0,T;L^{2}(\Omega))$ . Thus, passing to the limit $\tau\to 0$ in (3.6a) leads to (3.7a), and by the continuous embedding $L^{2}(0,T;H^{1}(\Omega))\cap H^{1}(0,T;(H^{1}(\Omega))^{*})\subset C^{0}([0,T];L^{2}(\Omega))$ , we have the attainment of the initial condition $\varphi(0)=\varphi_{0}$ . In the case where $f=\hat{f}(\varphi)$ is a bounded continuous function of $\varphi$ , we replace $f^{\tau}$ in (3.6a) with $f(\varphi^{\tau,-})$ . Then, invoking the a.e. convergence of $\varphi^{\tau,-}$ in $\Omega\times(0,T)$ and the dominated convergence theorem yield that $\int_{0}^{T}(f(\varphi^{\tau,-}),\phi)\,\mathrm{dt}\to\int_{0}^{T}(f(\varphi),\phi)\,\mathrm{dt}$ .

Passing to the limit $\tau\to 0$ in the terms of (3.6b) other than the potential term is straightforward, and so we omit the details. For arbitrary $\phi\in L^{2}(0,T;H^{1}(\Omega))$ we have $\frac{(F^{\prime}(\varphi),\phi)}{Q(\varphi)}\in L^{1}(0,T)$ and thus by (3.9)

		$\displaystyle\left\|\int_{0}^{T}\frac{r^{\tau,+}(F^{\prime}(\varphi^{\tau,-}),\phi)}{\varepsilon Q(\varphi^{\tau,-})}\mp\frac{r^{\tau,+}(F^{\prime}(\varphi^{\tau,-}),\phi)}{\varepsilon Q(\varphi)}\mp\frac{r^{\tau,+}(F^{\prime}(\varphi),\phi)}{\varepsilon Q(\varphi)}-\frac{q(F^{\prime}(\varphi),\phi)}{\varepsilon Q(\varphi)}\,\mathrm{dt}\right\|$		(3.10)
		$\displaystyle\quad\leq\int_{0}^{T}\frac{\|(F^{\prime}(\varphi^{\tau,-}),\phi)\|}{\varepsilon Q(\varphi^{\tau,-})Q(\varphi)}\frac{\|r^{\tau,+}\|}{Q(\varphi^{\tau,-})+Q(\varphi)}\\|F(\varphi^{\tau,-})-F(\varphi)\\|_{L^{1}}\,\mathrm{dt}$
		$\displaystyle\qquad+\int_{0}^{T}\frac{\|r^{\tau,+}\|}{\varepsilon Q(\varphi)}\\|\phi\\|_{L^{6}}\\|F^{\prime}(\varphi^{\tau,-})-F^{\prime}(\varphi)\\|_{L^{\frac{6}{5}}}\,\mathrm{dt}+\left\|\int_{0}^{T}(r^{\tau,+}-q)\frac{(F^{\prime}(\varphi),\phi)}{\varepsilon Q(\varphi)}\,\mathrm{dt}\right\|$
		$\displaystyle\quad\leq C\\|\varphi^{\tau,-}-\varphi\\|_{L^{2}(0,T;L^{2})}\\|\phi\\|_{L^{2}(0,T;H^{1})}+\Big{\|}\int_{0}^{T}(r^{\tau,+}-q)\frac{(F^{\prime}(\varphi),\phi)}{\varepsilon Q(\varphi)}\,\mathrm{dt}\Big{\|}\to 0,$

as $\tau\to 0$ , where we have applied the lower bound on $Q$ . This implies that

\displaystyle\lim_{\tau\to 0}\int_{0}^{T}\frac{r^{\tau,+}}{\varepsilon Q(\varphi^{\tau,-})}(F^{\prime}(\varphi^{\tau,-}),\phi)\,\mathrm{dt}=\int_{0}^{T}\frac{q}{\varepsilon Q(\varphi)}(F^{\prime}(\varphi),\phi)\,\mathrm{dt}.

(3.11)

It now remains to show $q=Q(\varphi)$ , so that we recover (3.7b) completely. We first consider the case $\zeta^{\tau}=\zeta\in[0,1)$ , i.e., for strict RSAV schemes. With the strong convergence $Q(\varphi^{\tau,+})\to Q(\varphi)$ in $L^{2}(0,T)$ and the weak $*$ convergence $r^{\tau,+}\to q$ in $L^{\infty}(0,T)$ , we see that

q^{\tau,+}=\zeta r^{\tau,+}+(1-\zeta)Q(\varphi^{\tau,+})\to\zeta q+(1-\zeta)Q(\varphi)\text{ weakly in }L^{2}(0,T).

On the other hand, the weak $*$ convergence $q^{\tau,+}\to q$ in $L^{\infty}(0,T)$ yields the identity

q=\zeta q+(1-\zeta)Q(\varphi),

which implies $q=Q(\varphi)$ as $1-\zeta>0$ . Then, we recover (3.7b) and the limit pair $(\varphi,\mu)$ constitutes a weak solution to (1.1).

We treat the remaining cases $\zeta^{\tau}=1$ and $\zeta^{\tau}=1-e_{\tau}\tau$ together, where the former can be obtained from the latter by setting $e_{\tau}=e_{*}=0$ . Adding and subtracting $q^{n}$ to (2.1c), dividing by $\tau$ and testing with an arbitrary test function $\kappa\in H^{1}(0,T)$ such that $\kappa(T)=0$ . This yields after an integration by parts:

		$\displaystyle\int_{0}^{T}(Q(\varphi_{0})-q^{\tau})\kappa^{\prime}\,\mathrm{dt}-\int_{0}^{T}\kappa\frac{(F^{\prime}(\varphi^{\tau,-}),\partial_{t}\varphi^{\tau})}{2\varepsilon Q(\varphi^{\tau,-})}\,\mathrm{dt}$
		$\displaystyle\quad=\int_{0}^{T}\kappa\frac{q^{\tau,+}-r^{\tau,+}}{\tau}\,\mathrm{dt}\,=e_{\tau}\int_{0}^{T}\kappa(Q(\varphi^{\tau,+})-r^{\tau,+})\,\mathrm{dt},$		(3.12)

where we used that $q^{\tau}(0)=Q(\varphi_{0})$ . Arguing as in [28, Proof of Thm. 4.3], we obtain

$\displaystyle\int_{0}^{T}\kappa\frac{(F^{\prime}(\varphi^{\tau,-}),\partial_{t}\varphi^{\tau})}{2\varepsilon Q(\varphi^{\tau,-})}\,\mathrm{dt}$	$\displaystyle=\int_{0}^{T}\frac{\kappa(F^{\prime}(\varphi^{\tau,-}),\partial_{t}\varphi^{\tau})}{2\varepsilon Q(\varphi^{\tau,-})Q(\varphi^{\tau})}\Big{(}\frac{(F(\varphi^{\tau})-F(\varphi^{\tau,-}),1)}{Q(\varphi^{\tau,-})+Q(\varphi^{\tau})}\Big{)}\,\mathrm{dt}$	(3.13)
	$\displaystyle\quad+\int_{0}^{T}\kappa\frac{(F^{\prime}(\varphi^{\tau,-})-F^{\prime}(\varphi^{\tau}),\partial_{t}\varphi^{\tau})}{2\varepsilon Q(\varphi^{\tau})}\,\mathrm{dt}$
	$\displaystyle\quad+\int_{0}^{T}\kappa\frac{d}{dt}Q(\varphi^{\tau})\,\mathrm{dt}.$

Using the strong convergence of $Q(\varphi^{\tau})$ in $L^{2}(0,T)$ , the third term on the right-hand side of (3.13) can be treated as follows after an integration by parts:

\int_{0}^{T}\kappa\frac{d}{dt}Q(\varphi^{\tau})\,\mathrm{dt}=-\int_{0}^{T}\kappa^{\prime}Q(\varphi^{\tau})\,\mathrm{dt}-\kappa(0)Q(\varphi^{\tau}(0))\to-\int_{0}^{T}\kappa^{\prime}Q(\varphi)\,\mathrm{dt}-\kappa(0)Q(\varphi_{0}).

With the help of (3.4) and (3.8), the second term on the right-hand side of (3.13) can be estimated as

	$\displaystyle\left\|\int_{0}^{T}\kappa\frac{(F^{\prime}(\varphi^{\tau,-})-F^{\prime}(\varphi^{\tau}),\partial_{t}\varphi^{\tau})}{2\varepsilon Q(\varphi^{\tau})}\,\mathrm{dt}\right\|\leq C\\|\kappa\\|_{L^{\infty}(0,T)}\int_{0}^{T}\\|F^{\prime}(\varphi^{\tau,-})-F^{\prime}(\varphi^{\tau})\\|\\|\partial_{t}\varphi^{\tau}\\|\,\mathrm{dt}$
	$\displaystyle\quad\leq C\int_{0}^{T}\Big{(}1+\\|\varphi^{\tau,-}\\|_{L^{6}}^{2}+\\|\varphi^{\tau}\\|_{L^{6}}^{2}\Big{)}\\|\varphi^{\tau,-}-\varphi^{\tau}\\|_{L^{6}}\\|\partial_{t}\varphi^{\tau}\\|\,\mathrm{dt}$
	$\displaystyle\quad\leq C\\|\varphi^{\tau,-}-\varphi^{\tau}\\|_{L^{2}(0,T;H^{1})}\\|\partial_{t}\varphi^{\tau}\\|_{L^{2}(0,T;L^{2})}\leq C\tau^{\frac{1}{4}}\to 0\quad\text{ as }\tau\to 0.$

Similarly, the first term on the right-hand side of (3.13) can be estimated as

	$\displaystyle\left\|\int_{0}^{T}\frac{\kappa(F^{\prime}(\varphi^{\tau,-}),\partial_{t}\varphi^{\tau})}{2\varepsilon Q(\varphi^{\tau,-})Q(\varphi^{\tau})}\Big{(}\frac{(F(\varphi^{\tau})-F(\varphi^{\tau,-}),1)}{Q(\varphi^{\tau,-})+Q(\varphi^{\tau})}\Big{)}\,\mathrm{dt}\right\|$
	$\displaystyle\quad\leq C\\|\kappa\\|_{L^{\infty}(0,T)}\\|F^{\prime}(\varphi^{\tau,-})\\|_{L^{\infty}(0,T;L^{2})}\\|\partial_{t}\varphi^{\tau}\\|_{L^{2}(0,T;L^{2})}\\|F(\varphi^{\tau,-})-F(\varphi^{\tau})\\|_{L^{2}(0,T;L^{1})}$
	$\displaystyle\quad\leq C\\|\partial_{t}\varphi^{\tau}\\|_{L^{2}(0,T;L^{2})}\\|\varphi^{\tau,-}-\varphi^{\tau}\\|_{L^{2}(0,T;L^{2})}\leq C\tau^{\frac{1}{4}}\to 0\quad\text{ as }\tau\to 0.$

Then, passing to the limit in (3) yields

\int_{0}^{T}(\kappa^{\prime}-e_{*}\kappa)(Q(\varphi)-q)\,\mathrm{dt}=0

holding for arbitrary $\kappa\in H^{1}(0,T)$ with $\kappa(T)=0$ . Choosing $\kappa$ as the solution to the ordinary differential equation $\kappa^{\prime}-e_{*}\kappa=Q(\varphi)-q$ with terminal condition $\kappa(T)=0$ yields that

\|Q(\varphi)-q\|_{L^{2}(0,T)}=0,

which provides the identification $q=Q(\varphi)$ . Hence, in the case $\zeta^{\tau}=1-e_{\tau}\tau$ we also recover (3.7b) with the limit pair $(\varphi,\mu)$ constituting a weak solution to (1.1). This completes the proof. ∎

4 Numerical discussions

4.1 Fully discrete finite element approximation

For a bounded, convex domain $\Omega\subset\mathbb{R}^{d}$ , $d\in\{1,2,3\}$ , let $\{\mathcal{T}_{h}\}_{h>0}$ denote a regular family of conformal quasiuniform triangulations that partition $\Omega$ into disjoint open simplices $K$ such that $\max_{K\in\mathcal{T}_{h}}\text{diam}(K)\leq h$ . Let $\mathcal{S}_{h}$ denote the finite element space of continuous and piecewise linear functions:

\mathcal{S}_{h}:=\big{\{}q_{h}\in C^{0}(\overline{\Omega})\,:\,q_{h}|_{K}\in\mathcal{P}_{1}\;\forall K\in\mathcal{T}_{h}\big{\}}\subset H^{1}(\Omega),\quad U_{h}:=\text{dim}(\mathcal{S}_{h}),

with the set of basis functions $\{\chi_{h,k}\}_{k=1}^{U_{h}}$ that forms a dual basis to the set of vertices $\{\bm{x}_{k}\}_{k=1}^{U_{h}}$ of $\mathcal{T}_{h}$ . The nodal interpolation operator $\mathrm{I}_{h}:C^{0}(\overline{\Omega})\to\mathcal{S}_{h}$ is defined by $(\mathrm{I}_{h}\eta)(\bm{x}_{k})=\eta(\bm{x}_{k})$ for all $k=1,\dots,U_{h}$ , and we introduce the semi-inner product on $C^{0}(\overline{\Omega})$ :

(\eta_{1},\eta_{2})^{h}:=\int_{\Omega}\mathrm{I}_{h}(\eta_{1}\eta_{2})\,\mathrm{dx}=\sum_{k=1}^{U_{h}}\eta_{1}(\bm{x}_{k})\eta_{2}(\bm{x}_{k})\int_{\Omega}\chi_{h,k}(\bm{x})\,\mathrm{dx}.

If $f=f(t,x)$ and $g=g(t,x)$ are given functions, they are approximated with $f^{k}_{h}(\cdot)=\mathrm{I}_{h}^{c}(f^{k}(\cdot))$ and $g^{k}_{h}(\cdot)=\mathrm{I}_{h}^{c}(g^{k}(\cdot))$ for $k=1,\dots,N_{\tau}$ , via the Clément operator $\mathrm{I}_{h}^{c}:L^{2}(\Omega)\to\mathcal{S}_{h}$ [8]. If they are functions of $\varphi$ , we take $f_{h}^{k}=f(\varphi_{h}^{k})$ and $g_{h}^{k}=g(\varphi_{h}^{k})$ .

For more regular initial data $\varphi_{0}\in H^{2}(\Omega)$ we consider the discrete initial data $\varphi^{0}_{h}:=\mathrm{I}_{h}(\varphi_{0})\in\mathcal{S}_{h}$ , otherwise we take $\varphi^{0}_{h}=\mathrm{I}_{h}^{c}(\varphi_{0})\in\mathcal{S}_{h}$ , with $q^{0}_{h}:=Q_{h}(\varphi^{0}_{h})=((\varepsilon^{-1}F(\varphi^{0}_{h}),1)^{h}+C_{0})^{1/2}$ . The fully discrete finite element scheme for (1.4) reads as follows: Given data $\varphi^{0}_{h}\in\mathcal{S}_{h}$ , $q^{0}_{h}\in\mathbb{R}$ , and $\{f^{k}_{h},g^{k}_{h}\}_{k=1}^{N_{\tau}}$ , for $n=1,\dots,N_{\tau}$ , find discrete solutions $\varphi^{n}_{h},\mu^{n}_{h}\in\mathcal{S}_{h}$ and $r^{n}_{h}\in\mathbb{R}$ that satisfy for all $j=1,\dots,U_{h}$ :


$\displaystyle 0$	$\displaystyle=(\varphi^{n}_{h}-\varphi^{n-1}_{h},\chi_{h,j})^{h}+\tau(\mathrm{I}_{h}(m(\varphi^{n-1}_{h}))\nabla\mu^{n}_{h},\nabla\chi_{h,j})-\tau(f^{n-1}_{h},\chi_{h,j})^{h},$	(4.1a)
$\displaystyle 0$	$\displaystyle=\varepsilon(\nabla\varphi^{n}_{h},\nabla\chi_{h,j})-(\mu^{n}+g^{n-1}_{h},\chi_{h,j})^{h}+\frac{r^{n}_{h}}{\varepsilon Q_{h}(\varphi^{n-1}_{h})}(F^{\prime}(\varphi^{n-1}_{h}),\chi_{h,j})^{h},$	(4.1b)
$\displaystyle 0$	$\displaystyle=r^{n}_{h}-q^{n-1}_{h}-\frac{1}{2\varepsilon Q_{h}(\varphi^{n-1}_{h})}(F^{\prime}(\varphi^{n-1}_{h}),\varphi^{n}_{h}-\varphi^{n-1}_{h})^{h},$	(4.1c)

and update $q^{n}_{h}=\zeta_{n}^{\tau}r_{h}^{n}+(1-\zeta_{n}^{\tau})Q_{h}(\varphi^{n}_{h})$ with a constant $\zeta_{n}^{\tau}\in[0,1]$ such that the analogue of (2.3) is fulfilled. This fully discrete scheme is linear with respect to the unknowns $(\varphi^{n}_{h},\mu^{n}_{h},r^{n}_{h},q^{n}_{h})\in\mathcal{S}_{h}\times\mathcal{S}_{h}\times\mathbb{R}\times\mathbb{R}$ . Let us define the following matrices

\mathbb{M}_{ij}=(\chi_{h,i},\chi_{h,j})^{h},\quad\mathbb{S}_{ij}=(\nabla\chi_{h,i},\nabla\chi_{h,j}),\quad(\mathbb{S}_{h}^{n-1})_{ij}=(\mathrm{I}_{h}(m(\varphi_{h}^{n-1}))\nabla\chi_{h,i},\nabla\chi_{h,j}),

where we note that $\mathbb{M}$ is diagonal with $\mathbb{M}_{ii}=(\chi_{h,i},1)$ ; vectors $\bm{\varphi}^{n}_{h}:=\mathbb{M}^{-1}[(\varphi^{n}_{h},\chi_{h,k})^{h}]_{k=1}^{U_{h}}=(\varphi_{h}^{n}(\bm{x}_{1}),\varphi_{h}^{n}(\bm{x}_{2}),\dots,\varphi_{h}^{n}(\bm{x}_{U_{h}}))^{\top}$ , likewise for $\bm{\mu}^{n}_{h}$ , $\bm{\varphi}^{n-1}_{h}$ ; vectors $\bm{f}^{n-1}_{h}:=\mathbb{M}^{-1}[(f^{n-1}_{h},\chi_{h,k})^{h}]_{k=1}^{U_{h}}=(f_{h}^{n-1}(\bm{x}_{1}),f_{h}^{n-1}(\bm{x}_{2}),\dots,f_{h}^{n-1}(\bm{x}_{U_{h}}))^{\top}$ , likewise for $\bm{g}^{n-1}_{h}$ ; and vector

\bm{b}^{n-1}:=\frac{F^{\prime}(\bm{\varphi}^{n-1}_{h})}{\varepsilon Q_{h}(\varphi^{n-1}_{h})},\text{ where }Q_{h}(\varphi^{n-1}_{h})=\Big{(}\frac{1}{\varepsilon}(F(\bm{\varphi}^{n-1}_{h}),1)^{h}+C_{0}\Big{)}^{1/2},

with the nonlinearities applied component-wise. Then, (4.1) can be expressed as

	$\displaystyle\bm{0}$	$\displaystyle=\mathbb{M}(\bm{\varphi}^{n}_{h}-\bm{\varphi}^{n-1}_{h}-\tau\bm{f}^{n-1}_{h})+\tau\mathbb{S}^{n-1}_{h}(\varepsilon\mathbb{M}^{-1}\mathbb{S}\bm{\varphi}^{n}_{h}+r^{n}\bm{b}^{n-1}-\bm{g}^{n-1}_{h}),$
	$\displaystyle 0$	$\displaystyle=r^{n}_{h}-q^{n-1}_{h}-\frac{1}{2}\mathbb{M}\bm{b}^{n-1}\cdot(\bm{\varphi}^{n}_{h}-\bm{\varphi}^{n-1}_{h}).$

Defining the invertible matrix $\mathbb{B}_{n-1}:=\mathbb{I}+\varepsilon\tau\mathbb{M}^{-1}\mathbb{S}^{n-1}_{h}\mathbb{M}^{-1}\mathbb{S}$ with identity matrix $\mathbb{I}$ , substituting the second equation into the first equation, applying $\mathbb{B}_{n-1}^{-1}$ and taking the product with $\mathbb{M}\bm{b}^{n-1}$ would first yield an expression for $d^{n-1}=\,\mathbb{M}\bm{b}^{n-1}\cdot\bm{\varphi}^{n}_{h}$ as

	$\displaystyle d^{n-1}$	$\displaystyle:=\mathbb{M}\bm{b}^{n-1}\cdot\bm{\varphi}^{n}_{h}=\frac{\mathbb{M}\bm{b}^{n-1}\cdot\mathbb{B}_{n-1}^{-1}\bm{c}^{n-1}}{1+\frac{\tau}{2}\mathbb{M}\bm{b}^{n-1}\cdot(\mathbb{B}_{n-1}^{-1}\mathbb{M}^{-1}\mathbb{S}^{n-1}_{h}\bm{b}^{n-1})},$
	$\displaystyle\bm{c}^{n-1}$	$\displaystyle:=\bm{\varphi}^{n-1}_{h}+\tau\bm{f}^{n-1}_{h}+\tau\mathbb{M}^{-1}\mathbb{S}^{n-1}_{h}(\bm{g}^{n-1}_{h}+[\tfrac{1}{2}\mathbb{M}\bm{b}^{n-1}\cdot\bm{\varphi}^{n-1}_{h}-q^{n-1}_{h}]\bm{b}^{n-1}).$

Then, $\bm{\varphi}^{n}_{h}$ , $r^{n}_{h}$ , $\bm{\mu}^{n}_{h}$ and $q^{n}_{h}$ can be computed via

	$\displaystyle\bm{\varphi}^{n}_{h}$	$\displaystyle=\mathbb{B}_{n-1}^{-1}\bm{c}^{n-1}-\frac{\tau}{2}d^{n-1}\mathbb{B}_{n-1}^{-1}\mathbb{M}^{-1}\mathbb{S}^{n-1}_{h}\bm{b}^{n-1},\quad$		$\displaystyle r^{n}_{h}=q^{n-1}_{h}+\frac{1}{2}\mathbb{M}\bm{b}^{n-1}\cdot(\bm{\varphi}^{n}_{h}-\bm{\varphi}^{n-1}_{h}),$
	$\displaystyle\bm{\mu}^{n}_{h}$	$\displaystyle=\varepsilon\mathbb{M}^{-1}\mathbb{S}\bm{\varphi}^{n}_{h}-\bm{g}^{n-1}_{h}+r^{n}_{h}\bm{b}^{n-1},\quad$		$\displaystyle q^{n}_{h}=\zeta_{n}^{\tau}r^{n}_{h}+(1-\zeta_{n}^{\tau})Q_{h}(\varphi^{n}_{h}).$

As noted in [34], the main computational cost in each step amounts to solving two linear systems involving $\mathbb{B}_{n-1}$ . Once $\mathbb{B}_{n-1}^{-1}\mathbb{M}^{-1}\mathbb{S}_{h}^{n-1}\bm{b}^{n-1}$ and $\mathbb{B}_{n-1}^{-1}\bm{c}^{n-1}$ are computed, the constant $d^{n-1}$ and update $\bm{\varphi}_{h}^{n}$ do not involve any matrix inverse operations.

4.2 Convergence test

In this section we report on numerical tests for the scheme (2.1) on an interval $\Omega=[0,1]$ discretized with equidistant nodal points and spatial step size $h=0.01$ . Let us first explore the optimal values of $\zeta_{n}^{\tau}$ as defined in (2.4) for various values of $\eta$ and $M$ . We choose the source function $f$ so that the analytical solution to (1.1) with $m(\varphi)=1$ , $g=0$ , $F^{\prime}(\varphi)=\varphi^{3}-\varphi$ and $\varepsilon=1$ is

\varphi^{*}(x,t)=\cos(\pi x)(1+t).

We fix $C_{0}=1$ , $\tau=0.01$ and $N_{\tau}=500$ (so that $T=N_{\tau}\tau=5$ ). In our experiments we observe that the discrete energy $\frac{1}{2}\mathbb{S}\bm{\varphi}^{n}_{h}\cdot\bm{\varphi}^{n}_{h}+Q(\bm{\varphi}_{h}^{n})^{2}-C_{0}$ is increasing in time, see Figure 1(a), and thus the equation does not exhibit a dissipative structure involving the Ginzburg–Landau functional.

Refer to caption — (a) $\varphi^{*}(x,t)=\cos(\pi x)(1+t)$

For $M,\eta\sim O(10^{-3})$ , the optimal $\zeta_{n}^{\tau}$ is $0$ for all $n=1,\dots,N_{\tau}$ . Moreover, we observed a switching behavior of the optimal $\zeta_{n}^{\tau}$ when either one of the parameter is zero, and the other parameter is small. For example, with $\eta=0$ the optimal $\zeta_{n}^{\tau}$ is $0$ for all $n=1,\dots,N_{\tau}$ as long as $M\geq 0.05$ . Decreasing $M$ yields a switching of the optimal value to $\zeta_{n}^{\tau}=1$ for later iterations, and the smaller the value of $M$ the earlier the switching, see Figure 2(a) and (b). On the other hand, with $M=0$ and $\eta=10^{-5}$ the optimal $\zeta_{n}^{\tau}$ is zero except for the first iteration, similar to [43, Fig. 2], while with $\eta=10^{-6}$ the optimal $\zeta_{n}^{\tau}$ hovers around 0.9 for $\sim 330$ iterations before switching to $0$ , see Figure 2(c).

From our experiments it is interesting to see that the optimal $\zeta_{n}^{\tau}$ is 0 for moderate values of $M$ and $\eta$ , meaning that it is possible to consider the idealized update $q^{n}=Q(\varphi^{n})$ which preserves the consistency between the modified and original Ginzburg–Landau energy. It also suggests some form of numerical stability is available for the choice $q^{n}=Q(\varphi^{n})$ , although this is not immediate from our stability analysis. When $\eta$ and $M$ are sufficiently small or even zero, the tendency for the optimal $\zeta_{n}^{\tau}$ to be $1$ is supported by the fact that $R_{n}(1)\leq 0$ is satisfied even when $\eta=M=0$ .

Next we report on the $L^{2}(0,5;L^{2}(\Omega))$ -error between the numerical solution and the analytical solution. The switching behavior for the optimal $\zeta_{n}^{\tau}$ for $\tau=0.01$ displayed in Figure 2 does not occur for smaller time steps with the same values of $\eta$ and $M$ . I.e., for $\tau=0.0001$ , $M=0$ and $\eta=10^{-6}$ the optimal $\zeta_{n}^{\tau}$ is $0$ for all $n=1,\dots,N_{\tau}$ . Hence, we compare the numerical errors in the $L^{2}(0,5;L^{2}(\Omega))$ -norm with fixed $\zeta_{n}^{\tau}=\zeta\in\{0,0.25,0.5,0.75,1\}$ and a finer spatial discretization with step size $h=0.001$ . Since different time steps $\tau\in\{0.1,0.05,0.025,0.0125,0.00625,0.003125,0.0015625\}$ are used, with a terminal time $T=5$ the number of iterations $N_{\tau}$ ranges from $100$ to $3200$ . The comparison is displayed in Table 1, from which we observe the reduction of numerical error as the time step decreases for all considered settings. We note that all numerical values are one order of magnitude smaller to those reported in [43, Fig. 4 (SAV/BDF1, R-SAV/BDF1)] for a standard Cahn–Hilliard equation discretized with a first order time discretization. While RSAV having smaller error values than SAV, unlike in [43], the magnitude of the errors are almost indistinguishable among all cases, and the ratio of RSAV errors to SAV errors approaches 1 as the time step decreases. Thus, in this test case we do not see any clear advantage of RSAV over SAV in terms of numerical accuracy. Furthermore, the idealized choice $\zeta_{n}^{\tau}=0$ may not always achieve the smallest error among RSAV schemes for certain time steps.

We repeat our error comparison with two more examples, where the source functions are chosen specifically so that with $m(\varphi)=1$ , $g=0$ , $F^{\prime}(\varphi)=\varphi^{3}-\varphi$ and $\varepsilon=1$ , the analytical solutions are

\varphi^{*}(x,t)=\exp(\cos(\pi x))\cos(t)\quad\text{ and }\quad\exp(\cos(t))(\sin(\pi x))^{2},

respectively. In Figure 1(b) and (c), we plot the corresponding discrete Ginzburg–Landau energy for the above two analytical solutions with $h=\tau=0.01$ . Then, we set $h=0.001$ in the error comparison between RSAV and SAV schemes for these two examples, which can be found in Tables 2 and 3, respectively. We again note that the modified Cahn–Hilliard equation does not exhibit a dissipative structure with the Ginzburg–Landau functional, RSAV have smaller (but similar in magnitude) errors than SAV, the ratio of RSAV errors to SAV errors approaches 1 as the time step decreases, and the idealized choice $\zeta_{n}^{\tau}=0$ may not achieve the smallest errors among RSAV schemes for certain time steps.

Our findings indicate that in contrast to dissipative systems of the form (1.5), for Cahn–Hilliard equations with mass sources such as (1.1), RSAV schemes perform only marginally better than the standard SAV scheme in terms of numerical accuracy.

$(\tau,\zeta)$	$1$	$0.75$	$0.5$	$0.25$	$0$
0.1	0.1235112595	0.1223488639	0.1221746557	0.1221087366	0.1220743681
0.05	0.0657475964	0.0654207105	0.0653961915	0.0653875327	0.0653831138
0.025	0.0370353831	0.0369482551	0.0369450147	0.0369439043	0.0369433434
0.0125	0.0227271731	0.0227046023	0.0227041869	0.0227040466	0.0227039746
0.00625	0.0155923393	0.0155909899	0.0155845443	0.0156222078	0.0155860488
0.003125	0.0120314862	0.0120300314	0.0120300093	0.0120300181	0.0120300069
0.0015625	0.0102552922	0.0102549247	0.0102549247	0.0102549230	0.0102549241

Table 1: Comparison of

L^{2}(0,5;L^{2}(\Omega))

-error for various choices of constant

\zeta

with

\eta=0.95

and

M=1

. The optimal

\zeta_{n}^{\tau}

computed via (2.4) is

\zeta_{n}^{\tau}=0

in all scenarios. The smallest value in each row is highlighted in bold. The exact solution is

\varphi^{*}(x,t)=\cos(\pi x)(1+t)

$(\tau,\zeta)$	$1$	$0.75$	$0.5$	$0.25$	$0$
0.1	0.1580897113	0.1573562578	0.1573021781	0.1572764459	0.1572614308
0.05	0.0785704830	0.0784190660	0.0784103517	0.0784068449	0.0784050015
0.025	0.0391241631	0.0390891339	0.0390878945	0.0390874601	0.0390872404
0.0125	0.0195471355	0.0195389811	0.0195358830	0.0195387801	0.0195387550
0.00625	0.0100802056	0.0100784782	0.0100784619	0.0100784607	0.0100784602
0.003125	0.0057771876	0.0057769012	0.0057769104	0.0057769112	0.0057769075
0.0015625	0.0042246587	0.0042246609	0.0042246027	0.0042246434	0.0042246623

Table 2: Comparison of

L^{2}(0,5;L^{2}(\Omega))

-error for various choices of constant

\zeta

with

\eta=0.95

and

M=1

. The optimal

\zeta_{n}^{\tau}

computed via (2.4) is

\zeta_{n}^{\tau}=0

in all scenarios. The smallest value in each row is highlighted in bold. The exact solution is

\varphi^{*}(x,t)=\exp(\cos(\pi x))\cos(t)

$(\tau,\zeta)$	$1$	$0.75$	$0.5$	$0.25$	$0$
0.1	0.0743044854	0.0742724251	0.0742580052	0.0742503247	0.0742457608
0.05	0.0362455345	0.0362330117	0.0362305278	0.0362294871	0.0362289303
0.025	0.0175288174	0.0175263383	0.0175260131	0.0175258949	0.0175258344
0.0125	0.0085019366	0.0085014574	0.0085014301	0.0085014207	0.0085014149
0.00625	0.0046007411	0.0046007312	0.0046007299	0.0046007307	0.0046007395
0.003125	0.0034729641	0.0034730059	0.0034729916	0.0034729941	0.0034729946
0.0015625	0.0034123532	0.0034123725	0.0034123619	0.0034123592	0.0034123702

Table 3: Comparison of

L^{2}(0,5;L^{2}(\Omega))

-error for various choices of constant

\zeta

with

\eta=0.95

and

M=1

. The optimal

\zeta_{n}^{\tau}

computed via (2.4) is

\zeta_{n}^{\tau}=0

in all scenarios. The smallest value in each row is highlighted in bold. The exact solution is

\varphi^{*}(x,t)=\exp(\cos(t))(\sin(\pi x))^{2}

5 Applications

In this section we apply the SAV framework to various Cahn–Hilliard models that do not exhibit an obvious Lyapunov or dissipative structure due to the presence of source terms. In preliminary investigations we did not observe significant visual differences between the numerical solutions obtained from the standard SAV approach and the RSAV approach. Hence, to simplify the presentation, we propose time discretizations based on the standard SAV approach and demonstrate their stability. The corresponding RSAV-based time discretizations and their stability can be obtained with a straightforward modification. In turn, the convergence as $\tau\to 0$ to the corresponding weak solutions can be inferred similarly as in the proof of Theorem 3.1. We then display some qualitative simulations on a square domain $\Omega=[0,1]^{2}$ discretized with a standard Friedrichs–Keller triangulation and spatial step size $h=0.01$ . Unless specified, in all simulations below we take the quartic potential $F(s)=\frac{1}{4}(s^{2}-1)^{2}$ , and set constant $C_{0}=1$ in (1.3).

5.1 Microphase separation of diblock copolymer

A phase field model proposed to describe the dynamics of microphase separation of diblock copolymers is the following Cahn–Hilliard–Oono equation [19, 32]

\displaystyle\partial_{t}\varphi+\eta(\varphi-c)=\Delta\mu,\quad\mu=-\varepsilon\Delta\varphi+\varepsilon^{-1}F^{\prime}(\varphi),

(5.1)

furnished with homogeneous Neumann conditions. In the above $\eta>0$ is a constant related to the chain length of the copolymer and $c\in\mathbb{R}$ is a prescribed relative mass average. If $c=\langle\varphi_{0}\rangle_{\Omega}=\frac{1}{|\Omega|}\int_{\Omega}\varphi_{0}\,\mathrm{dx}$ equals the mean value of the initial condition, then the system is also known as the Ohta–Kawasaki equation [31], which is the conserved gradient flow of the Ohta–Kawaski functional

E(\varphi)=\int_{\Omega}\frac{\varepsilon}{2}|\nabla\varphi|^{2}+\frac{1}{\varepsilon}F(\varphi)\,\mathrm{dx}+\int_{\Omega\times\Omega}\eta G(x-y)(\varphi(x)-\langle\varphi\rangle_{\Omega})(\varphi(y)-\langle\varphi\rangle_{\Omega})\,\mathrm{dx}\,\mathrm{dy},

where $G$ is the Green function associated to the Neumann–Laplacian. A formal integration of the first equation in (5.1) yields a differential equation for the mean value $\langle\varphi\rangle_{\Omega}$ :

\displaystyle\frac{d\langle\varphi\rangle_{\Omega}}{dt}+\eta(\langle\varphi\rangle_{\Omega}-c)=0\quad\implies\quad\langle\varphi\rangle_{\Omega}(t)=c+e^{-\eta t}(\langle\varphi_{0}\rangle_{\Omega}-c)\;\forall t>0,

(5.2)

where we note the conservation of mass $\langle\varphi\rangle_{\Omega}(t)=\langle\varphi_{0}\rangle_{\Omega}$ occurs only when $c=\langle\varphi_{0}\rangle_{\Omega}$ .

A SAV-based time discretization of (5.1) reads as


$\displaystyle\varphi^{n}-\varphi^{n-1}$	$\displaystyle=\tau\,\Delta\mu^{n}+\tau\eta(c-\varphi^{n-1})$	(5.3a)
$\displaystyle\mu^{n}$	$\displaystyle=-\varepsilon\Delta\varphi^{n}+\frac{q^{n}}{\varepsilon Q(\varphi^{n-1})}F^{\prime}(\varphi^{n-1}),$	(5.3b)
$\displaystyle q^{n}-q^{n-1}$	$\displaystyle=\frac{1}{2\varepsilon Q(\varphi^{n-1})}(F^{\prime}(\varphi^{n-1}),\varphi^{n}-\varphi^{n-1}),$	(5.3c)

furnished with homogeneous Neumann conditions. An interesting feature of this SAV discretization is that the source term need not be a bounded function, as we have the following discrete analogue of (5.2):

\langle\varphi^{n}\rangle_{\Omega}=(1-\tau\eta)\langle\varphi^{n-1}\rangle_{\Omega}+\tau\eta c\quad\implies\quad\langle\varphi^{n}\rangle_{\Omega}=\langle\varphi_{0}\rangle_{\Omega}(1-\tau\eta)^{n}+c(1-(1-\tau\eta)^{n}).

For $\tau<\eta^{-1}$ we have the a priori estimate $|\langle\varphi^{n}\rangle_{\Omega}|\leq C$ for all $n\in\mathbb{N}$ . Then, testing (5.3a) with $\mu^{n}$ and $\varphi^{n}$ , (5.3b) with $\varphi^{n}-\varphi^{n-1}$ , (5.3c) with $2q^{n}$ , and upon summing we obtain for $G^{k}$ defined as in (2.8):

	$\displaystyle G^{n}-G^{n-1}+\frac{1}{2}\big{(}\\|\varphi^{n}-\varphi^{n-1}\\|^{2}+\varepsilon\\|\nabla(\varphi^{n}-\varphi^{n-1})\\|^{2}\big{)}+\|q^{n}-q^{n-1}\|^{2}+\tau\\|\nabla\mu^{n}\\|^{2}$
	$\displaystyle\quad=\tau\eta\Big{(}(c-\langle\varphi^{n-1}\rangle_{\Omega},\mu^{n}+\varphi^{n})+(\langle\varphi^{n-1}\rangle_{\Omega}-\varphi^{n-1},\mu^{n}+\varphi^{n})\Big{)}$
	$\displaystyle\quad\leq C\tau\big{(}\\|\mu^{n}\\|_{L^{1}}+\\|\varphi^{n}\\|\big{)}+C\tau\\|\nabla\varphi^{n-1}\\|(\\|\nabla\mu^{n}\\|+\\|\varphi^{n}\\|),$

on account of the Poincaré inequality and

(\langle\varphi^{n-1}\rangle_{\Omega}-\varphi^{n-1},\mu^{n})=(\langle\varphi^{n-1}\rangle_{\Omega}-\varphi^{n-1},\mu^{n}-\langle\mu^{n}\rangle_{\Omega})\leq C\|\nabla\varphi^{n-1}\|\|\nabla\mu^{n}\|.

Upon using (2.10) and (2.14) with $m_{0}=1$ , $g=0$ , $\zeta=0$ , $r^{n}=q^{n}$ , we infer

		$\displaystyle G^{n}-G^{n-1}+\frac{\min(1,\varepsilon)}{2}\\|\varphi^{n}-\varphi^{n-1}\\|_{H^{1}}^{2}+\frac{3}{4}\|q^{n}-q^{n-1}\|^{2}+\frac{\tau}{2}\\|\nabla\mu^{n}\\|^{2}$		(5.4)
		$\displaystyle\quad\leq C\tau(1+G^{n-1}+G^{n}),$		(5.4)

which yields stability of the SAV scheme for $\tau$ sufficiently small. Notice that above calculation relies on the uniform control of the mean value of $\varphi^{n}$ , which is possible here due to the specific linear structure of the source term $f(\varphi)=\eta(c-\varphi)$ .

We consider the scheme (5.3) with parameter values $\eta=0.001$ and $\varepsilon=0.01$ , time step $\tau=0.01$ , and $c=\langle\varphi^{0}\rangle_{\Omega}$ is the mean value of the initial condition $\varphi^{0}$ . Taking $\varphi^{0}$ as the perturbation of $-0.5$ (Figure 3 top row) and of $-0.1$ (Figure 3 bottom row) with uniformly distributed random numbers in $(0,0.2)$ , we display the discrete solution $\varphi^{n}$ at $n=0$ (initial condition), $500$ , $10000$ and $50000$ in Figure 3. The behavior of the discrete solution for both settings is similar to the simulations in [24, Sec. 3].

5.2 Image segmentation

In [41], the authors proposed an image segmentation model with the following Cahn–Hilliard equation

\displaystyle\partial_{t}\varphi=\Delta\mu-\frac{\eta\Big{(}\lambda_{1}(I(x)-c_{1})^{2}-\lambda_{2}(I(x)-c_{2})^{2}\Big{)}}{\pi(\eta^{2}+(\varphi-0.5)^{2})},\quad\mu=-\varepsilon\Delta\varphi+\varepsilon^{-1}F^{\prime}(\varphi),

(5.5)

furnished with homogeneous Neumann conditions, $F(s)=2s^{2}(s-1)^{2}$ , positive constants $\lambda_{1},\lambda_{2}$ , grayscale image function $I$ rescaled to the range $[0,1]$ , and $c_{1}$ and $c_{2}$ represent averaged pixel intensities inside and outside the segmented regions. The source term can be interpreted as a $L^{2}$ -gradient of the energy functional

M(u)=\lambda_{1}\int_{\Omega}(I(x)-c_{1})^{2}H_{\eta}(u-\tfrac{1}{2})\,\mathrm{dx}+\lambda_{2}\int_{\Omega}(I(x)-c_{2})^{2}[1-H_{\eta}(u-\tfrac{1}{2})]\,\mathrm{dx}

where $H_{\eta}(z)=\frac{1}{2}(1+\frac{2}{\pi}\arctan(z/\eta))$ is a $C^{\infty}$ -regularization of the Heaviside function, and $\delta_{\eta}(z)=\eta/(\pi(\eta^{2}+z^{2}))=H_{\eta}^{\prime}(z)$ approximates the Dirac measure. The average intensities $c_{1}$ and $c_{2}$ are defined as

\displaystyle c_{1}(\varphi)=\frac{(I,H_{\eta}(\varphi-\frac{1}{2}))}{(1,H_{\eta}(\varphi-\frac{1}{2}))},\quad c_{2}(\varphi)=\frac{(I,1-H_{\eta}(\varphi-\frac{1}{2}))}{(1,1-H_{\eta}(\varphi-\frac{1}{2}))}.

(5.6)

Due to the bounded source term, a SAV-based discretization of (5.5) reads as

$\displaystyle\varphi^{n}-\varphi^{n-1}$	$\displaystyle=\tau\Delta\mu^{n}-\tau\frac{\eta\Big{(}\lambda_{1}(I(x)-c_{1}^{n-1})^{2}-\lambda_{2}(I(x)-c_{2}^{n-1})^{2}\Big{)}}{\pi(\eta^{2}+(\varphi^{n-1}-0.5)^{2})},$	(5.7)
$\displaystyle\mu^{n}$	$\displaystyle=-\varepsilon\Delta\varphi^{n}+\frac{q^{n}}{\varepsilon Q(\varphi^{n-1})}F^{\prime}(\varphi^{n-1}),$
$\displaystyle q^{n}-q^{n-1}$	$\displaystyle=\frac{1}{2\varepsilon Q(\varphi^{n-1})}(F^{\prime}(\varphi^{n-1}),\varphi^{n}-\varphi^{n-1}),$

furnished with homogeneous Neumann conditions. Then, the average intensities are updated to $c_{1}^{n}:=c_{1}(\varphi^{n})$ and $c_{2}^{n}:=c_{2}(\varphi^{n})$ via (5.6). It is straightforward to infer the estimate (5.4) for $G^{k}$ , leading to stability of the SAV scheme for sufficiently small $\tau$ .

We apply (5.5) to segment images of a cow (Figure 4 top row) and of a blood vessel (Figure 4 bottom row), where for initialization we take $c_{1}=1$ and $c_{2}=0$ . Similar to [41], in these experiments we choose parameters $\eta=0.1$ , $\lambda_{1}=0.65$ , $\lambda_{2}=1$ and time step $\tau=0.001$ . A two-stage procedure is used where we first solve (5.7) with a large value of $\varepsilon=80$ , and once a quasi-steady state has been achieved (approximately at $n=5000$ ), we reduce the value of $\varepsilon$ to $0.01$ and resume the iterative process until a new steady state is reached. In Figure 4 we display the $1/2$ -level set of discrete solution $\varphi^{n}$ in red overlayed with the original image at iterations $n=1000,3000$ and $10000$ , as well as the discrete solution itself at $n=10000$ . We note that the segmentation task is well-performed in both experiments, and for a comparison between the Cahn–Hilliard approach (5.5) with the classical Chan–Vese approach we refer the reader to [41].

5.3 Image inpainting

In [3] the authors proposed the following Cahn–Hilliard model to restore missing or damaged details in a measurable subdomain $D\subset\Omega$ for a binary image function $I:\Omega\to\{-1,1\}$ :

\displaystyle\partial_{t}\varphi=\Delta\mu+\lambda(x)(I(x)-\varphi),\quad\mu=-\varepsilon\Delta\varphi+\varepsilon^{-1}F^{\prime}(\varphi),\quad\lambda(x)=\lambda_{0}\chi_{\Omega\setminus D}(x),

(5.8)

furnished with homogeneous Neumann conditions. In the above $\chi_{A}(x)$ denotes the characteristic function of the set $A$ , and the fidelity term $\lambda(x)(I(x)-\varphi)$ with large constant $\lambda_{0}$ demands the solution $\varphi$ to be close to the data $I(x)$ in the undamaged region $\Omega\setminus D$ . In contrast to (5.1) there is no analogue to (5.2) for the mean value $\langle\varphi\rangle_{\Omega}$ due to the function $\lambda(x)=\lambda_{0}\chi_{\Omega\setminus D}(x)$ (see [7] for an inequality estimate) and thus our stability analysis does not apply to a standard SAV discretization of (5.8). On account of the purpose of the fidelity term, a reasonable modification of (5.8) is

\displaystyle\partial_{t}\varphi=\Delta\mu+\lambda(x)(I(x)-\mathcal{T}(\varphi)),\quad\mathcal{T}(\varphi)=\max(-1,\min(1,\varphi)).

(5.9)

For the modified equation (5.9), the source term is bounded and so a SAV-based discretization reads as

$\displaystyle\varphi^{n}-\varphi^{n-1}$	$\displaystyle=\tau\Delta\mu^{n}+\tau\lambda(x)(I(x)-\mathcal{T}(\varphi^{n-1})),$	(5.10)
$\displaystyle\mu^{n}$	$\displaystyle=-\varepsilon\Delta\varphi^{n}+\frac{q^{n}}{\varepsilon Q(\varphi^{n-1})}F^{\prime}(\varphi^{n-1}),$
$\displaystyle q^{n}-q^{n-1}$	$\displaystyle=\frac{1}{2\varepsilon Q(\varphi^{n-1})}(F^{\prime}(\varphi^{n-1}),\varphi^{n}-\varphi^{n-1}),$

furnished with homogeneous Neumann conditions. It is straightforward to infer the estimate (5.4) for $G^{k}$ , leading to stability of the SAV scheme for sufficiently small $\tau$ .

We apply (5.10) to perform inpainting of a double stripe shown in Figure 5(a). Similar to image segmentation, we adopt a two-stage procedure where we first solve (5.10) with parameters $\lambda_{0}=10$ , $\varepsilon=100$ , $\tau=0.1$ , and after $n=3000$ iterations we then switch a new set of parameters $\lambda_{0}=0.1$ , $\varepsilon=5$ , $\tau=1$ . Figure 5(d) displays the discrete solution $\varphi^{n}$ at $n=4000$ and is an acceptable reconstruction result, cf. [24].

5.4 Tumor growth dynamics

In [9, 17, 40] the following model is used to describe tumor growth dynamics:


$\displaystyle\partial_{t}\varphi$	$\displaystyle=\mathrm{div}\,(m(\varphi)\nabla\mu)+P(\varphi,\mu,\sigma),$	(5.11a)
$\displaystyle\mu$	$\displaystyle=-\varepsilon\Delta\varphi+\varepsilon^{-1}F^{\prime}(\varphi)-\chi\sigma,$	(5.11b)
$\displaystyle\partial_{t}\sigma$	$\displaystyle=\mathrm{div}\,(n(\varphi)\nabla(\chi_{\sigma}\sigma-\eta\varphi))-S(\varphi,\mu,\sigma).$	(5.11c)

In the above, $\varphi\in[-1,1]$ represents the difference in volume fraction between tumor and healthy tissues; $\mu$ is the associated chemical potential; $\sigma$ denotes the concentration of a nutrient for the tumor; $m(\varphi)$ and $n(\varphi)$ denote positive mobilities for $\varphi$ and $\sigma$ , bounded above and below by constants $m_{1},m_{0}$ and $n_{1},n_{0}$ , respectively; parameters $\chi_{\sigma}>0$ , $\chi>0$ and $\eta\geq 0$ model the nutrient diffusivity, chemotaxis and active transport mechanisms respectively; $P$ accounts for tumor proliferation and apoptosis; $S$ models nutrient consumption. When $\chi=\eta$ , and furnishing with homogeneous Neumann conditions, the model admits the following energy identity

	$\displaystyle\frac{d}{dt}\int_{\Omega}\frac{\varepsilon}{2}\|\nabla\varphi\|^{2}+\frac{1}{\varepsilon}F(\varphi)+\frac{1}{2}\|\sigma\|^{2}-\chi\varphi\sigma\,\mathrm{dx}+\int_{\Omega}m(\varphi)\|\nabla\mu\|^{2}+n(\varphi)\|\nabla(\chi_{\sigma}\sigma-\chi\varphi)\|^{2}\,\mathrm{dx}$
	$\displaystyle\quad=\int_{\Omega}P(\varphi,\mu,\sigma)\mu-S(\varphi,\mu,\sigma)(\sigma-\eta\varphi)\,\mathrm{dx}.$

The forms of the source terms $P$ and $S$ determine whether $E(\varphi,\sigma)=\int_{\Omega}\frac{\varepsilon}{2}|\nabla\varphi|^{2}+\frac{1}{\varepsilon}F(\varphi)+\frac{1}{2}|\sigma|^{2}-\chi\varphi\sigma\,\mathrm{dx}$ is a Lyapunov functional. One example proposed in [22] is to consider $P(\varphi,\mu,\sigma)=S(\varphi,\mu,\sigma)=Q(\varphi)(\sigma-\chi\varphi-\mu)=Q(\varphi)(\frac{\delta E}{\delta\sigma}-\frac{\delta E}{\delta\varphi})$ as a difference of chemical potentials weighted by a bounded, nonnegative function $Q$ . This choice would yield $\frac{d}{dt}E(\varphi,\sigma)\leq 0$ , and the standard SAV schemes can be applied, see [36]. Another, more phenomenological in nature, is the following used in e.g. [17, 40]:

\displaystyle P(\varphi,\sigma)=h(\varphi)(\mathcal{P}k(\sigma)-\mathcal{A}),\quad S(\varphi,\sigma)=\mathcal{C}h(\varphi)\sigma,

(5.12)

with constant proliferation rate $\mathcal{P}$ , apoptosis rate $\mathcal{A}$ , nutrient consumption rate $\mathcal{C}$ and bounded, nonnegative functions $h$ and $k$ . One example is $h(s)=\max(-1,\min(1,\tfrac{1}{2}(1+s)))$ and $k(s)=\max(0,\min(\sigma_{\infty},s))$ with fixed external supply concentration $\sigma_{\infty}$ .

We consider a reformulation of the model with a new variable $\hat{\mu}:=\mu+\chi\sigma$ , where a SAV-based time discretization of (5.11) with source terms of the form (5.12) is


$\displaystyle\sigma^{n}-\sigma^{n-1}$	$\displaystyle=\tau\,\mathrm{div}\,(n(\varphi^{n-1})\nabla(\chi_{\sigma}\sigma^{n}-\eta\varphi^{n-1}))-\tau S(\varphi^{n-1},\sigma^{n}),$	(5.13a)
$\displaystyle\varphi^{n}-\varphi^{n-1}$	$\displaystyle=\tau\,\mathrm{div}\,(m(\varphi^{n-1})\nabla(\hat{\mu}^{n}-\chi\sigma^{n}))+\tau P(\varphi^{n-1},\sigma^{n})$	(5.13b)
$\displaystyle\hat{\mu}^{n}$	$\displaystyle=-\varepsilon\Delta\varphi^{n}+\frac{q^{n}}{\varepsilon Q(\varphi^{n-1})}F^{\prime}(\varphi^{n-1}),$	(5.13c)
$\displaystyle q^{n}-q^{n-1}$	$\displaystyle=\frac{1}{2\varepsilon Q(\varphi^{n-1})}(F^{\prime}(\varphi^{n-1}),\varphi^{n}-\varphi^{n-1}),$	(5.13d)

furnished with homogeneous Neumann conditions. Note that we allow $\eta$ and $\chi$ to be different constants. Given $(\varphi^{n-1},q^{n-1},\sigma^{n-1})$ , we first solve for $\sigma^{n}$ with (5.13a), and then solve for $(\varphi^{n},\hat{\mu}^{n},q^{n})$ with (5.13b)-(5.13d). As $P(\varphi,\sigma)$ is bounded and $|S(\varphi,\sigma)|\leq C|\sigma|$ , by testing (5.13b) with $\hat{\mu}^{n}$ and $\varphi^{n}$ , (5.13c) with $\varphi^{n}-\varphi^{n-1}$ , (5.13d) with $2q^{n}$ and (5.13a) with $K\sigma^{n}$ for some constant $K>0$ yet to be specified, upon summing we obtain for

G^{k}:=\frac{1}{2}\|\varphi^{k}\|^{2}+\frac{\varepsilon}{2}\|\nabla\varphi^{k}\|^{2}+|q^{k}|^{2}+\frac{K}{2}\|\sigma^{k}\|^{2},

the estimate

	$\displaystyle G^{n}-G^{n-1}+\frac{1}{2}\big{(}\\|\varphi^{n}-\varphi^{n-1}\\|^{2}+\varepsilon\\|\nabla(\varphi^{n}-\varphi^{n-1})\\|^{2}+K\\|\sigma^{n}-\sigma^{n-1}\\|^{2}\big{)}$
	$\displaystyle\qquad+\|q^{n}-q^{n-1}\|^{2}\ +\tau\int_{\Omega}m(\varphi^{n-1})\|\nabla\hat{\mu}^{n}\|^{2}+K\chi_{\sigma}n(\varphi^{n-1})\|\nabla\sigma^{n}\|^{2}\,\mathrm{dx}$
	$\displaystyle\quad=\tau(P(\varphi^{n-1},\sigma^{n-1}),\hat{\mu}^{n}+\varphi^{n})-\tau(m(\varphi^{n-1})\nabla\hat{\mu}^{n},\nabla(\varphi^{n}-\chi\sigma^{n}))$
	$\displaystyle\qquad+\tau(\nabla\sigma^{n},K\eta n(\varphi^{n-1})\nabla\varphi^{n-1}+\chi m(\varphi^{n-1})\nabla\varphi^{n})-K\tau(S(\varphi^{n-1},\sigma^{n}),\sigma^{n})$
	$\displaystyle\quad\leq C\tau\\|\hat{\mu}^{n}\\|_{L^{1}}+\frac{\tau m_{0}}{4}\\|\nabla\hat{\mu}^{n}\\|^{2}+\tau\Big{(}\frac{K\chi_{\sigma}n_{0}}{4}+c\Big{)}\\|\nabla\sigma^{n}\\|^{2}$
	$\displaystyle\qquad+C\tau\big{(}\\|\varphi^{n}\\|_{H^{1}}^{2}+\\|\nabla\varphi^{n-1}\\|^{2}+\\|\sigma^{n}\\|^{2}\big{)},$

where $c>0$ is independent of $K$ . Choosing $K$ sufficiently large so that $c<\frac{K\chi_{\sigma}n_{0}}{4}$ , and using (2.10) and (2.14) for $\|\hat{\mu}^{n}\|_{L^{1}}$ we arrive at

	$\displaystyle G^{n}-G^{n-1}+\frac{\min(1,\varepsilon)}{2}\\|\varphi^{n}-\varphi^{n-1}\\|_{H^{1}}^{2}+\frac{K}{2}\\|\sigma^{n}-\sigma^{n-1}\\|^{2}$
	$\displaystyle\qquad+\|q^{n}-q^{n-1}\|^{2}\ +\frac{\tau m_{0}}{2}\\|\nabla\hat{\mu}^{n}\\|^{2}+\frac{\tau K\chi_{\sigma}n_{0}}{2}\\|\nabla\sigma^{n}\\|^{2}$
	$\displaystyle\quad\leq C\tau\big{(}1+G^{n}+G^{n-1}\big{)},$

which ensures the stability of the SAV scheme for sufficiently small $\tau$ .

Remark 5.1.

If $\eta=0$ and $\sigma^{0}$ is bounded and non-negative, then by applying a comparison principle to (5.13a) we can deduce that $\sigma^{n}$ is also bounded and non-negative for all $n$ . Thus, we can replace $k(\sigma)$ in (5.12) with simply $\sigma$ .

We consider a similar setting to [18] where we set $\chi_{\sigma}=25$ , $\eta=\chi=5$ , constant mobilities $n(\varphi)=m(\varphi)=1$ , $\varepsilon=0.01$ , $\tau=0.001$ , with initial conditions $\sigma^{0}=1$ and

\varphi^{0}(\bm{x})=-\tanh\Big{(}\frac{|\bm{x}|-(0.05+0.02\cos(2\theta))}{\sqrt{2}\times 0.01}\Big{)},\quad\bm{x}=|\bm{x}|(\cos\theta,\sin\theta)^{\top}.

Such a choice for the initial condition $\varphi^{0}$ means the tumor is initially bounded by a slightly perturbed circle. The source functions are chosen as $P(\varphi,\sigma)=S(\varphi,\sigma)=\frac{1}{2}k(\sigma)(1+\varphi)$ with $\sigma_{\infty}=1$ . In Figure 6 we display the discrete solution $\varphi^{n}$ at $n=0$ (initial condition), $n=8000$ , $n=13000$ and $n=18000$ , where the tumor undergoes morphological instabilities and develop fingers towards regions of high nutrient concentration. Not shown here are plots of the discrete solution $\sigma^{n}$ , which are visually similar to Figure 6, namely the nutrient concentration $\sigma^{n}$ is lower in the tumor region $\{\varphi^{n}=1\}$ and higher in the tissue region $\{\varphi^{n}=-1\}$ . We note that these seem to resemble previous simulations found in [17, 18, 40].

6 Conclusion

In this work we investigate the stability of a time discretization based on the relaxed scalar auxiliary variable (RSAV) approach of [15] for Cahn–Hilliard systems with bounded mass source. In general these systems may not exhibit dissipative structures, and so the stability of such SAV-based schemes are not immediate. Our proofs rely on a key estimate of a product term between the scalar auxiliary variable and the cubic nonlinearity. For the convergence of time discrete solutions our analysis covers two choices of the interpolating parameter $\zeta_{n}^{\tau}$ in (2.2). Numerical simulations are supportive of our results and we are able to replicate the expected solution behavior for Cahn–Hilliard systems in tumor growth, image segmentation and image inpainting.

Acknowledgments

The authors gratefully acknowledge the support by the Research Grants Council of the Hong Kong Special Administrative Region, China [Project No.: HKBU 14302319] and Hong Kong Baptist University [Project No.: RC-OFSGT2/20-21/SCI/006].

References

[1] A.C. Aristotelous, O.A. Karakashian and S.M. Wise, Adaptive, second-order in time, primitive-variable discontinuous Galerkin schemes for a Cahn–Hilliard equation with a mass source, IMA J. Numer. Anal. 35 (2015), 1167–1198.
[2] A. Bertozzi, S. Esdoḡlu and A. Gillette, Analysis of a two-scale Cahn–Hilliard model for binary image inpainting, Multiscale Model. Simul. 6 (2007), 913–936.
[3] A. Bertozzi, S. Esdoḡlu and A. Gillette, Inpainting of binary images using the Cahn–Hilliard equation, IEEE Trans. Image Process. 16 (2007), 285–291.
[4] H. Chen, J. Mao and J. Shen, Optimal error estimates for the scalar auxiliary variable finite-element schemes for gradient flows, Numer. Math. 145 (2020), 167–196
[5] Q. Cheng, C. Liu and J. Shen, A new Lagrange multiplier approach for gradient flows, Comput. Method Appl. Math. Engrg. 367 (2020), 113070.
[6] Q. Cheng, X. Yang and J. Shen, Efficient and accurate numerical schemes for a hydro-dynamically coupled phase field diblock copolymer model, J. Comput. Phys. 341 (2017), 44–60.
[7] L. Cherfils, H. Fakih and A. Miranville, Finite-dimensional attractors for the Bertozzi-Esedoglu-Gillette-Cahn-Hilliard equation in image inpainting, Inv. Probl. Imag. 9 (2015, 105–125.
[8] P. Clément, Approximation by finite element functions using local regularization, ESAIM: M2NA 9 (1975), 77–84.
[9] V. Cristini and J. Lowengrub, Multiscale Modeling of Cancer: An Integrated Experimental and Mathematical Modeling Approach, Cambridge University Press, 2010.
[10] C.M. Elliott and A.M. Stuart, The global dynamics of discrete semilinear parabolic equations, SIAM J. Numer. Anal. 30 (1993), 1622–1663.
[11] D.J. Eyre, Unconditionally gradient stable time marching the Cahn-Hilliard equation, Mater. Res. Soc. Symp. Proc. 529 (1998), 39–46.
[12] H. Fakih, R. Mghames and N. Nasreddine, On the Cahn–Hilliard equation with mass source for biological applications, Commun. Pure Appl. Anal. 20 (2021), 495–510.
[13] F. Huang and J. Shen, A new class of implicit-explicit BDFk SAV schemes for general dissipative systems and their error analysis, Comput. Methods Appl. Mech. Engrg. 392 (2022), 114718.
[14] F. Huang, J. Shen and Z. Yang, A highly efficient and accurate new scalar auxiliary variable approach for gradient flows, SIAM J. Sci. Comput. 42 (2020), A2514–A2536.
[15] M. Jiang, Z. Zhang and J. Zhao, Improving the accuracy and consistency of the scalar auxiliary variable (sav) method with relaxation, J. Comput. Phys. (2022), 110954.
[16] H. Garcke and K.F. Lam, Well-posedness of a Cahn–Hilliard system modelling tumour growth with chemotaxis and active transport, Eur. J. Appl. Math. 28 (2017), 284–316.
[17] H. Garcke, K.F. Lam, E. Sitka and V. Styles, A Cahn–Hilliard–Darcy model for tumour growth with chemotaxis and active transport, Math. Models Methods Appl. Sci. 26 (2016), 1095–1148.
[18] H. Garcke, and D. Trautwein, Numerical analysis for a Cahn–Hilliard system modelling tumour growth with chemotaxis and active transport, J. Numer. Math. 30 (2022), 295–324.
[19] A. Giorgini, M. Grasselli and A. Miranville, The Cahn–Hilliard–Oono equation with singular potential, Math. Models Methods Appl. Sci. 27 (2017), 2485–2510.
[20] P. Grisvard, Elliptic Problems in Nonsmooth Domains, SIAM, Philadelphia, 2011.
[21] F. Guillén-González and G. Tierra, On linear schemes for a Cahn–Hilliard diffuse interface model, J. Comput. Phys. 234 (2013), 140–171.
[22] A. Hawkins-Daarud, K.G. van der Zee and J.T. Oden, Numerical simulation of a thermodynamically consistent four-species tumor growth model, Int. J. Numer. Method Biomed. Eng. 28 (2012), 3–24.
[23] D. Hou, M. Azaiex and C. Xu, A variant of scalar auxiliary variable approaches for gradient flows, J. Comput. Phys. 395 (2019), 307–332.
[24] J. Kim, S. Lee, Y. Choi, S.-M. Lee and D. Jeong, Basic Principles and Practical Applications of the Cahn–Hilliard Equation, Math. Probl. Eng. 2016 (2016), Article ID 9532608, 11 pages.
[25] P. Knopf, K.F. Lam, C. Liu and S. Metzger, Phase-field dynamics with transfer of materials: The Cahn–Hilliard equation with reaction rate dependent dynamical boundary conditions, ESAIM: M2AN 55 (2021), 229–282.
[26] X. Li and J. Shen, On a SAV-MAC scheme for the Cahn–Hilliard–Navier–Stokes phase field model and its error analysis for the corresponding Cahn–Hilliard–Stokes case, Math. Models Methods Appl. Sci. 30 (2020), 2263–2297.
[27] Z. Liu and X. Li, The exponential scalar auxiliary variable (E-SAV) approach for phase field models and its explicit computing, SIAM J. Sci. Comput. 42 (2020), pp. B630–B655.
[28] S. Metzger, A convergent SAV scheme for Cahn–Hilliard equations with dynamic boundary conditions, IMA J. Numer. Anal., (2023), doi:10.1093/imanum/drac078.
[29] L. Li, A. Miranville and R. Guillevin, Cahn–Hilliard models for glial cells, Appl. Math. Optim. 84 (2021), 1821–1842.
[30] A. Miranville, E. Rocca and G. Schimperna, On the long time behavior of a tumor growth model, J. Differential Equations 267 (2019), 2616–2642.
[31] T. Ohta and K. Kawasaki, Equilibrium morphology of block copolymer melts, Marcomolecules 19 (1986), 2621–2632.
[32] Y. Oono and S. Puri, Computationally efficient modeling of ordering of quenched phases, Phys. Rev. Lett. 58 (1987), 836–839.
[33] J. Shen and J. Xu, Convergence and error analysis of the scalar auxiliary variable (SAV) schemes to gradient flows, SIAM J. Numer. Anal. 56 (2018), 2895–2912.
[34] J. Shen, J. Xu and J. Yang, The scalar auxiliary variable (SAV) approach for gradient flows, J. Comput. Phys. 353 (2018), 407–416.
[35] J. Shen, J. Xu and J. Yang, A new class of efficient and robust energy stable schemes for gradient flows, SIAM Review 61 (2019), 474–506.
[36] X. Shen, L. Wu, J. Wen and J. Zhang, SAV Fourier-spectral method for diffuse-interface tumor-growth model, Comput. Math. with Appl. (2022) DOI: 10.1016/j.camwa.2022.09.031.
[37] J. Shen and X. Yang, Numerical approximations of Allen-Cahn and Cahn-Hilliard equations, Discrete Contin. Dyn. Syst. 28 (2010), 1669–1691.
[38] J. Shen and X. Yang, The IEQ and SAV approaches and their extensions for a class of highly nonlinear gradient flow systems, Contemp. Math. 754 (2020), 217–245.
[39] J. Simon, Compact Sets in the Space $L^{p}(0,T;B)$ , Ann. di Mat. Pura ed Appl. (IV) 146 (1987), 65–96.
[40] S.M. Wise, J.S. Lowengrub, H.B. Frieboes and V. Cristini, Three-dimensional multispecies nonlinear tumor growth - I: Model and numerical method, J. Theor. Biol. 253 (2008), 524–543.
[41] W. Yang, Z. Huang and W. Zhu, Image segmentation using the Cahn–Hilliard equation, J. Sci. Comput. 79 (2019), 1057–1077.
[42] Z. Yang and S. Dong, A roadmap for discretely energy-stable schemes for dissipative systems based on a generalized auxiliary variable with guaranteed positivity, J. Comput. Phys. 404 (2020), 109121.
[43] Y. Zhang and J. Shen, A generalized SAV approach with relaxation for dissipative systems, J. Comput. Phys. 464 (2022), 111311.
[44] N. Zheng and X. Li, Error analysis of the SAV Fourier-spectral method for the Cahn-Hilliard-Hele-Shaw system, Adv. Comput. Math. 47 (2021), 71.

	$\displaystyle\|(f^{n-1},\mu^{n})\|$	$\displaystyle\leq C\\|\mu^{n}\\|_{L^{1}}\leq C\\|\nabla\mu^{n}\\|_{L^{2}}+C\|(\mu^{n},1)\|$		(2.10)
		$\displaystyle\leq C+\frac{(1-\eta)m_{0}}{4}\\|\nabla\mu^{n}\\|^{2}+C\|(\mu^{n},1)\|.$		(2.10)

$\displaystyle C\frac{\|r^{n}\|\|(F^{\prime}(\varphi^{n-1}),1)\|}{\varepsilon Q(\varphi^{n-1})}$	$\displaystyle\leq C\frac{\|r^{n}\|((F(\varphi^{n-1}),1)^{\frac{3}{4}}+1)}{Q(\varphi^{n-1})}\leq C\|r^{n}\|\,(Q(\varphi^{n-1})^{\frac{1}{2}}+1)$	(2.13)
	$\displaystyle\leq\frac{1}{8}\|r^{n}\|^{2}+C\big{(}1+Q(\varphi^{n-1})\big{)}\leq\frac{1}{8}\|r^{n}\|^{2}+C\big{(}1+(F(\varphi^{n-1}),1)^{1/2}\big{)}$
	$\displaystyle\leq\frac{1}{8}\|r^{n}\|^{2}+C\big{(}1+\\|\varphi^{n-1}\\|_{L^{4}}^{2}\big{)}\leq\frac{1}{8}\|r^{n}\|^{2}+C\big{(}1+\\|\varphi^{n-1}\\|_{H^{1}}^{2}\big{)}.$

	$\displaystyle C\|(\mu^{n},1)\|$	$\displaystyle\leq C+C\\|g^{n-1}\\|^{2}+\frac{1}{8}\|r^{n}\|^{2}+C\\|\varphi^{n-1}\\|_{H^{1}}^{2}$		(2.14)
		$\displaystyle\leq C+C\\|g^{n-1}\\|^{2}+\frac{1}{4}\|r^{n}-q^{n-1}\|^{2}+C\|q^{n-1}\|^{2}+C\\|\varphi^{n-1}\\|_{H^{1}}^{2},$		(2.14)

	$\displaystyle\tau\|(f^{n-1},\mu^{n})\|$	$\displaystyle\leq C\tau+\frac{(1-\eta)\tau m_{0}}{4}\\|\nabla\mu^{n}\\|^{2}+\frac{\tau}{4}\|r^{n}-q^{n-1}\|^{2}$		(2.15)
		$\displaystyle\quad+C\tau\\|g^{n-1}\\|^{2}+C\tau\|q^{n-1}\|^{2}+C\tau\\|\varphi^{n-1}\\|_{H^{1}}^{2}.$		(2.15)

	$\displaystyle\frac{\varepsilon}{2}\\|\Delta\varphi^{n}\\|^{2}$	$\displaystyle\leq C\\|g^{n-1}\\|^{2}+C\frac{\|r^{n}\|^{2}}{Q(\varphi^{n-1})^{2}}\\|F^{\prime}(\varphi^{n-1})\\|^{2}$		(2.19)
		$\displaystyle\quad+C\\|\nabla\mu^{n}\\|\\|\nabla\varphi^{n}\\|\leq C+C\\|\nabla\mu^{n}\\|,$		(2.19)

Stability and convergence of relaxed scalar auxiliary variable schemes for Cahn–Hilliard systems with bounded mass source

Abstract

1 Introduction

2 RSAV scheme

2.1 First order time discretization

Remark 2.1.

2.2 Stability estimates

Remark 2.2 (Unique solvability).

Remark 2.3.

Remark 2.4.

Theorem 2.1 (Stability).

Proof.

Remark 2.5 (Sharpness of (A5\mathrm{A5})).

2.3 Discussion on the polynomial growth of FF

Remark 2.6.

3 Convergence of the time discrete solutions

Theorem 3.1 (Convergence).

Proof.

4 Numerical discussions

4.1 Fully discrete finite element approximation

4.2 Convergence test

5 Applications

5.1 Microphase separation of diblock copolymer

5.2 Image segmentation

5.3 Image inpainting

5.4 Tumor growth dynamics

Remark 5.1.

6 Conclusion

Acknowledgments

References

Remark 2.5 (Sharpness of ( $\mathrm{A5}$ )).

2.3 Discussion on the polynomial growth of $F$