High-order nonuniform time-stepping and MBP-preserving linear schemes for the time-fractional Allen-Cahn equation

TL;DR

This paper introduces high-order, nonuniform time-stepping linear schemes for the time-fractional Allen-Cahn equation that preserve energy stability and the maximum-bound principle, with improved MBP preservation for large time steps.

Contribution

The paper develops a new prediction strategy and auxiliary functional to design unconditionally energy-stable, MBP-preserving linear schemes with enhanced large-time-step MBP preservation.

Findings

The L1 scheme unconditionally preserves the discrete MBP.

The L2-1σ scheme requires a mild time-step restriction.

An improved L2-1σ scheme enhances MBP preservation for large time steps.

Abstract

In this paper, we present a class of nonuniform time-stepping, high-order linear stabilized schemes that can preserve both the discrete energy stability and maximum-bound principle (MBP) for the time-fractional Allen-Cahn equation. To this end, we develop a new prediction strategy to obtain a second-order and MBP-preserving predicted solution, which is then used to handle the nonlinear potential explicitly. Additionally, we introduce an essential nonnegative auxiliary functional that enables the design of an appropriate stabilization term to dominate the predicted nonlinear potential, and thus to preserve the discrete MBP. Combining the newly developed prediction strategy and auxiliary functional, we propose two unconditionally energy-stable linear stabilized schemes, L1 and L2-$1_\sigma$ schemes. We show that the L1 scheme unconditionally preserves the discrete MBP, whereas the L2-$1_\sigma$ scheme requires a mild time-step restriction. Furthermore, we develop an improved L2-$1_\sigma$ scheme with enhanced MBP preservation for large time steps, achieved through a novel unbalanced stabilization term that leverages the boundedness and monotonicity of the auxiliary functional. Representative numerical examples validate the accuracy, effectiveness, and physics-preserving of the proposed methods.

Tables

Table 1. Table 1: Time accuracy of the L L 1-sESAV scheme with α = 0.4 \alpha=0.4 and ι = 0.4 \iota=0.4

potential	$N$	$\gamma =2$		$\gamma =3$		$\gamma =4$ $\left(=\left(2-\alpha \right)/\iota \right)$
		$e\left(N\right)$	Order	$e\left(N\right)$	Order	$e\left(N\right)$	Order
	20	$5.06\times {10}^{-2}$	—	$1.68\times {10}^{-2}$	—	$1.25\times {10}^{-2}$	—
double-well	40	$2.91\times {10}^{-2}$	0.80	$7.64\times {10}^{-3}$	1.13	$4.35\times {10}^{-3}$	1.52
	80	$1.67\times {10}^{-2}$	0.80	$3.40\times {10}^{-3}$	1.17	$1.50\times {10}^{-3}$	1.54
	160	$9.59\times {10}^{-3}$	0.80	$1.48\times {10}^{-3}$	1.20	$5.07\times {10}^{-4}$	1.56
	20	$5.06\times {10}^{-2}$	—	$1.61\times {10}^{-2}$	—	$1.23\times {10}^{-2}$	—
Flory–Huggins	40	$2.91\times {10}^{-2}$	0.80	$7.50\times {10}^{-3}$	1.10	$4.12\times {10}^{-3}$	1.58
	80	$1.67\times {10}^{-2}$	0.80	$3.37\times {10}^{-3}$	1.15	$1.43\times {10}^{-3}$	1.52
	160	$9.59\times {10}^{-3}$	0.80	$1.48\times {10}^{-3}$	1.19	$4.91\times {10}^{-4}$	1.54
$\min \{2-\alpha ,\gamma \iota \}$			0.80		1.20		1.60

Table 2. Table 2: Time accuracy of the L L 1-sESAV scheme with α = 0.8 \alpha=0.8 and ι = 0.4 \iota=0.4

potential	$N$	$\gamma =2$		$\gamma =3$ $\left(=\left(2-\alpha \right)/\iota \right)$		$\gamma =4$
		$e\left(N\right)$	Order	$e\left(N\right)$	Order	$e\left(N\right)$	Order
	20	$1.38\times {10}^{-1}$	—	$1.00\times {10}^{-1}$	—	$1.05\times {10}^{-1}$	—
double-well	40	$7.91\times {10}^{-2}$	0.80	$4.59\times {10}^{-2}$	1.13	$4.72\times {10}^{-2}$	1.15
	80	$4.55\times {10}^{-2}$	0.80	$2.09\times {10}^{-2}$	1.14	$2.09\times {10}^{-2}$	1.18
	160	$2.61\times {10}^{-2}$	0.80	$9.28\times {10}^{-3}$	1.17	$9.15\times {10}^{-3}$	1.19
	20	$1.37\times {10}^{-1}$	—	$1.00\times {10}^{-1}$	—	$1.04\times {10}^{-1}$	—
Flory–Huggins	40	$7.91\times {10}^{-2}$	0.80	$4.58\times {10}^{-2}$	1.13	$4.71\times {10}^{-2}$	1.15
	80	$4.55\times {10}^{-2}$	0.80	$2.08\times {10}^{-2}$	1.14	$2.08\times {10}^{-2}$	1.18
	160	$2.61\times {10}^{-2}$	0.80	$9.27\times {10}^{-3}$	1.17	$9.14\times {10}^{-3}$	1.19
$\min \{2-\alpha ,\gamma \iota \}$			0.80		1.20		1.20

Table 3. Table 3: CPU times and the total number of time steps yielded by the L L 1-sESAV scheme

	double-well potential		Flory–Huggins potential
time-stepping strategy	$N$	CPU times	$N$	CPU times
uniform step $\tau =2$	281	4.81 s	281	4.82 s
adaptive step	577	11.17 s	598	12.05 s
uniform step $\tau =0.02$	25007	1 h 39 m 24 s	25007	1 h 40 m 24 s

Table 4. Table 4: Time accuracy of the two L ​ 2 L2 - 1 σ 1_{\sigma} -sESAV schemes with α = 0.8 \alpha=0.8 and γ = 2 / ι \gamma=2/\iota : the double-well potential

method	$N$	$\iota =0.3$		$\iota =0.5$		$\iota =0.8$
		$e\left(N\right)$	Order	$e\left(N\right)$	Order	$e\left(N\right)$	Order
	20	$9.74\times {10}^{-3}$	—	$2.67\times {10}^{-3}$	—	$3.35\times {10}^{-4}$	—
$L2$-$1_{\sigma }$-sESAV	40	$2.49\times {10}^{-3}$	1.97	$6.88\times {10}^{-4}$	1.96	$8.46\times {10}^{-5}$	1.98
scheme	80	$6.25\times {10}^{-4}$	2.00	$1.73\times {10}^{-4}$	1.99	$2.21\times {10}^{-5}$	1.94
	160	$1.59\times {10}^{-4}$	1.98	$4.33\times {10}^{-5}$	2.00	$5.54\times {10}^{-6}$	1.99
	20	$9.74\times {10}^{-3}$	—	$2.67\times {10}^{-3}$	—	$3.54\times {10}^{-4}$	—
unbalanced $L2$-$1_{\sigma }$	40	$2.49\times {10}^{-3}$	1.97	$6.88\times {10}^{-4}$	1.96	$8.46\times {10}^{-5}$	2.06
sESAV scheme	80	$6.25\times {10}^{-4}$	2.00	$1.73\times {10}^{-4}$	1.99	$2.21\times {10}^{-5}$	1.94
	160	$1.59\times {10}^{-4}$	1.98	$4.33\times {10}^{-5}$	2.00	$5.54\times {10}^{-6}$	1.99
$\min \{2,\gamma \iota \}$			2.00		2.00		2.00

Table 5. Table 5: Time accuracy of the two L ​ 2 L2 - 1 σ 1_{\sigma} -sESAV schemes with α = 0.8 \alpha=0.8 and γ = 2 / ι \gamma=2/\iota : the Flory–Huggins potential

method	$N$	$\iota =0.3$		$\iota =0.5$		$\iota =0.8$
		$e\left(N\right)$	Order	$e\left(N\right)$	Order	$e\left(N\right)$	Order
	20	$9.59\times {10}^{-3}$	—	$2.59\times {10}^{-3}$	—	$2.95\times {10}^{-4}$	—
$L2$-$1_{\sigma }$-sESAV	40	$2.49\times {10}^{-3}$	1.95	$6.82\times {10}^{-4}$	1.93	$8.26\times {10}^{-5}$	1.84
scheme	80	$6.25\times {10}^{-4}$	1.99	$1.73\times {10}^{-4}$	1.98	$2.18\times {10}^{-5}$	1.92
	160	$1.59\times {10}^{-4}$	1.98	$4.33\times {10}^{-5}$	2.00	$5.53\times {10}^{-6}$	1.98
	20	$9.64\times {10}^{-3}$	—	$2.59\times {10}^{-3}$	—	$2.95\times {10}^{-4}$	—
unbalanced $L2$-$1_{\sigma }$	40	$2.49\times {10}^{-3}$	1.95	$6.82\times {10}^{-4}$	1.93	$8.26\times {10}^{-5}$	1.84
sESAV scheme	80	$6.25\times {10}^{-4}$	1.99	$1.73\times {10}^{-4}$	1.98	$2.18\times {10}^{-5}$	1.92
	160	$1.59\times {10}^{-4}$	1.98	$4.33\times {10}^{-5}$	2.00	$5.53\times {10}^{-6}$	1.98
$\min \{2,\gamma \iota \}$			2.00		2.00		2.00

Equations

$${}_0^C{D}_t^{\alpha }\varphi =\mathcal{M}({\epsilon }^2\Delta \varphi +f(\varphi )),t>0,\mathbf{x}\in \Omega ;\varphi (\mathbf{x},0)={\varphi }_{\text{init}}(\mathbf{x}),\mathbf{x}\in \bar{\Omega },$$

$${}_0^C{D}_t^{\alpha }v={\int }_0^t{\omega }_{1-\alpha }(t-s){\partial }_{s}v(s)ds,{\omega }_{\mu }(t):=\frac{t^{\mu -1}}{\Gamma (\mu )}.$$

$$E[\varphi ](t)\le E[\varphi ](s),\forall t>s\text{with}E[\varphi ](t):={\int }_{\Omega }(\frac{{\epsilon }^2}{2}|\nabla \varphi (\mathbf{x}){|}^2+F(\varphi (\mathbf{x})))d\mathbf{x},$$

$$\underset{\mathbf{x}\in \bar{\Omega }}{\max }|{\varphi }_{\text{init }}(\mathbf{x})|\le \beta ⟹\underset{\mathbf{x}\in \bar{\Omega }}{\max }|\varphi (\mathbf{x},t)|\le \beta ,\forall t>0,$$

$$E[\varphi ](t)\le E[\varphi ](0),\forall t>0,$$

$$f(\varphi )>0,\forall \varphi \in (-\infty ,-\beta )\cap \text{Dom}(f);f(\varphi )<0,\forall \varphi \in (\beta ,\infty )\cap \text{Dom}(f).$$

$$F(\varphi )=\frac{1}{4}(1-{\varphi }^2{)}^2,f(\varphi )=-F^{\prime }(\varphi )=\varphi -{\varphi }^3,$$

$$F(\varphi )=\frac{\theta }{2}[(1+\varphi )\ln (1+\varphi )+(1-\varphi )\ln (1-\varphi )]-\frac{{\theta }_c}{2}{\varphi }^2,f(\varphi )=\frac{\theta }{2}\ln \frac{1-\varphi }{1+\varphi }+{\theta }_c\varphi ,$$

$${}_0^C{D}_t^{\alpha }v(t_{\ast })\le (\ge )0.$$

$$\underset{\mathbf{x}\in \bar{\Omega }}{\max }|\varphi (\mathbf{x},t)|\le \beta ,\forall t>0.$$

$${}_0^C{D}_t^{\alpha }\varphi ({\mathbf{x}}_{\ast },t_{\ast })\le 0.$$

$$\Delta \varphi ({\mathbf{x}}_{\ast },t_{\ast })\ge 0.$$

$$0\ge {}_0^C{D}_t^{\alpha }\varphi ({\mathbf{x}}_{\ast },t_{\ast })-\mathcal{M}{\epsilon }^2\Delta \varphi ({\mathbf{x}}_{\ast },t_{\ast })=\mathcal{M}f(\varphi ({\mathbf{x}}_{\ast },t_{\ast }))>0,$$

$$t_{\ast }(ϵ)=\max \{0<t\le T|\underset{\mathbf{x}\in \bar{\Omega }}{\max }|\varphi (\mathbf{x},t)|=\beta +ϵ\text{and}\underset{\mathbf{x}\in \bar{\Omega }}{\max }|\varphi (\mathbf{x},s)|<\beta +ϵ,\forall s\in [0,t)\}.$$

$${}_0^C{D}_t^{\alpha }\varphi ({\mathbf{x}}_{\ast }(ϵ),t_{\ast }(ϵ))\le 0\text{and}\Delta \varphi ({\mathbf{x}}_{\ast }(ϵ),t_{\ast }(ϵ))\ge 0.$$

$$0\ge {}_0^C{D}_t^{\alpha }\varphi ({\mathbf{x}}_{\ast }(ϵ),t_{\ast }(ϵ))-\mathcal{M}{\epsilon }^2\Delta \varphi ({\mathbf{x}}_{\ast }(ϵ),t_{\ast }(ϵ))=\mathcal{M}f(\varphi ({\mathbf{x}}_{\ast }(ϵ),t_{\ast }(ϵ)))>0,$$

$$\mathcal{E}[\varphi ,R]:=\frac{{\epsilon }^2}{2}{\int }_{\Omega }|\nabla \varphi (\mathbf{x}){|}^{2}d\mathbf{x}+R,$$

$${}_0^C{D}_t^{\alpha }\varphi =\mathcal{M}({\epsilon }^2\Delta \varphi +V(g(\varphi ,R))f(\varphi )),R_t=-V(g(\varphi ,R))\left(f(\varphi ),{\varphi }_t\right),$$

$${\mathbb{V}}_h:=\{v\mid v=\{v_{i,j}{\}}_{i,j=1}^M\text{and}v\text{is periodic}\}.$$

$$⟨v,w⟩=h^2\sum _{i,j=1}^M{v}_{ij}{w}_{ij},\Vert v\Vert =\sqrt{⟨v,v⟩},\Vert v{\Vert }_{\infty }=\underset{1\le i,j\le M}{\max }|v_{ij}|$$

$${\mathbb{V}}_{\beta }:=\left\{v\mid v\in {\mathbb{V}}_h\text{with}\Vert v{\Vert }_{\infty }\le \beta \right\}.$$

$${\Delta }_h{v}_{ij}=\frac{1}{h^2}\left(v_{i+1,j}+v_{i-1,j}+v_{i,j+1}+v_{i,j-1}-4v_{ij}\right),1\le i,j\le M,$$

$${\nabla }_h{v}_{ij}={\left(\frac{v_{i+1,j}-v_{ij}}{h},\frac{v_{i,j+1}-v_{ij}}{h}\right)}^{\top },1\le i,j\le M.$$

$${\mathbb{D}}_{\tau }^{\alpha }{v}^n=\sum _{k=1}^n{\int }_{t_{k-1}}^{t_k}\frac{{\omega }_{1-\alpha }(t_n-s)}{{\tau }_k}{\nabla }_{\tau }{v}^{k}ds=\sum _{k=1}^n{A}_{n-k}^{(n)}{\nabla }_{\tau }{v}^k,$$

$${\mathbb{D}}_{\tau }^{\alpha }{\hat{\varphi }}^1:=A_0^{(1)}({\hat{\varphi }}^1-{\varphi }^0)=\mathcal{M}({\epsilon }^2{\Delta }_h{\hat{\varphi }}^1+f({\hat{\varphi }}^1)),$$

$${\mathbb{D}}_{\tau }^{\alpha }{\varphi }^1=\mathcal{M}({\epsilon }^2{\Delta }_h{\varphi }^1+V(g_h({\hat{\varphi }}^1,R^0))f({\hat{\varphi }}^1)-\kappa V(g_h({\hat{\varphi }}^1,R^0))({\varphi }^1-{\hat{\varphi }}^1)),$$

$${\mathbb{D}}_{\tau }{R}^1=V(g_h(\widehat{{\varphi }^1},R^0))⟨-f({\hat{\varphi }}^1)+\kappa ({\varphi }^1-{\hat{\varphi }}^1),{\mathbb{D}}_{\tau }{\varphi }^1⟩.$$

$${\hat{\varphi }}^n=\min \left\{\max \left\{(1+r_n){\varphi }^{n-1}-r_n{\varphi }^{n-2},-\beta \right\},\beta \right\},$$

$${\mathbb{D}}_{\tau }^{\alpha }{\varphi }^n=\mathcal{M}({\epsilon }^2{\Delta }_h{\varphi }^n+V(g_h({\hat{\varphi }}^n,R^{n-1}))f({\hat{\varphi }}^n)-\kappa V(g_h({\hat{\varphi }}^n,R^{n-1}))({\varphi }^n-{\hat{\varphi }}^n)),$$

$${\mathbb{D}}_{\tau }{R}^n=V(g_h(\widehat{{\varphi }^n},R^{n-1}))⟨-f({\hat{\varphi }}^n)+\kappa ({\varphi }^n-{\hat{\varphi }}^n),{\mathbb{D}}_{\tau }{\varphi }^n⟩.$$

Full text

††journal: Journal of Computational Physics

1 Introduction

As a diffuse interface model, the classic phase-field model has been widely applied across various research areas, including material sciences Acta_Allen_1979 , JCP_CH_1958 , hydrodynamics ARFM_Anderson_1998 , PRE_Qian_2003 , biology and tumor growth JCP_Du_2004 , IJNMBE_Oden_2012 , JTB_Wise_2008 , CMAME_Huang_2024 . More recently, there has been growing interest in nonlocal phase-field models SISC_Tang_2019 , JCP_Wang_2017 , JDE_Akagi_2016 , SINUM_Ainthworth_2017 , SINUM_Du_2019 , PASMA_Inc_2018 , CMA_Wang_2019 . The incorporation of nonlocal operators into phase-field equations has been shown to substantially alter the diffusive dynamics. For example, in space-fractional Allen–Cahn phase field model, the fractional order can regulate the sharpness of the interface CMAME_Xu_2016 ; and in time-fractional phase-field models, the fractional order significantly influences the coarsening behavior JCP_Wang_2017 , SISC_Tang_2019 , CPC_Chen_2019 . Moreover, fractional models have been shown to provide more accurate descriptions of complex phenomena involving anomalous diffusion and memory effects than traditional integer-order models; see, e.g., JCP_Wang_2017 , SISC_Hou_2021 , PRL_Golding_2006 , Nature_Kirchner_2000 , PSS_Nigmatullin_1986 . In heterogeneous porous media, diffusive transport of solute particles is influenced by the strong interactions between fluids and solids, where a large quantity of solute particles may get absorbed to the solid formation MS_Sharma_2015 . As a result, the travel times of these adsorbed particles can differ substantially from those of particles that move freely within the bulk phase JCP_Zhokn_2017 , giving rise to anomalous subdiffusive transport. This regime is characterized by a sublinear growth of the particle mean square displacement $<r^{2}>$ with respect to the time $t$ PRL_Chepizhko_2013 . Such anomalous transport behavior can be naturally modeled using time-fractional diffusion models PR_Metzler_2000 , JMAA_Wang_2019 , where the fractional-order parameter quantifies the strength of subdiffusive behavior.

In this paper, we focus on the following phase-field model involving special time nonlocaity, namely the time-fractional Allen–Cahn (tFAC) equation

[TABLE]

where the spatial domain is denoted by $\Omega\subset\mathbb{R}^{d}$ with $d\leq 3$, and the Caputo fractional-order derivative operator ${}_{0}^{C}D^{\alpha}_{t}$ of order $\alpha\in(0,1)$ is defined as

[TABLE]

In (1.1), $\phi(\mathbf{x},t)$ is the unknown phase-field function, $\mathcal{M}$ is a positive constant mobility, $\varepsilon>0$ is the interaction length that describes the thickness of the transition boundary between materials, and $f(\phi)=-F^{\prime}(\phi)$ is a continuously differentiable nonlinear potential function.

It is well known that when $\alpha\rightarrow 1$, the tFAC equation recovers the classical Allen–Cahn equation, which satisfies the so-called energy dissipation law, i.e.,

[TABLE]

and the maximum-bound principle (MBP), i.e.,

[TABLE]

see SIREV_Du_2021 for details. Regarding the energy stability of the tFAC model, the authors SISC_Tang_2019 provided the first theoretical conclusion, i.e.,

[TABLE]

which differs from the classical energy dissipation law (1.2). Moreover, it was demonstrated in JSC_Du_2020 that the double-well tFAC equation also admits the MBP (1.3).

Another key feature of the tFAC equation is that its evolution process often requires a long time to reach a steady state and typically involves multiple time scales (see JCP_Liao_2020 , SISC_Liao_2021 , JSC_Liao_2024 ). Therefore, an adaptive time-stepping strategy SISC_2011_Qiao , JCP_Gomez_2011 , JCP_Wodo_2011 serves as a heuristic and available approach to enhance computational efficiency without compromising accuracy. Moreover, it is well known that solutions to subdiffusion problems, including the tFAC equation, exhibit weak singularities near the initial time, although they would be smooth away from $t=0$, see IMA_Jin_2016 , SINUM_Martin_2017 . The initial weak singularity may cause loss of convergence order in numerical simulations based upon a uniform temporal mesh. However, this issue can be alleviated by employing a special nonuniform temporal mesh, e.g., the graded mesh SINUM_Martin_2017 , SINUM_Liao_2018 . These two characteristics of the tFAC equation jointly motivate the development of effective and accurate nonuniform time-stepping numerical methods.

To achieve stable numerical simulations and avoid nonphysical solutions, it is hence necessary to develop numerical schemes that can preserve both the discrete energy stability and MBP for the tFAC model (1.1). In SISC_Tang_2019 , Tang et al. presented a first-order method combining the uniform $L1$ formula JCP_Xu_2007 and the stabilization technique SINUM_Xu_2006 , which preserves the energy stability (1.4) and MBP (1.3) at the discrete level. By discretizing the fractional derivative using backward Euler convolution quadrature, three energy-stable and MBP-preserving schemes of order $O(\tau^{\alpha})$ were proposed and analyzed in JSC_Du_2020 . To achieve high-order (more than first-order) and MBP-preserving numerical methods, the $L1$ scheme ACM_Ji_2020 , the $L1_{R}$ scheme SISC_Liao_2021 , the $L2$-$1_{\sigma}$ scheme JCP_Liao_2020 , JSC_Liao_2024 , and the SFTR (shifted fractional trapezoidal rule) scheme JSC_Huang_2023 were developed. However, there are basically two limitations to the above mentioned high-order numerical schemes: (i) they treat the nonlinear term either fully or partially implicit, which implies that their unique solvability is uncertain, especially for Flory–Huggins potential case, and a nonlinear iteration must be implemented at each time step; (ii) these schemes always require certain restrictions on the time-step to preserve the discrete MBP. These motivate us to develop a novel class of high-order linear schemes that are more efficient and can preserve the discrete MBP unconditionally. Although inherently nonlocal, the scalar auxiliary variable (SAV) method JCP_Shen_2018 , SINUM_Shen_2018 has been widely applied to phase-field models SISC_Hou_2021 , SINUM_Qiao_2022 , SISC_Zhao_2024 , JCP_Quan_2022 , SISC_Ji_2020 , owing to its ability to facilitate the design of linear and unconditionally ’modified’ energy-stable schemes. To make the modified energy closer to the original one, several improved SAV-type techniques have been proposed recently, including the relaxation methods JCP_Shen_2022 , JCP_Jiang_2022 and energy-optimized (EOP) strategies JSC_Liu_2024 , ANM_2025_Zhang . Moreover, due to the uncertainty in the signs of the nonlinear term coefficients (typically functions with respect to the SAV), the construction of MBP-preserving SAV-type schemes are nontrivial and compelling, particularly for high-order time discretization methods or time-fractional model problems SISC_Hou_2021 , JCP_Quan_2022 .

An inspiring approach is the combinations of prediction-correction method and stabilization technique proposed in JSC_Qiao_2022 , SINUM_Qiao_2022 , MOC_Hou_2023 for classical Allen–Cahn type equations. They first adopted a first-order scheme to provide a predicted solution with local second-order accuracy and preservation of the MBP, and then presented their second-order schemes by controlling the nonlinear term using an artificial stabilization term. Through this strategy, second-order linear schemes with conditional JSC_Qiao_2022 , MOC_Hou_2023 or unconditional SINUM_Qiao_2022 preservation of the MBP were successfully established for the classical Allen–Cahn equation. However, it should be noted that this approach is difficult to directly apply to the tFAC model (1.1), as a first-order method no longer provides sufficiently accurate predicted solutions, although this is feasible in the integer-order case. Thus, in the current work, we shall propose a general prediction strategy capable of providing satisfactory predicted solutions. We argue that it can also serve as a more efficient alternative approach for the prediction step in the numerical methods developed in JSC_Qiao_2022 , SINUM_Qiao_2022 , MOC_Hou_2023 .

In this paper, we aim to design a class of linear, energy-stable, and MBP-preserving high-order nonuniform time-stepping schemes for the tFAC equation (1.1), by leveraging the stabilized exponential SAV approach, namely the sESAV. To achieve this, we introduce a novel nonnegative auxiliary functional in Section 2, which serves twofold functions: (i) it ensures that the first-order approximation of the SAV does not compromise the temporal accuracy of $\phi$; and (ii) it enables the development of an effective stabilization term, allowing us to establish an MBP-preserving scheme. In combination with the newly developed prediction strategy and auxiliary functional, we first propose a linear stabilized $L1$-sESAV nonuniform time-stepping scheme, which possesses the following remarkable advantages:

it is unconditionally energy-stable;

• 2.

it can unconditionally preserve the discrete MBP, which seems to be the first time such a linear, MBP-preserving, and $(2-\alpha)$th order $L1$ scheme is established;

• 3.

it is computationally efficient, in fact, only systems of linear algebraic equations require to be solved for $n\geq 2$.

Then, a linear stabilized $L2$-$1_{\sigma}$-sESAV scheme is also investigated and proven to be unconditionally energy-stable, but conditionally MBP-preserving. Moreover, to further improve the $L2$-$1_{\sigma}$-sESAV scheme for large time-step computation, we present an unbalanced $L2$-$1_{\sigma}$-sESAV scheme by introducing a novel stabilization term based on the boundedness and monotonicity of the novel auxiliary functional. To the best of our knowledge, so far no such linear second-order scheme with provable energy stability and MBP-preserving in the discrete setting has been developed. Furthermore, we demonstrate that the proposed method can be coupled with the EOP technique to enhance the modified energy, while simultaneously retaining the discrete energy stability and the discrete MBP; see Remark 5 and Example 3.2.

The remainder of the paper is organized as follows. In Section 2, we first prove the MBP for the tFAC equation under either double-well or Flory–Huggins potential. Then, we propose and analyze the fully-discrete linear stabilized $L1$ type scheme, including the unique solvability, unconditional MBP-preservation and energy-stability. Numerical experiments are conducted in Section 3 to illustrate the performance of the proposed method. In Section 4, we develop two linear stabilized $L2$-$1_{\sigma}$-ESAV schemes, incorporating two different stabilization terms: one balanced and one unbalanced, and unconditional energy-stability and conditional MBP-preservation are also discussed and tested. Some concluding remarks are finally drawn in the last section.

2 Maximum-bound principle and stabilized $L1$ scheme

For simplicity, the whole discussion below is confined to the two-dimensional (2D) square domain $\Omega=(0,L)^{2}$ with periodic boundary conditions. It is worth noting that the extensions to the three-dimensional (3D) case and/or homogeneous Neumann boundary conditions present no significant difficulties. In this section, we first discuss the MBP for the tFAC equation (1.1). Assume that there exists a positive constant $\beta$ such that $f(\beta)=f(-\beta)=0$ and $F(\phi)$ satisfies the monotone conditions away from $(-\beta,\beta)$, i.e.,

[TABLE]

Two such typical potentials are widely used in the classical Allen–Cahn equation: one is the double-well potential

[TABLE]

where $\text{Dom}(f)=\mathbb{R}$ and $\beta=1$, and the other is the Flory–Huggins potential

[TABLE]

with $\theta_{c}>\theta>0$, where $\text{Dom}(f)=(-1,1)$ and $\beta$ is the positive root of $f(\rho)=0$. For instance, if we set $\theta=0.8$ and $\theta_{c}=1.6$, we have $\beta\approx 0.9575$.

The following lemma is originally presented in Lemma 2.3 of Ref. JSC_Du_2020 , we provide a detailed proof in A for completeness.

Lemma 2.1**.**

Assume that $v\in C[0,T]\cap C^{1}(0,T]$ and attains its minimum (maximum) at $t_{*}\in(0,T]$. Then there holds

[TABLE]

Then, based on Lemma 2.1, we now proceed to establish the main MBP conclusion for the tFAC equation (1.1).

Theorem 2.2**.**

Suppose that $\phi\in C([0,T];C^{2}(\bar{\Omega}))\cap C^{1}((0,T],C(\bar{\Omega}))$ is a solution to the model problem (1.1) with double-well or Flory–Huggins potential. If the initial value satisfies $\max_{\mathbf{x}\in\bar{\Omega}}|\phi_{\text{init }}(\mathbf{x})|\leq\beta$, there holds

[TABLE]

Proof.

The conclusion for the double-well potential with $\beta=1$ can be found in [JSC_Du_2020, , Theorem 2.3]. For completeness, we briefly revisit its proof here. Assume that the maximum of $\phi$ is smaller than $-1$ and is achieved at $(\mathbf{x}_{*},t_{*})\in\Omega\times(0,T]$. Then, by the regularity of $\phi$ and Lemma 2.1, we have

[TABLE]

Moreover, owing to the periodic boundary conditions, the following estimate holds by Taylor expansion

[TABLE]

Thus, we obtain

[TABLE]

which is a contradiction. The upper bound of the solution can be proved analogously.

It remains to prove the MBP (2.4) for the Flory–Huggins potential (2.3). To this aim, for any $\epsilon\in(0,\frac{1-\beta}{2})$ such that $\beta+\epsilon<1$, we define

[TABLE]

It is clear that $t_{*}(\epsilon)$ is well-defined by the continuity assumption on $\phi$. Without loss of generality, assume that $\phi(\mathbf{x},t_{*}(\epsilon))$ achieves the maximum absolute value at $\mathbf{x}_{*}(\epsilon)\in\bar{\Omega}$ and $\phi(\mathbf{x}_{*}(\epsilon),t_{*}(\epsilon))=-\beta-\epsilon$. Similar to (2.5)–(2.6), we have

[TABLE]

Therefore, we arrive at the result

[TABLE]

where we have used (2.1) in the last inequality. This contradiction shows that no point $(\mathbf{x}_{*}(\epsilon),t_{*}(\epsilon))\in\bar{\Omega}\times(0,T]$ satisfies $\phi(\mathbf{x}_{*}(\epsilon),t_{*}(\epsilon))=-\beta-\epsilon$ for any given $T$ and $\epsilon$. A similar conclusion holds for the upper bound of $\phi$. Thus, the MBP (2.4) is proved for the Flory–Huggins potential, as both $T>0$ and $\epsilon\in(0,\frac{1-\beta}{2})$ are arbitrary. ∎

Introduce an auxiliary variable $R(t)=E_{1}[\phi]:=\int_{\Omega}F(\phi)d\mathbf{x}$, and define the modified energy

[TABLE]

which is equivalent to the original energy $E[\phi]$ defined in (1.2) at the continuous level. Moreover, let $g(\phi,R):=\frac{\exp\left\{R\right\}}{\exp\left\{E_{1}[\phi]\right\}}$, which is always equal to 1 at the continuous level. Then, model (1.1) can be rewritten into the following equivalent system:

[TABLE]

where $V(\cdot)$ is an auxiliary functional satisfying the following three assumptions:

• (A1).

The auxiliary functional $V(\cdot)\in C^{1}(\mathbb{R})\cap W^{2,\infty}(\mathbb{R})$ such that $V(1)=1$, $V^{\prime}(1)=0$, and there exists a positive constant $K_{1}$ such that $|V^{\prime}\left(\cdot\right)|\leq K_{1}$;

• (A2).

The auxiliary functional is nonnegative and bounded by a positive constant $K_{2}$, i.e., $0\leq V\left(\cdot\right)\leq K_{2}$;

• (A3).

The auxiliary functional satisfies the following monotonicity, i.e., for any $z_{1},z_{2}\in\mathbb{R}$ satisfy $|z_{1}-1|\leq|z_{2}-1|$, we have $|V(z_{1})-1|\leq|V(z_{2})-1|$.

By [SISC_Tang_2019, , Corollary 2.1], it is straightforward to verify that the new SAV system (2.8) satisfies energy stablity with respect to $\mathcal{E}[\phi,R]$, i.e., $\mathcal{E}[\phi,R](t)\leq\mathcal{E}[\phi,R](0)$ for $t>0$.

2.1 Stabilized $L1$-ESAV scheme

Given a positive integer $M$, let $h=L/M$ be the spatial grid length and set $\Omega_{h}:=\{\mathbf{x}_{h}=(ih,jh)\mid 0\leq i,j\leq M\}$. Let $\mathbb{V}_{h}$ be the set of all $M$-periodic real-valued grid functions on $\Omega_{h}$, i.e.,

[TABLE]

Notice that $\mathbb{V}_{h}$ is a finite-dimensional linear space, thus any grid function in $\mathbb{V}_{h}$ and any linear operator $P:\mathbb{V}_{h}\rightarrow\mathbb{V}_{h}$ can be treated as a vector in $\mathbb{R}^{M^{2}}$ and a matrix in $\mathbb{R}^{M^{2}\times M^{2}}$, respectively.

Define the following discrete inner product, and discrete $L^{2}$ and $L^{\infty}$ norms

[TABLE]

for any $v,w\in\mathbb{V}_{h}$. Moreover, we define a maximum-norm function space with the positive constant $\beta$ for the underlying MBP problem, that is

[TABLE]

Next, we consider a central finite difference discretization for the spatial differential operators. For any $v\in\mathbb{V}_{h}$, the discrete Laplace operator $\Delta_{h}$ is defined by

[TABLE]

and the discrete gradient operator $\nabla_{h}$ is defined by

[TABLE]

Consider (generally nonuniform) time levels $0=t_{0}<t_{1}<\cdots<t_{N}=T$ with time steps $\tau_{k}:=t_{k}-t_{k-1}$ for $1\leq k\leq N$. Denote by $r_{k}:=\tau_{k}/\tau_{k-1}(k\geq 2)$ the adjacent time-step ratio. Let $v^{k}=v\left(t_{k}\right)$ and the temporal difference operator $\nabla_{\tau}v^{k}=v^{k}-v^{k-1}$ and $\mathbb{D}_{\tau}v^{k}=\nabla_{\tau}v^{k}/\tau_{k}$ for $k\geq 1$. The nonuniform $L1$ formula of the Caputo derivative is given by

[TABLE]

where the time-dependent convolution kernels $A_{n-k}^{(n)}:=\frac{1}{\tau_{k}}\int_{t_{k-1}}^{t_{k}}{\omega_{1-\alpha}(t_{n}-s)}ds$ satisfy the following two main properties (see SINUM_Liao_2018 , SCM_Liao_2024 ):

(P1). The discrete kernels are positive, i.e., $A_{n-k}^{(n)}>0$ for $1\leq k\leq n$;

(P2). The discrete kernels are monotone, i.e., $A_{n-k+1}^{(n)}<A_{n-k}^{(n)}$ for $1\leq k\leq n$.

Now, let $\{\phi^{n},R^{n}\}$ be the fully discrete approximations of the exact solutions $\{\phi(t_{n}),R(t_{n})\}$ to the original continuous problem (2.8). We apply the nonuniform $L1$ formula (2.10) and central difference approximation (2.9) to derive the following stabilized $(2-\alpha)$th order ESAV scheme (denoted as $L1$-sESAV):

Step 1

Let $\hat{\phi}^{0}=\phi^{0}:=\phi_{\text{init }}$. For $n=1$, first solve a predicted solution $\hat{\phi}^{1}$ from the fully-implicit difference scheme:

[TABLE]

and then, find $\{\phi^{1},R^{1}\}\in\mathbb{V}_{h}\times\mathbb{R}$ by

[TABLE]

[TABLE]

where $\kappa\geq 0$ is a stabilizing constant, and $g_{h}$, a discrete version of $g$, denoted by $g_{h}(v,w):=\frac{\exp\left\{w\right\}}{\exp\left\{E_{1h}[v]\right\}}$ with $E_{1h}$ the discrete counterpart of $E_{1}$, i.e., $E_{1h}[v]:=\langle F(v),1\rangle$.

Step 2

For $n\geq 2$, given $\phi^{n-1}$ and $\phi^{n-2}$, first predict a MBP-preserving solution $\hat{\phi}^{n}$ by pre-processing the standard linear extrapolation, i.e.,

[TABLE]

and then, find $\{\phi^{n},R^{n}\}\in\mathbb{V}_{h}\times\mathbb{R}$ such that

[TABLE]

[TABLE]

Remark 1*.*

In the peoposed $L1$-sESAV scheme (2.11)–(2.16), the assumption (A1) on the functional $V(\cdot)$ ensures that the first-order approximation of the auxiliary variable $R$ does not affect the temporal convergence of the phase-field function $\phi$. Actually, a direct application of Taylor expansion, together with the assumption (A1), gives us

[TABLE]

which implies that if $z$ is a first-order approximation to $1$, then $V(z)$ will be a second-order approximation to $1$. Moreover, the non-negativity and boundedness assumption (A2) allows us to develop effective stabilization term to dominate the nonlinear potential term, leading to MBP-preserving numerical methods, see Theorems 2.9 and 4.5 for details.

Remark 2*.*

To effectively solve the nonlinear scheme (2.11) in Step 1 for $n=1$, we propose the following simple iteration scheme with stabilization:

[TABLE]

where the subscript $\mathfrak{s}$ denotes the iteration index, and the initial iteration value is taken as $\hat{\phi}^{1}_{(0)}=\phi^{0}$. In the next section, we will show that the iteration scheme (2.17) unconditionally preserves the discrete MBP, and meanwhile, the solution $\hat{\phi}^{1}_{(\mathfrak{s})}$ converges to the unique solution $\hat{\phi}^{1}$ of the nonlinear scheme (2.11). The detailed implementation of the $L1$-sESAV scheme is summarized as below.

In the end of this subsection, we discuss the linear extrapolation error $\mathcal{R}^{n}[v]:=v(t_{n})-\hat{v}(t_{n})$ for $n\geq 2$ under the nonuniform mesh, where $\hat{v}(t_{n}):=(1+r_{n})v(t_{n-1})-r_{n}v(t_{n-2})$ denotes the linear extrapolation. It is well known that if $v$ is sufficiently smooth, the extrapolation can achieve second-order temporal accuracy. However, for the tFAC model that exhibits an initial weak singularity, the temporal accuracy may deteriorate. The following lemma provides an estimate of $\mathcal{R}^{n}[v]$, see also B for detailed proof.

Lemma 2.3**.**

Assume that $v\in C^{2}(0,T]$ and $|v^{\prime\prime}(t)|\leq C(1+t^{\iota-2})$ for any regularity parameter $0<\iota<1$. Then, there exists a constant $C_{v}$ that depends only on $v$ such that

[TABLE]

Moreover, if the graded temporal grids $t_{k}=T(k/N)^{\gamma}$ with grading parameter $\gamma\geq 1$ are employed, it follows that

[TABLE]

where $C_{v,\gamma}$ is a constant depending only on $v$ and $\gamma$.

Remark 3*.*

It is well known that under the initial weak singularity assumption

[TABLE]

the global consistency error, which is commonly used in error analysis, attained by the $L1$ and $L2$-$1_{\sigma}$ formulas achieves the order of $\min{\{2-\alpha,\gamma\iota\}}$ and $\min{\{2,\gamma\iota\}}$ on graded temporal grids, respectively; see SINUM_Liao_2018 , CiCP_Liao_2021 for more details. Together with the interpolation error $\mathcal{R}^{n}[v]$ established in Lemma 2.3, this indicates that, to compensate for the accuracy reduction caused by the weak singularity of the solution to the tFAC model (1.1), it is sufficient to select the mesh grading parameter $\gamma=(2-\alpha)/\iota$ for the $L1$ scheme and $\gamma=2/\iota$ for the $L2$-$1_{\sigma}$ scheme introduced in Section 4. The numerical results presented in Examples 3.1 and 4.4.1 have confirmed the validity of this statement.

2.2 MBP-preservation and unique solvability of the prediction strategy

In this subsection, we shall show that the predicted solution $\hat{\phi}^{n}$ yielded by Step 1 or Step 2 is MBP-preserving and uniquely solvable. From (2.14), it can be seen that the conclusion is clearly valid for $n\geq 2$. It remains to consider the case $n=1$. The following two lemmas will be used later.

Lemma 2.4** (JCM_Tang_2016 ).**

For any $\lambda>0$, we have $\|(\lambda I-\Delta_{h})^{-1}\|_{\infty}\leq\lambda^{-1}$, where $I$ represents the identity operator.

Lemma 2.5** (SIREV_Du_2021 ).**

If $\kappa\geq\|f^{\prime}\|_{C[-\beta,\beta]}$ holds for some positive constant $\kappa$, then we have $|f(\xi)+\kappa\xi|\leq\kappa\beta$ for any $\xi\in[-\beta,\beta]$.

Next, we use the well-known Banach’s fixed point theorem to prove that the nonlinear scheme (2.11) is uniquely solvable and preserves the discrete MBP.

Theorem 2.6**.**

If $\kappa\geq\|f^{\prime}\|_{C[-\beta,\beta]}$, the simple iterative scheme (2.17) unconditionally preserves the MBP for $\{\hat{\phi}^{1}_{(\mathfrak{s})}\}$, that is,

[TABLE]

Furthermore, if the first level time-step $\tau_{1}$ is sufficiently small such that

[TABLE]

then (i) the iterative scheme (2.17) converges to the unique solution $\hat{\phi}^{1}$ of the nonlinear scheme (2.11) in the maximum norm such that $\|\hat{\phi}^{1}\|_{\infty}\leq\beta$, and (ii) the $L1$-sESAV scheme (2.11)–(2.16) is uniquely solvable.

Proof.

For any $\mathfrak{s}\geq 1$, we assume that $\|\hat{\phi}^{1}_{(\mathfrak{s}-1)}\|_{\infty}\leq\beta.$ Rewrite (2.17) equivalently as

[TABLE]

which together with Lemmas 2.4–2.5 and $\|\phi^{0}\|_{\infty}\leq\beta$ yields

[TABLE]

Thus, the conclusion (2.21) is proved.

Next, we discuss the MBP-preservation and unique solvability of the nonlinear scheme (2.11). To this aim, for any grid function $v\in\mathbb{V}_{\beta}$, we define the mapping $\mathcal{T}:=\mathbb{V}_{\beta}\rightarrow\mathbb{V}_{\beta}$ through (2.17) by $\mathcal{T}[v]=w\in\mathbb{V}_{\beta}$, i.e.,

[TABLE]

We shall show that $\mathcal{T}$ is a contraction mapping under the time-step restriction (2.22) in what follows. For any $v_{1},v_{2}\in\mathbb{V}_{\beta}$, let $w_{1}=\mathcal{T}[v_{1}],w_{2}=\mathcal{T}[v_{2}]$, and define $\delta w=w_{1}-w_{2}$. Then it follows from (2.23) that

[TABLE]

Further, Lemma 2.4 and the condition $\kappa\geq\|f^{\prime}\|_{C[-\beta,\beta]}$ gives us

[TABLE]

As a consequence, it holds that

[TABLE]

Since $A_{0}^{(1)}=\frac{1}{\tau_{1}}\int_{t_{0}}^{t_{1}}\omega_{1-\alpha}(t_{1}-s)ds=\frac{\tau_{1}^{-\alpha}}{\Gamma(2-\alpha)}$, we conclude that if $\tau_{1}\leq\sqrt[\alpha]{1/(\kappa\mathcal{M}\Gamma(2-\alpha))}$, it follows that $\ell_{\mathcal{T}}<1$, ensuring that $\mathcal{T}$ is a contractive mapping. By Banach’s fixed point theorem, the mapping $\mathcal{T}$ admits a unique fixed point $\hat{\phi}^{1}\in\mathbb{V}_{\beta}$. Consequently, the iterative scheme (2.17) is convergent in the maximum norm, and the fixed point $\hat{\phi}^{1}$ is just the unique solution of the nonlinear scheme (2.11).

For the simplicity of presentation, below we denote $V^{n}:=V(g_{h}(\hat{\phi}^{n},R^{n-1}))$ for $n\geq 1$. Then, equations (2.12)–(2.13) and (2.15)–(2.16) can be equivalently and uniformly rewritten as:

[TABLE]

[TABLE]

for $n\geq 1$. It follows from $A^{(n)}_{0}>0$ and assumption (A2) that the coefficient matrix of (2.25) is symmetric and positive definite, and thus (2.25)–(2.26) is uniquely solvable, which further implies the uniquely solvablity of the $L1$-sESAV scheme (2.11)–(2.16). The proof is completed. ∎

Remark 4*.*

Note that a useful estimate (2.24) for the contraction factor $\ell_{\mathcal{T}}$ is provided in the proof of Theorem 2.6, indicating that the iteration error decreases by at least a factor of $\ell_{\mathcal{T}}$ per iteration. Let the first step iteration error be denoted by $e_{(1)}:=\|\hat{\phi}^{1}_{(1)}-\phi^{0}\|_{\infty}\leq 2\beta$. Then, for a given termination tolerance $tol$, a maximum number of iterations $\mathfrak{s}^{*}=\big\lceil\frac{\ln(tol/e_{(1)})}{\ln(\ell_{\mathcal{T}})}\big\rceil+1\leq\big\lceil\frac{\ln(tol/(2\beta))}{\ln(\ell_{\mathcal{T}})}\big\rceil+1$ is sufficient to ensure $\|\hat{\phi}^{1}_{(\mathfrak{s}^{*}+1)}-\hat{\phi}^{1}_{(\mathfrak{s}^{*})}\|_{\infty}\leq tol$. Moreover, a smaller stabilization constant $\kappa$ leads to faster convergence of the simple iterative scheme (2.17).

2.3 Discrete energy stability and MBP of the $L1$-sESAV scheme

Define the discrete version energy that corresponding to (2.7):

[TABLE]

Next, we will establish the energy stability of the $L1$-sESAV scheme.

Lemma 2.7** (SCM_Liao_2024 ).**

For fixed $n\geqslant 1$, the convolution kernels $\{A_{n-k}^{(n)}\}$ are positive semi-definite in the sense that

[TABLE]

Theorem 2.8**.**

Assume that the conditions in Theorem 2.6 hold, the $L1$-sESAV scheme (2.11)–(2.16) is unconditionally energy-stable in the sense that $\mathcal{E}_{h}[\phi^{n},R^{n}]\leq\mathcal{E}_{h}[\phi^{0},R^{0}]$ for $n\geq 1$.

Proof.

Let $n=k$ in (2.15), and then taking the inner product with $\nabla_{\tau}\phi^{k}=\phi^{k}-\phi^{k-1}$, we have

[TABLE]

which, together with (2.26), leads to

[TABLE]

for $k\geq 2.$ Similarly, it follows from (2.12) and (2.26) that (2.27) also holds for $k=1$.

Now, sum the above inequality from $k=1$ to $n$, we obtain from (2.10) and Lemma 2.7 that

[TABLE]

which completes the proof. ∎

In the following theorem, we are ready to establish the discrete MBP for the proposed $L1$-sESAV scheme.

Theorem 2.9**.**

Assume that the conditions in Theorem 2.6 hold, the $L1$-sESAV scheme (2.11)–(2.16) unconditionally preserves the MBP for $\left\{\phi^{n}\right\}$, that is,

[TABLE]

Proof.

This claim will be verified by the complete mathematical induction argument. Assume that $\|\phi^{k}\|_{\infty}\leq\beta$ for $0\leq k\leq n-1$. From (2.14) and Theorem 2.6, we know that $\|\hat{\phi}^{n}\|_{\infty}\leq\beta$, which together with (2.25) and Lemmas 2.4–2.5 leads to

[TABLE]

where the assumption (A2) has been used. Due to the positivity and monotonicity of the discrete kernels $\{A^{(n)}_{n-k}\}$, it holds that

[TABLE]

Therefore, inserting (2.29) into (2.28) yields the conclusion for $k=n$, and thereby the proof is completed. ∎

Remark 5*.*

The linear $L1$-sESAV scheme developed in this section for the tFAC model (1.1) is MBP-preserving. However, as noted in the introduction, the SAV approach only guarantees a modified form of energy stability. To further reduce the gap between the modified and original energies, motivated by the EOP strategy JSC_Liu_2024 , an improved strategy, referred to as the $L1$-sESAV-EOP scheme, is proposed as follows: Let $\tilde{R}^{n}$ be the solution of (2.16) in Step 2, the auxiliary variable $R^{n}$ is then updated by

[TABLE]

which is a solution to the following linear optimization problem

[TABLE]

It is clear that the resulting $L1$-sESAV-EOP scheme remains energy-stable with respect to a modified energy, which, notably, closely approximates the original energy as $\eta_{0}\in[0,1]$. In particular, when $R^{n}=E_{1h}[\phi^{n}]$ in (2.30), we have $\eta_{0}=0$ in (2.31). In this case, the modified energy coincides exactly with the original energy. Moreover, by the same argument as in Theorem 2.9, the improved scheme also preserves the discrete MBP.

3 Numerical experiments

In this section, ample numerical tests are provided to illustrate the accuracy and structure-preserving properties of the $L$1-sESAV scheme for the tFAC model (1.1). Note that the auxiliary functional introduced in Section 2 plays a crucial role in the construction and analysis of the MBP-preserving $L$1-sESAV scheme. In the numerical implementation, we choose the following piecewise polynomial function that satisfies the assumptions (A1)–(A3) with $K_{1}=49/24$ and $K_{2}=1$:

[TABLE]

In all the following tests, we always set the stabilization constant $\kappa=\|f^{\prime}\|_{C[-1,1]}=2$ for the double-well potential (2.2), and $\kappa=\|f^{\prime}\|_{C[-\beta,\beta]}\approx 8.02$ for the Flory–Huggins potential (2.3), for which the parameters are chosen as $\theta=0.8$ and $\theta_{c}=1.6$.

Moreover, to address the initial weak singularity, we adopt two types of mixed nonuniform temporal meshes composed of a graded mesh in the initial subinterval. Specifically, for given $\hat{T}$ and $\hat{N}$, we split the time interval $[0,T]$ into two parts $[0,\hat{T}]$ and $[\hat{T},T]$ with total $N$ subintervals, and employ the graded mesh to subdivide $[0,\hat{T}]$ as follows

[TABLE]

where $\gamma\geq 1$ is the mesh grading parameter. For the remainder interval $[\hat{T},T]$, the following two different types of temporal girds are considered

uniform mesh:

[TABLE]

• 2.

adaptive mesh based on the energy variation SISC_Liao_2021 , SISC_2011_Qiao :

[TABLE]

where $\tau_{\max},\tau_{\min}$ are predetermined maximum and minimum time steps, and $\eta$ is a tunable parameter.

Remark 6*.*

The energy-based adaptive time-stepping strategy (3.4) employed in this work adjusts time steps based on variations in the discrete energy SISC_Liao_2021 , SISC_2011_Qiao , which significantly enhances the computational efficiency for long-term simulations as shown in Table 3. In practice, the error-based adaptive time-stepping strategy JCP_Wodo_2011 , JCP_Gomez_2011 , JCP_Liao_2020 also offers an effective alternative for the efficient simulation of phase-field models. Note that both strategies inherently generate nonuniform temporal meshes. Since this paper focuses on the development and analysis of MBP-preserving numerical schemes for general nonuniform time grids, our theoretical results provide a unified foundation to ensure their reliable application.

Furthermore, owing to the well-structured profile of the phase variable—taking distinct (close-to-constant) values in the bulk phases and exhibiting large gradients only within the thin diffuse interface JSC_Feng_2005 —the development and analysis of numerical schemes on nonuniform spatial meshes are also highly meaningful, as they furnish theoretical support for space adaptivity JCP_Provatas_1999 , JCP_Mackenzie_2002 , JCP_Feng_2006 . In addition, combining the proposed numerical method with $h$-refinement and local time-stepping (LTS) techniques—where each grid point is allowed to adopt its own locally time steps SISC_Fernando_2022 , SISC_Mehlin_2015 —may further enhance computational efficiency and accuracy. Nevertheless, preserving key physical properties, such as the discrete energy stability and the discrete MBP, within nonuniform spatial meshes and LTS frameworks remains highly challenging and requires further investigation in future.

In the following numerical tests, the error is measured in the spatio-temporal maximum-norm that $e(N):=\max_{1\leq n\leq N}\|\phi(t_{n})-\phi^{n}\|_{\infty}$.

3.1 Temporal convergence

Consider the exterior-forced tFAC model

[TABLE]

with parameters $\mathcal{M}=0.01$ and $\varepsilon=1$. Here, the linear part $g(\mathbf{x},t)$ is chosen such that the exact solution is given by $\phi(\mathbf{x},t)=\omega_{1+\iota}(t)\sin(x)\sin(y)$ with parameter $\iota\in(0,1)$. It is clear that the solution exhibits a weak singularity as indicated in (2.20).

For a given mesh grading parameter $\gamma\geq 1$, we set $\hat{T}=\min\{1/\gamma,T\}$ and $\hat{N}=\lceil\frac{1}{T+1+\gamma^{-1}}\rceil$ in (3.2). The uniform temporal mesh (3.3) is employed in the remainder interval $[\hat{T},T]$. The spatial domain is partitioned into $M=400$ uniform elements along each spatial direction. We test the proposed $L$1-sESAV scheme for solving (3.5) with different nonlinear potentials and a fixed regularity parameter $\iota=0.4$. For fractional orders $\alpha=0.4$ and $0.8$, different mesh grading parameters $\gamma=2,3,4$ are considered, respectively. The corresponding numerical results are presented in Tables 1–2, from which we can observe that the expected temporal accuracy $\mathcal{O}(N^{-\min\{2-\alpha,\gamma\iota\}})$ is achieved.

3.2 Unconditional preservation of MBP and energy stability

In this subsection, we numerically verify the discrete MBP and energy stability of the proposed $L$1-sESAV scheme by simulating the spinodal decomposition of a homogeneous mixture into two coexisting phases governed by the tFAC model (1.1) with $\alpha=0.5$, $\mathcal{M}=1$, $\varepsilon=0.01$ in $\Omega=(0,1)^{2}$. For the forthcoming simulation, we set $M=128$, and the initial phase field is generated using uniformly distributed random data in the range $[-0.8,0.8]$, which is highly oscillating, see Figure 1.

Although a modified energy is introduced as an approximation of the original one to facilitate the proof of energy stability, we primarily focus on the behavior of the original (discrete) energy, which reflects the real physical mechanism of the dynamic process. In Figure 2, we show the time evolutions of the maximum norms and original energies of the numerical solutions computed using the proposed $L$1-sESAV scheme with various time steps $\tau=2,0.2,0.02,0.002$ for the case of double-well potential. It shows that the $L$1-sESAV scheme perfectly preserves both the MBP and the energy stability, even for large time steps (e.g., $\tau=2$ and $0.2$), which indicates the unconditional preservation of the discrete energy stability and MBP as proved in Theorems 2.8–2.9. Next, we also consider the simulation of Flory–Huggins potential case and corresponding numerical results are displayed in Figure 3. Similar conclusions can also be observed.

Furthermore, we compare the modified and original energies yielded by the $L1$-sESAV scheme and $L1$-sESAV-EOP scheme, respectively, with $\tau=0.002$ for the double-well potential. As seen from Figure 4, a discrepancy exists between the modified energy and the original energy for the original $L1$-sESAV scheme, however, by applying the EOP technique introduced in Remark 5, this discrepancy is substantially reduced for the $L1$-sESAV-EOP scheme. Moreover, Figure 5 demonstrates that the $L1$-sESAV-EOP scheme continues to preserve both the discrete energy stability and MBP.

3.3 Application of adaptive time-stepping strategy

It is well known that the dynamics process governed by the tFAC model involves multiple time scales, and typically requires a long time to reach the steady state. Hence, an adaptive time-stepping strategy (e.g., (3.4)) is an effective means for improving the computational efficiency without sacrificing accuracy. In our computations, the graded mesh with $\hat{T}=0.5$, $\hat{N}=30$ and $\gamma=(2-\alpha)/\alpha$ is employed in the initial interval $[0,\hat{T}]$, and the remainder $[\hat{T},T]$ is partitioned using three different types of temporal meshes, i.e., uniform large time step $\tau=2$ and small time step $\tau=0.02$, and adaptive time step based on (3.4), with the parameters $\tau_{\max}=2,\tau_{\min}=0.02,\eta=10^{6}$. The other settings are the same as those used in Example 3.2.

Figure 6 presents a comparison of the solution snapshots obtained using three types of temporal meshes up to $T=500$ for the model with double-well potential. It is observed that using a large uniform time step $\tau=2$ results in an inaccurate solution $\phi$, while the adaptive time-stepping strategy yields coarsening patterns consistent with those obtained from the small uniform time step $\tau=0.02$. Besides, the time evolutions of the maximum-norm of $\phi$ and the energy are depicted in Figures 7(a)–7(b), from which we can observe that (i) the proposed $L$1-sESAV scheme is energy-stable and MBP-preserving; (ii) the evolution curves of the maximum-norm of $\phi$ and energy obtained by the adaptive strategy match very well with those generated by the small uniform time step. Moreover, Table 3 indicates the high efficiency of the proposed scheme combined with the adaptive strategy (3.4). For instance, it takes more than one and a half hours for the implementation with uniform time step $\tau=0.02$, whereas the adaptive method requires only about 12 seconds! In fact, the total number of adaptive time steps is only 577, which is significantly fewer than in the small uniform time step case. This improvement is mainly due to the frequent use of large time steps in the adaptive method. As shown in Figure 7(c), only a few small time steps are employed when the energy decays rapidly. For the Flory–Huggins potential, similar conclusions can also be drawn from Figures 8–9 and Table 3.

3.4 Dynamical behaviour governed by the tFAC model

In order to discover the influence of fractional order $\alpha$ on the dynamical behaviour, we consider the 2D tFAC model (1.1) with $\mathcal{M}=1$, $\varepsilon=0.01$, and double-well potential. The initial state is chosen as

[TABLE]

The computational domain $\Omega=(0,1)^{2}$ is uniformly divided into 128 elements along each spatial direction.

In the adaptive time-stepping strategy (3.4), we set $\tau_{\max}=2$, $\tau_{\min}=0.02$, and the parameter $\eta=10^{7}$. The time evolutions of the phase-field function with different fractional orders $\alpha=0.9,0.7$ and $0.4$ at different time instants are shown in Figure 10. It is clearly seen that the dynamical behaviour is significantly affected by the fractional order: the bigger the $\alpha$ is, the faster the dynamics evolutes. Thus, smaller $\alpha$ requires a longer time to reach equilibrium. This observation is further confirmed by the energy curves shown in Figure 11. In fact, for $\alpha=0.9$, the dynamical process reaches the steady state around $t=600$ and the corresponding energy is nearly zero; while for $\alpha=0.7$ and $0.4$, even at $t=700$, it has not reached the steady state.

3.5 3D bubble merging

We finally consider the 3D tFAC model with $\mathcal{M}=1$, $\varepsilon=0.03$, and $\alpha=0.5$. In this example, the double-well potential (2.2) is considered again. The adaptive $L$1-sESAV scheme is employed to simulate the bubble merging with an initial condition $\phi_{\text{init}}(x,y,z)=\max\{\phi_{1},\phi_{2}\}$ such that

[TABLE]

In this test, we use $80\times 80\times 80$ uniform meshes in space to discretize the computational domain $\Omega=(-0.5,0.5)^{3}$, and choose $\tau_{\max}=1$, $\tau_{\min}=0.01$, and $\eta=10^{7}$ in the adaptive time-stepping strategy (3.4). Figure 12 shows the iso-surfaces (value [math]) of the numerical solution obtained by the adaptive $L$1-sESAV scheme at different time instants, clearly illustrating that the two balls gradually shrink and merge into one smaller ball. Moreover, the discrete energy stability and MBP are well preserved at all times, as shown in Figures 13(a)–13(b). Finally, Figure 13(c) displays the history curves of the adaptive time steps, demonstrating the high efficiency of the time-adaptivity technique.

4 Stabilized $L2$-$1_{\sigma}$ type ESAV schemes

In this section, we extend the ideas and derivations presented in Section 2 to develop a second-order in time MBP-preserving scheme using the nonuniform $L2$-$1_{\sigma}$ formula JCP_Alikhanov_2015 for the Caputo derivative. For a given grid function $\{w^{k}\}$ with $k\geq 1$ and $0<\varsigma<1$, we define the off-set time level $t_{k-\varsigma}:=(1-\varsigma)t_{k}+\varsigma t_{k-1}$ and the weighted operator $w^{k-\varsigma}:=(1-\varsigma)w^{k}+\varsigma w^{k-1}$. Following JCP_Alikhanov_2015 , CiCP_Liao_2021 , we introduce the $L2$-$1_{\sigma}$ formula on a general nonuniform grid to approximate the Caputo derivative at $t=t_{n-\varsigma}$ with $\varsigma=\alpha/2$:

[TABLE]

where the discrete convolution kernels $\{B_{n-k}^{(n)}\}$ are defined as follows: $B^{(1)}_{0}=a_{0}^{(1)}$ if $n=1$, and

[TABLE]

for $n\geq 2$, with the discrete coefficients given by

[TABLE]

[TABLE]

It is easy to check that $a_{0}^{(n)}=\frac{1}{\tau_{n}}\int^{t_{n-\varsigma}}_{t_{n-1}}\omega_{1-\alpha}(t_{n-\varsigma}-s)ds=\frac{(1-\varsigma)^{1-\alpha}}{\Gamma(2-\alpha)\tau_{n}^{\alpha}}>0$. Besides, it is shown in CiCP_Liao_2021 , JCP_Liao_2020 that the discrete convolution kernels in (4.2) satisfy the following properties.

Lemma 4.1** (JCP_Liao_2020 ).**

Assume that the time-step ratio $r_{k}\geq 4/7$ for $2\leq k\leq n$ holds. Then

• (i)

the discrete kernels $B_{n-k}^{(n)}$ fulfill $B_{0}^{(n)}\leq\frac{24}{11\tau_{n}}\int_{t_{n-1}}^{t_{n}}\omega_{1-\alpha}\left(t_{n}-s\right)ds$ and

[TABLE]

• (ii)

the discrete kernels $B_{n-k}^{(n)}$ are monotone, i.e., $B_{n-k+1}^{(n)}<B_{n-k}^{(n)}$ for $1\leq k\leq n$;

• (iii)

and the first kernel $B_{0}^{(n)}$ is appropriately larger than the second one, i.e.,

[TABLE]

The following lemma is presented in [JSC_Liao_2024, , Lemmas 2.1–2.2 & Theorem 2.1], which is used to prove the discrete energy stability of the scheme proposed in next subsection.

Lemma 4.2** (JSC_Liao_2024 ).**

Let $r_{*}=r_{*}(\alpha)$ be the unique positive root of the nonlinear equation

[TABLE]

Then, it holds that $0.3865\approx r_{*}(0)<r_{*}(\alpha)<r_{*}(1)\approx 0.4037$ for $\alpha\in(0,1)$. Furthermore, assume that the adjacent time-step ratio $r_{k}\geq r_{*}(\alpha)$ for $2\leq k\leq n$, then it holds

[TABLE]

where the nonnegative functionals $\mathcal{G}$ is defined by

[TABLE]

with $C_{0}^{(1)}=4(1-\alpha)a_{0}^{(1)}/(2-\alpha)$, $C_{0}^{(n)}=4(1-\alpha)a_{0}^{(n)}/(2-\alpha)+2b_{1}^{(n)}/r_{n}$ for $n\geq 2$, and $C_{n-k}^{(n)}=B_{n-k}^{(n)}$ for $1\leq k\leq n-1$.

4.1 Stabilized $L2$-$1_{\sigma}$-ESAV scheme

Following the concept of the $L$1-sESAV scheme developed in Section 2, we apply the nonuniform $L2$-$1_{\sigma}$ formula (4.1) for the approximation of the Caputo derivative in (2.8) to derive the stabilized second-order ESAV scheme (denoted as $L2$-$1_{\sigma}$-sESAV):

Step 1

Let $\hat{\phi}^{0}=\phi^{0}$. For $n=1$, first solve a predicted solution $\hat{\phi}^{1}$ from the fully-implicit difference scheme:

[TABLE]

where $\hat{\phi}^{1-\varsigma}:=(1-\varsigma)\hat{\phi}^{1}+\varsigma\phi^{0}$, and let $V^{1-\varsigma}:=V(g_{h}(\hat{\phi}^{1-\varsigma},R^{0}))$, then find $\{\phi^{1},R^{1}\}\in\mathbb{V}_{h}\times\mathbb{R}$ by

[TABLE]

[TABLE]

Step 2

For $n\geq 2$, given $\phi^{n-1}$ and $\phi^{n-2}$, first predict a MBP-preserving solution $\hat{\phi}^{n}$ by (2.14) and let $\hat{\phi}^{n-\varsigma}:=(1-\varsigma)\hat{\phi}^{n}+\varsigma\phi^{n-1}$, $V^{n-\varsigma}:=V(g_{h}(\hat{\phi}^{n-\varsigma},R^{n-1}))$, then find $\{\phi^{n},R^{n}\}\in\mathbb{V}_{h}\times\mathbb{R}$ such that

[TABLE]

[TABLE]

Remark 7*.*

As mentioned in Remark 2, the nonlinear scheme (4.3) in Step 1 of the $L2$-$1_{\sigma}$-sESAV scheme can be solved by performing the following simple iterative scheme with a given termination error $tol$:

[TABLE]

where $\hat{\phi}^{1-\varsigma}_{(\ell)}=(1-\varsigma)\hat{\phi}^{1}_{(\ell)}+\varsigma\phi^{0}$ with $\ell=\mathfrak{s}-1,\mathfrak{s}$, and the initial iteration value is set as $\hat{\phi}^{1}_{(0)}=\phi^{0}$. Therefore, the practical implementation of the $L2$-$1_{\sigma}$-sESAV scheme can refer to that of the $L1$-sESAV scheme described in Algorithm 1. The following Theorem 4.3 demonstrates that the nonlinear scheme (4.3) is MBP-preserving and uniquely solvable. Following similar mathematical induction argument as in Theorem 2.6, we can prove Theorem 4.3, see C for details.

Theorem 4.3**.**

Assume that the condition in Lemma 4.1 holds and $\kappa\geq\left\|f^{\prime}\right\|_{C[-\beta,\beta]}$. If the first level time-step $\tau_{1}$ is sufficiently small such that

[TABLE]

the iterative scheme (4.8) preserves the MBP for $\{\hat{\phi}^{1}_{(\mathfrak{s})}\}$. Furthermore, if $\tau_{1}$ satisfies

[TABLE]

then the iterative scheme (4.8) converges to the unique solution $\hat{\phi}^{1}$ of the nonlinear scheme (4.3) in the maximum norm satisfying $\|\hat{\phi}^{1}\|_{\infty}\leq\beta$.

Remark 8*.*

The unique solvability of the $L2$-$1_{\sigma}$-sESAV scheme can be discussed similarly as Theorem 2.6, by equivalently rewriting (4.4)–(4.5) and (4.6)–(4.7) in a unified form:

[TABLE]

[TABLE]

Since $B^{(n)}_{0}>0$, $0<\varsigma<1$ and assumption (A2), it can be verified that the coefficient matrix of (4.11) is symmetric and positive definite. In fact, it is also strictly diagonally dominant.

4.2 Discrete energy stability and MBP

The following two theorems establish the discrete energy stability and MBP for the proposed $L2$-$1_{\sigma}$-sESAV scheme (4.3)–(4.7).

Theorem 4.4**.**

Assume that the conditions in Lemma 4.2 and Theorem 4.3 hold, the $L2$-$1_{\sigma}$-sESAV scheme (4.3)–(4.7) is unconditionally energy-stable in the sense that $\mathcal{E}_{h}[\phi^{n},R^{n}]\leq\mathcal{E}_{h}[\phi^{0},R^{0}]$ for $n\geq 1$.

Proof.

Let $n=k$ in (4.6) and take the inner product with $\nabla_{\tau}\phi^{k}=\phi^{k}-\phi^{k-1}$. We then combine the resulting equation with (4.12) to derive

[TABLE]

It is straightforward to verify that

[TABLE]

Thus, by taking $\sigma=\varsigma=\alpha/2$ with $\alpha\in(0,1)$ in (4.14), the third left-hand side term of (4.13) can be estimated below by

[TABLE]

Inserting (4.15) into (4.13), we get

[TABLE]

Similarly, it follows from (4.4) and (4.12) that (4.16) also holds for $k=1$.

Now, we sum the above inequality (4.16) from $k=1$ to $n$ to derive

[TABLE]

where Lemma 4.2 has been used in the first inequality. The proof is completed. ∎

Theorem 4.5**.**

Assume that the conditions in Theorem 4.3 hold. Furthermore, if the time-step satisfies

[TABLE]

where $K_{2}$ is the upper bound of $V(\cdot)$, the $L2$-$1_{\sigma}$-sESAV scheme (4.3)–(4.7) preserves the MBP for $\left\{\phi^{n}\right\}$, that is,

[TABLE]

Proof.

For $n=1$, it follows from (4.11) that

[TABLE]

with $\mathbf{M}_{2}^{1}:=(B^{(1)}_{0}-\kappa\varsigma\mathcal{M}V^{1-\varsigma})I+\varsigma\mathcal{M}\varepsilon^{2}\Delta_{h}$. Theorem 4.3 demonstrates that $\|\hat{\phi}^{1}\|_{\infty}\leq\beta$ when $\tau_{1}$ satisfies (4.9)–(4.10), thereby ensuring that $\|\hat{\phi}^{1-\varsigma}\|_{\infty}\leq\beta$. Then, the application of Lemmas 2.4–2.5 leads to

[TABLE]

Moreover, if (4.17) holds for $n=1$, all elements of $\mathbf{M}_{2}^{1}$ are nonnegative, and thus

[TABLE]

which yields $\|\mathbf{M}_{2}^{1}\phi^{0}\|_{\infty}\leq\|\mathbf{M}_{2}^{1}\|_{\infty}\|\phi^{0}\|_{\infty}\leq B_{0}^{(1)}\beta-\kappa\varsigma\mathcal{M}V^{1-\varsigma}\beta$. Inserting this estimate into (4.19), we obtain $\|\phi^{1}\|_{\infty}\leq\beta$.

For $2\leq n\leq N$, suppose that $\|\phi^{k}\|_{\infty}\leq\beta$ for $0\leq k\leq n-1$. It follows from (2.14) and the definition of $\hat{\phi}^{n-\varsigma}$ that $\|\hat{\phi}^{n-\varsigma}\|_{\infty}\leq\beta$. Thus, it remains to verify that $\|\phi^{n}\|_{\infty}\leq\beta$. Noting that (4.1) can be expressed as

[TABLE]

Then, equation (4.6) can be rewritten as

[TABLE]

with $\mathbf{M}_{2}^{n}:=(B_{0}^{(n)}-B_{1}^{(n)}-\kappa\varsigma\mathcal{M}V^{n-\varsigma})I+\varsigma\mathcal{M}\varepsilon^{2}\Delta_{h}.$

Under the time-step condition (4.17) and Lemma 4.1 (i) and (iii), we see

[TABLE]

which further means that all elements of $\mathbf{M}_{2}^{n}$ are nonnegative and

[TABLE]

Consequently, the hypothesis $\|\phi^{n-1}\|_{\infty}\leq\beta$ yields

[TABLE]

Moreover, the positivity and monotonicity of the discrete kernels $\{B^{(n)}_{n-k}\}$, along with the assumption that $\|\phi^{k}\|_{\infty}\leq\beta$ for $0\leq k\leq n-1$, lead to

[TABLE]

Therefore, collecting the estimates (4.21)–(4.22) and using Lemmas 2.4–2.5, it follows from (4.20) that

[TABLE]

This immediately implies $\|\phi^{n}\|_{\infty}\leq\beta$ and completes the proof. ∎

4.3 Novel unbalanced stabilized $L2$-$1_{\sigma}$-sESAV scheme

In cases that the auxiliary functional $V(\cdot)$ satisfies the assumptions (A1)–(A3) with $K_{2}=1$, we can further improve the $L2$-$1_{\sigma}$-sESAV scheme by introducing a novel unbalanced stabilization term

[TABLE]

The imbalance, reflected by $V^{n-\varsigma}\neq 1$ for $n\geq 1$ generally, is beneficial for the preservation of MBP especially for long-term simulations with a large time step, see Remark 9 and Figures 14–15 in Example 4.4.2. Using the modified stabilization term (4.24), a novel unbalanced $L2$-$1_{\sigma}$-sESAV scheme is proposed as follows:

Step 1

Let $\hat{\phi}^{0}=\phi^{0}$. For $n=1$, first solve a predicted solution $\hat{\phi}^{1}$ from the fully-implicit difference scheme (4.3), then find $\{\phi^{1},R^{1}\}\in\mathbb{V}_{h}\times\mathbb{R}$ by

[TABLE]

[TABLE]

Step 2

For $n\geq 2$, given $\phi^{n-1}$ and $\phi^{n-2}$, first predict a MBP-preserving solution $\hat{\phi}^{n}$ by (2.14), then find $\{\phi^{n},R^{n}\}\in\mathbb{V}_{h}\times\mathbb{R}$ such that

[TABLE]

[TABLE]

Following the proofs of Theorems 4.4–4.5, together with Lemmas 4.1–4.3, the preservation of discrete energy stability and MBP for the unbalanced $L2$-$1_{\sigma}$-sESAV scheme can be established in a similar manner.

Theorem 4.6**.**

Assume that the conditions in Lemma 4.2 and Theorem 4.3 hold, the unbalanced $L2$-$1_{\sigma}$-sESAV scheme (4.25)–(4.28) is unconditionally energy-stable in the sense that $\mathcal{E}_{h}[\phi^{n},R^{n}]\leq\mathcal{E}_{h}[\phi^{0},R^{0}]$ for $n\geq 1$.

Theorem 4.7**.**

Assume that the conditions in Theorem 4.3 hold. If the time-step satisfies

[TABLE]

the unbalanced $L2$-$1_{\sigma}$-sESAV scheme (4.25)–(4.28) preserves the MBP for $\left\{\phi^{n}\right\}$, that is,

[TABLE]

Proof.

Together with a similar treatment to (4.18)–(4.23), we conclude from (4.25) and (4.27) that

[TABLE]

under the ime-step condition (4.29), where $\mathcal{P}^{n}:=\kappa\mathcal{M}(V^{n-\varsigma}-1)\beta$. The assumption (A2) with $K_{2}=1$ ensures that $\mathcal{P}^{n}\leq 0$, thereby further guarantees that the unbalanced $L2$-$1_{\sigma}$-sESAV scheme preserves the discrete MBP. ∎

Remark 9*.*

The result (4.30) directly leads to

[TABLE]

From assumption (A3), we know that $g_{h}(\hat{\phi}^{n-\varsigma},R^{n-1})$ and $V^{n-\varsigma}$ may deviate significantly from $1$ for a relatively larger time step. In this sense, the negative term ${\mathcal{P}^{n}}/{(B^{(n)}_{0}+\kappa(1-\varsigma)\mathcal{M})}$ can be considered as a penalty term to generate solutions with smaller maximum-norm. However, no such term exists in the standard $L2$-$1_{\sigma}$-sESAV scheme, indicating that, compared to the standard one, the proposed unbalanced $L2$-$1_{\sigma}$-sESAV scheme is more effective in preserving the discrete MBP, see also the numerical results in Figures 14–15 for comparisons.

4.4 Numerical experiments

This subsection is devoted to numerical tests of both the $L2$-$1_{\sigma}$-sESAV and unbalanced $L2$-$1_{\sigma}$-sESAV schemes. Note that for the $L2$-$1_{\sigma}$ type methods, a time-step ratio constraint $r_{k}\geq 4/7$ is required to preserve the discrete energy stability (cf. Theorems 4.4 and 4.6) and MBP (cf. Theorems 4.5 and 4.7). Thus, in the numerical simulations, the following modified adaptive time-stepping strategy

[TABLE]

is adopted, where $r_{\min}\geq 4/7$ is a user-defined constant.

4.4.1 Temporal convergence

First, we verify the temporal accuracy of both the balanced and unbalanced $L2$-$1_{\sigma}$-sESAV schemes, and the numerical results are obtained by solving the exterior-forced problem (3.5) with a fixed fractional order $\alpha=0.8$ on $200\times 200$ uniform spatial grids. The exact solution is prescribed as $\phi(\mathbf{x},t)=0.5\,\omega_{1+\iota}(t)\sin x\sin y$ through an appropriate choice of the linear term $g(\mathbf{x},t)$. For different regularity parameters $\iota=0.3,0.5$ and $0.8$, the mesh grading parameter is set to be $\gamma=2/\iota$ such that the theoretical convergence order is $\min\{2,\gamma\iota\}=2$. The other setup is the same as in Example 3.1. Numerical results for different $\iota$ and potentials are summarized in Tables 4–5, from which we can clearly observe that (i) both methods exhibit second-order temporal accuracy as expected; and (ii) they generate numerical solutions with the same magnitude accuracy.

4.4.2 Effectiveness of the unbalanced stabilization term

The coarsening dynamics of the tFAC model (1.1) is examined with $\alpha=0.9$, $\mathcal{M}=1$, and $\varepsilon=0.01$. We adopt a $128\times 128$ uniform spatial mesh to discretize the spatial domain $\Omega=(0,1)^{2}$, and give the initial value by random numbers uniformly distributed between $-0.8$ and $0.8$ at each grid point, see Figure 1.

For special choice $V(\cdot)$ in (3.1), it deduces from Theorems 4.5 and 4.7 that $\tau\leq 0.1044$ for the double-well potential and $\tau\leq 0.0646$ for the Flory–Huggins potential are required to preserve the discrete MBP on uniform temporal grids. Note that these conditions are only sufficient but not necessary. Therefore, in the following tests, three different types of temporal girds are employed for each potential, two of which violate the time-step restriction while the other satisfies it. Specifically, we take $\tau=1,0.5,0.1$ for the double-well potential and $\tau=1/3,1/6,1/30$ for the Flory–Huggins potential. The coarsening dynamics are tested using the two $L2$-$1_{\sigma}$-sESAV schemes with different stabilization terms. Numerical results for the double-well and Flory–Huggins potentials are depicted in Figures 14 and 15, respectively. From these results, we observe that: (i) when the time steps satisfy the restriction (4.29) (or (4.17) with $K_{2}=1$), both schemes preserve the discrete MBP, see Figures 14(c) and 15(c); (ii) when the time steps moderately violate the constraint, both schemes can still preserve the discrete MBP, indicating that the conditions in Theorems 4.5 and 4.7 are only sufficient but not necessary, see Figures 14(b) and 15(b); (iii) however, when the given largest time steps are used for both schemes, the unbalanced $L2$-$1_{\sigma}$-sESAV scheme can still preserve the discrete MBP, whereas the standard $L2$-$1_{\sigma}$-sESAV scheme yields solutions exceeding the corresponding maximum-bound $\beta$, see Figures 14(a) and 15(a). This demonstrates the effectiveness and excellence of the proposed unbalanced stabilization term in preserving the discrete MBP.

4.4.3 Long-term coarsening dynamics simulations

The unbalanced $L2$-$1_{\sigma}$-sESAV scheme has demonstrated superior reliability in preserving the MBP compared to the standard scheme, making it the preferred choice for simulating long-term coarsening dynamics governed by the tFAC model (1.1), particularly for extended time scales up to $T=2000$ in the following test. Moreover, to improve the computational efficiency without sacrificing accuracy in the simulation, a mixed adaptive time-stepping strategy (3.2) and (4.31) with $\hat{T}=0.5$ and $\hat{N}=30$ is employed. Here we set $\mathcal{M}=1$, $\varepsilon=0.01$, $r_{\min}=4/7$, $\eta=10^{7}$, and $(\tau_{\min},\tau_{\max})=(0.1,1)$ for the double-well potential and $(\tau_{\min},\tau_{\max})=(1/30,1/3)$ for the Flory–Huggins potential. The random initial values and spatial partition are the same as aforementioned.

For the double-well potential, the profiles of coarsening dynamics with different fractional orders $\alpha=0.9,0.7$ and $0.4$ are depicted in Figure 16, where snapshots are taken at times $t=20,100,500$, and $2000$, respectively. It is clearly observed that the coarsening speed governed by the tFAC equation depends heavily on the fractional order $\alpha$ and smaller values of $\alpha$ require much longer time to reach the steady state. Moreover, Figure 17(a)–17(b) demonstrate that the proposed unbalanced $L2$-$1_{\sigma}$-sESAV scheme with the adaptive time-stepping strategy successfully preserves both the discrete energy stability and MBP, and Figure 17(c) further shows that the time steps are accurately and effectively selected to capture the energy variations in the adaptive method. For the Flory–Huggins potential, the corresponding numerical results are presented in Figures 18–19, and similar conclusions can also be drawn.

5 Concluding remarks

In this paper, we first established the MBP for the tFAC model (1.1) with two different nonlinear potentials. Subsequently, we developed and analyzed a linear stabilized nonuniform time-stepping scheme, referred to as the $L1$-sESAV scheme, in which an essential auxiliary functional was introduced to ensure that the first-order approximation of the SAV does not compromise the temporal accuracy of the phase function. This treatment is critical for deriving unconditional energy stability of nonuniform time-stepping linear scheme. Moreover, the non-negativity of this novel auxiliary functional also allows for the design of a stabilization term that effectively controls the nonlinear term, which together with a newly developed prediction strategy guarantees the discrete MBP-preservation. To the best of our knowledge, this is the first $(2-\alpha)$th order $L1$-type linear scheme with unconditional preservation of the energy stability and MBP. Building on these ideas, we further established a linear second-order scheme, termed the $L2$-$1_{\sigma}$-sESAV scheme, which can preserve the energy stability unconditionally and MBP conditionally. To enhance the discrete MBP preservation for large time steps, we also introduced an improved unbalanced stabilization term by imposing appropriate boundedness and monotonicity assumptions on the auxiliary functional, and then an unbalanced $L2$-$1_{\sigma}$-sESAV scheme was developed and discussed. Numerical simulations for the 2D and 3D tFAC models with either the double-well potential or the Flory–Huggins potential suggest the satisfactory and high effectiveness of the proposed methods. Moreover, the ideas and derivations presented in this paper enable the development of MBP-preserving numerical methods based on other discrete Caputo derivatives, such as the $L1_{R}$ formula SISC_Liao_2021 and the fast versions of $L1$ ACM_Ji_2020 and $L2$-$1_{\sigma}$ JCP_Liao_2020 formulas. Currently, rigorous error analyses of SAV-type methods for time-fractional phase-field models under weak singularity assumptions remain scarce, representing an important direction for future research.

CRediT authorship contribution statement

Bingyin Zhang: Methodology, Formal analysis, Software, Writing-Original draft, Funding acquisition. Hongfei Fu: Methodology, Conceptualization, Supervision, Writing-Reviewing and Editing, Funding acquisition.

Data availability

A free github repository contains the Matlab code employed, which can be accessed through the following link: https://github.com/HFu20/time-Fractional-Allen-Cahn/tree/sESAV-MBP.

Declaration of competing interest

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

Acknowledgements

B. Zhang was supported in part by the Fundamental Research Funds for the Central Universities (No. 202461103). H. Fu was supported in part by the National Natural Science Foundation of China (Nos. 11971482, 12131014), by the Shandong Provincial Natural Science Foundation (No. ZR2024MA023), by the Fundamental Research Funds for the Central Universities (No. 202264006), by the OUC Scientific Research Program for Young Talented Professionals.

Appendix A Proof of Lemma 2.1

For convenience, we set $w(t)=v(t)-v(t_{*})$. Then, we have $\partial_{t}w(t)=\partial_{t}v(t)$ and $w(t_{*})=0$, which further gives

[TABLE]

Moreover, since $v\in C[0,T]\cap C^{1}(0,T]$, one can easily verify that $\lim\limits_{s\rightarrow t_{*}}(t_{*}-s)^{-\alpha}w(s)=0$. Thus, applying integration by parts yields

[TABLE]

If $v$ attains its minimum at $t_{*}\in(0,T]$, then $w(s)\geq 0$ for all $s\in[0,T]$. Consequently, we obtain ${}_{0}^{C}D^{\alpha}_{t}v(t_{*})\leq 0.$ An analogous argument yields the corresponding inequality when $v$ attains its maximum.

Appendix B Proof of Lemma 2.3

The integral version of Taylor’s theorem gives

[TABLE]

Under the regularity assumption, one has

[TABLE]

Specifically, when $n=2$, the second right-hand side term of (B.1) can be estimated by

[TABLE]

Inserting the above estimates into (B.1) completes the proof of (2.18).

In particular, when the graded temporal grids are used, we recall the conclusion in SINUM_Mustapha_2014 that there exists a constant $C_{\gamma}>0$ such that $\tau_{k}\leq C_{\gamma}t_{k}^{1-1/\gamma}N^{-1}$ and $t_{k}\leq C_{\gamma}t_{k-1}$ for $2\leq k\leq N$, which further implies $\tau_{1}+\tau_{2}\leq 2\tau_{2}\leq C_{\gamma}^{\prime}N^{-\gamma}$ and thus $(\tau_{1}+\tau_{2})^{\iota}\leq C_{\gamma}^{\prime}N^{-\gamma\iota}$. Denote $\zeta=\min\{2,\gamma\iota\}$, then we have $\tau_{\ell}^{\zeta}\leq C_{\gamma}t_{\ell}^{\zeta-\zeta/\gamma}N^{-\zeta}$ with $\ell=n,n-1$. Thus, it holds that

[TABLE]

Finally, by inserting the above estimates into (2.18), we complete the proof of (2.19).

Appendix C Proof of Theorem 4.3

For any $\mathfrak{s}\geq 1$, we assume that $\|\hat{\phi}^{1}_{(\mathfrak{s}-1)}\|_{\infty}\leq\beta$, which further implies $\|\hat{\phi}^{1-\varsigma}_{(\mathfrak{s}-1)}\|_{\infty}\leq\beta$. Equivalently, equation (4.8) can be rewritten as

[TABLE]

where $\mathbf{M}_{1}:=(B^{(1)}_{0}-\kappa\varsigma\mathcal{M})I+\varsigma\mathcal{M}\varepsilon^{2}\Delta_{h}$. Due to $B_{0}^{(1)}\geq\frac{4}{11}\int_{t_{0}}^{t_{1}}\frac{\omega_{1-\alpha}\left(t_{1}-s\right)}{\tau_{1}}\mathrm{~d}s=\frac{4\tau_{1}^{-\alpha}}{11\Gamma(2-\alpha)}$ (see Lemma 4.1 (i)), it then follows from (4.9) that all elements of $\mathbf{M}_{1}$ are nonnegative, and thus, we have

[TABLE]

which together with Lemmas 2.4–2.5 gives us

[TABLE]

Consequently, $\|\hat{\phi}^{1}_{(\mathfrak{s})}\|_{\infty}\leq\beta$ holds for any $\mathfrak{s}\geq 1$.

Furthermore, equation (4.8) is also equivalent to

[TABLE]

which defines a mapping $\mathcal{T}[v]=w$ from $\mathbb{V}_{\beta}$ to $\mathbb{V}_{\beta}$, i.e.,

[TABLE]

Similar as in Theorem 2.6, it is straightforward to verify that $\mathcal{T}$ is a contractive mapping under the time-step condition (4.10). Consequently, the iterative scheme (4.8) admits a unique fixed point in $\mathbb{V}_{\beta}$, denoted by $\hat{\phi}^{1-\varsigma}$, and $\hat{\phi}^{1}=\frac{\phi^{1-\varsigma}-\varsigma\phi^{0}}{1-\varsigma}$ is just the unique solution of the nonlinear scheme (4.3). The proof of Theorem 4.3 is completed.