<figure><img src="image_1.png" /></figure>

HAL
open science

# A bound-preserving scheme for the Allen-Cahn equation

Zhaonan Dong, Alexandre Ern, Zuodong Wang

► To cite this version:

Zhaonan Dong, Alexandre Ern, Zuodong Wang. A bound-preserving scheme for the Allen-Cahn equation. Computers & Mathematics with Applications, In press. hal-04931694v2

HAL Id: hal-04931694

https://hal.science/hal-04931694v2

Submitted on 26 Jun 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# A bound-preserving scheme for the Allen–Cahn equation

Zhaonan Dong†*, Alexandre Ern*†, and Zuodong Wang*†

June 26, 2025

## Abstract

We propose and analyze a bound-preserving scheme for the Allen–Cahn equation. The key idea is to apply a bound-preserving nonlinear stabilization technique to the implicit Euler time-stepping method coupled with the continuous finite element method. To our best knowledge, this is the first scheme which theoretically preserves the maximum principle and has an error estimate that is optimal in the $L^2(H^1)$-seminorm and with a polynomial dependence on $\epsilon^{-1}$ at the same time. The proof of the error estimate combines a nonlinear Ritz projection together with a special Grönwall inequality. Numerical experiments are conducted to compare the performance of our scheme with a bound-preserving operator-splitting scheme.

**Keywords.** bound-preserving scheme, error estimate, Allen–Cahn equation, stiff source, phase field model.

**MSC.** 65M12, 65M60

## 1 Introduction

We consider the Allen–Cahn equation introduced in [2]:

$$ egin{cases} \partial_t u - \Delta u + \epsilon^{-2} f(u) = 0, & \forall (x, t) \in D \times J, \\ u(x, 0) = u_0(x), & \forall x \in D, \\ \partial_n u(x, t) = 0, & \forall (x, t) \in \partial D \times J, \end{cases} \quad (1) $$

where $D$ denotes a bounded domain in $\mathbb{R}^d$, $d \in \{1, 2, 3\}$, with a Lipschitz boundary, and $J := (0, T)$ is the time interval with a given final time $T > 0$. We assume the initial condition $u_0 \in H^1(D; [-1, 1])$. The function $f(v) := v(v^2 - 1)$ represents the forcing term, and its anti-derivative, $F(v) := \frac{1}{4}(v^2 - 1)^2$, is referred to as the double-well potential. $\epsilon \in \mathbb{R}^+ := (0, +\infty)$ represents the scaling of the reaction term. The Allen–Cahn equation constitutes one of the simplest phase-field models. We refer to [19] and the references therein for a detailed description of general phase-field models.

The principal challenge in the numerical approximation of the Allen–Cahn equation lies in resolving the reaction source term when the parameter $\epsilon$ is very small compared to the characteristic time and viscous scales of the equation. To fix the ideas, we assume that $T = O(1)$, the diameter of $D$ is $O(1)$ and $\epsilon \ll 1$, i.e., we focus on the regime where the time horizon is large enough and reaction dominates diffusion. This regime is commonly referred to as the stiff regime. Without loss of generality, we assume that $\epsilon \le 1$ in what follows.

<footnote>*CERMICS, ENPC, Institut Polytechnique de Paris, 77455 Marne-la-Vallee Cedex 2, France</footnote>

<footnote>†SERENA, Centre Inria de Paris, 48 rue Barrault, 75647, France</footnote>
It is well-known that two main numerical challenges of the Allen–Cahn equation are the bound-preserving property (i.e., the numerical solution should remain within the interval $[-1, 1]$) and the exponential $\epsilon$-dependence in the error analysis. The first challenge arises from the requirement of a physically meaningful solution, while the second stems from bounding the nonlinear source term $\epsilon^{-2}f(u)$. The bound-preserving property has been well studied through the application of various techniques, such as the cut-off method used in Yang et al. [28], Li et al. [22], the minimization-based post-processing proposed in Liu et al. [23], the monotonicity theory applied in Shen and Zhang [25], and the contractive semigroup theory considered in Du et al. [14].

However, the issue of exponential $\epsilon$-dependence of the error bound is not explicitly addressed in the above works. This second challenge arises from the use of the Grönwall inequality, which overestimates the contribution of the source term. To the best of the authors' knowledge, the only existing approach to mitigate this issue involves using the principal eigenvalue of a linearized eigenvalue problem. For details on the derivation and basic properties of this problem, we refer to Chen [11]. Subsequent numerical studies have observed that the bound of this eigenvalue is related to topological changes in the zero-level set of the solution to the Allen–Cahn equation, see, for example, Bartels and Müller [8], Bartels et al. [9]. Here, a topological change refers to a situation where the number of connected components of the zero-level set changes, see Bartels [7, Section 6.1.4]. To leverage this eigenvalue and reduce the $\epsilon$-dependence, various techniques have been proposed: (1) A generalized Grönwall inequality introduced in Feng and Prohl [18, 19], Bartels et al. [9], Bartels [7], Feng and Li [17], (2) A well-designed dual problem avoiding the Grönwall-type argument in Chrysafinos [12], and (3) A truncation technique for the stiff source term in Akrivis and Li [1]. However, the bound-preserving property is not explicitly considered in these works.

The main novelty of the present work is the development of a scheme that enjoys at the same time a bound-preserving property and an error bound with only polynomial $\epsilon$-dependence. In our scheme, the nonlinear stabilization technique from Barrenechea et al. [4], Amiri et al. [3], Barrenechea et al. [5] is employed to ensure the bound-preserving property, while the generalized Grönwall inequality proposed in Bartels [6, Lemma 6.2] is utilized to mitigate the $\epsilon$-dependence. The choice of the above nonlinear stabilization technique is motivated by the following reasons: (1) Effectiveness in damping oscillations: Unlike the cut-off technique, this technique suppresses small oscillations. (2) Error analysis compatibility: This technique enables optimal error analysis in the energy norm, which remains a challenge for the minimization-based post-processing approach.

The rest of the paper is organized as follows. The model problem is detailed in §2. The discrete setting, the proposed scheme together with its fundamental properties are discussed in §3. In §4, the error estimate is derived, the main results being Theorem 4.4 and Corollary 4.5. Finally, §5 provides implementation details of our scheme, presents numerical tests, and compares our scheme with another simple bound-preserving scheme from the literature.

## 2 Model problem

We denote the dual space of $H^1(D)$ by $H^1(D)'$. The $L^2$-inner product is abbreviated as $\langle \cdot, \cdot \rangle$, the associated norm as $\|\cdot\|$, and the duality pairing between $H^1(D)'$ and $H^1(D)$ as $\langle \cdot, \cdot \rangle$. We omit the dependence in time and space for inner products, duality pairs and norms when there is no ambiguity, e.g., $\|\cdot\|_{L^2(L^2)} := \|\cdot\|_{L^2(J;L^2(D))}$. In various bounds, we denote by $C$ any generic constant (its value can change at each occurrence) that is independent of the mesh size, the time step and the reaction parameter $\epsilon$.

The model problem (1) can be posed in the following weak sense: Find $u \in L^2(J; H^1(D)) \cap$
75 $H^1(J; H^1(D)')$ such that

$$ \begin{cases} \int_J [\langle \partial_t u, v \rangle + (\nabla u, \nabla v) + \epsilon^{-2}(f(u), v)] = 0, & \forall v \in L^2(J; H^1(D)), \\ u(\cdot, 0) = u_0 \in H^1(D; [-1, 1]). \end{cases} \tag{2} $$

76 Notice that $(f(u), v)$ is meaningful owing to the Sobolev embedding $H^1(D) \hookrightarrow L^6(D)$.

77 We recall the following result on the solution to (1): The well-posedness from [18, Proposition 1], the maximum principle and the regularity result from [18, Proposition 2].

79 **Lemma 2.1** (Well-posedness and maximum principle). There exists a unique solution $u$ to the model problem (1). This solution belongs to $L^\infty(J; H^1(D)) \cap H^1(J; L^2(D))$ and satisfies the following maximum principle: For almost all $(x, t) \in D \times J$, $u(x, t) \in [-1, 1]$. In addition, we have the following energy estimates:

$$ \|u\|_{L^\infty(L^2)}^2 + 2\|\nabla u\|_{L^2(L^2)}^2 \le \|u_0\|^2 + C\epsilon^{-2}, \tag{3a} $$

$$ \sup_{t \in J} \mathcal{J}_\epsilon(u(\cdot, t)) + \int_J \|\partial_t u\|^2 \le \mathcal{J}_\epsilon(u_0), \tag{3b} $$

83 where the total energy $\mathcal{J}_\epsilon: H^1(D) \to \mathbb{R}^+$ is defined as

$$ \mathcal{J}_\epsilon(v) := \frac{1}{2} \|\nabla v\|^2 + \frac{1}{\epsilon^2} \int_D F(v), \quad \forall v \in H^1(D). \tag{3c} $$

84 **Lemma 2.2** (Higher regularity). Assume that the Laplace operator with homogeneous Neumann boundary condition is $H^2$-regular (i.e., $\|\nabla^2 u\| \le C \|\Delta u\|$, see [21, Theorem 3.2.1.3].), $u_0 \in H^2(D)$ and $\lim_{t \to 0^+} \|\partial_t \nabla u(t)\| < \infty$. Then $u \in W^{1,\infty}(J; H^1(D)) \cap H^2(J; L^2(D)) \cap H^1(J; H^2(D))$.

87 As mentioned in the introduction, the notion of principal eigenvalue is the key ingredient for removing the exponential $\epsilon$-dependence in the error analysis. We now give its definition which was initially proposed in [11].

90 **Definition 2.3** (Principal eigenvalue). Let $u$ be the solution to the model problem (1). The principal eigenvalue $\lambda: J \to \mathbb{R}$ is defined for almost all $t \in J$ as

$$ \lambda(t) := \max \left\{ 0, - \inf_{v \in H^1(D) \setminus \{0\}} \frac{(\epsilon^{-2} f'(u(\cdot, t)), v^2) + \|\nabla v\|^2}{\|v\|^2} \right\}. \tag{4} $$

92 Based on numerical [9, 8] and theoretical [6] studies, we make the following assumption which will be used only in §4. The other results presented in this work do not require this assumption.

94 **Assumption 2.4** (Temporal integrability of principal eigenvalue). We assume that the principal eigenvalue defined in (4) satisfies $\lambda \in L^1(J)$ and

$$ \int_J \lambda(t) dt \le C_\lambda + \log(\epsilon^{-\kappa}), \tag{5} $$

96 where the constants $\kappa$ and $C_\lambda > 0$ are independent of $\epsilon$.

97 **Remark 2.5** (Topological changes). The constant $\kappa$ represents the total number of topological changes of the zero-level set of $u$ during the evolution in the time interval $J$.

99 **3 Discrete scheme**

100 In this section, we present the discrete setting, our numerical scheme, and some of its fundamental properties.
## 3.1 Time and space discrete setting

The time discretization is defined by using the discrete time nodes $t^n$ for all $n \in \{0:N\}$, with $t^0 = 0$ and $t^N = T$. The time step $\tau^n$ is defined as $\tau^n := t^n - t^{n-1}$ for $n \in \mathcal{N} := \{1:N\}$. Let $J^n := (t^{n-1}, t^n]$ be the local time interval for all $n \in \mathcal{N}$. To simplify the notation, we omit the superscript $n$ and denote the time step by $\tau$.

For the space discretization, we consider either a simplicial or a tensor-product conforming mesh $\{\mathcal{T}_h\}_{h>0}$ of $D$, with mesh size $h$. A generic mesh cell is denoted by $K \in \mathcal{T}_h$. It is assumed to be an open set, and its diameter is denoted by $h_K$. We also define $h_{\max} := \max_{K \in \mathcal{T}_h} h_K$ and $h_{\min} := \min_{K \in \mathcal{T}_h} h_K$. We only assume that the mesh is shape-regular, instead of the acute angle assumption which guarantees the discrete maximum principle in the literature (see for example [10]). We focus for simplicity on continuous, piecewise affine finite elements. Thus, each degree of freedom (dof) $i$ is the value at the mesh vertex $\mathbf{x}_i$ for all $i \in \mathcal{V} := \{1: I\}$, where $I$ is the number of mesh vertices. The global shape functions are denoted by $\{\varphi_i\}_{i \in \mathcal{V}}$. The corresponding $H^1$-conforming FEM space is denoted by

$$V_h := \operatorname{span}\{\varphi_i\}_{i \in \mathcal{V}}. \qquad (6)$$

The basis functions take non-negative values and satisfy the following partition of unity property:

$$\varphi_i(\mathbf{x}) \ge 0, \quad \forall \mathbf{x} \in \bar{D}, \quad \forall i \in \mathcal{V}, \qquad \sum_{i \in \mathcal{V}} \varphi_i(\mathbf{x}) = 1, \quad \forall \mathbf{x} \in \bar{D}. \qquad (7)$$

Moreover, for each dof $i \in \mathcal{V}$, we define its ‘mesh diameter’ $h_i$ as the diameter of $\omega(i) := \{\mathbf{x} \in D \mid \varphi_i(\mathbf{x}) > 0\}$.

We define $V_h^+ \subset V_h$ as the subset composed of bound-preserving trial functions in $V_h$, i.e.,

$$V_h^+ := \{v_h := \sum_{i \in \mathcal{V}} V_i \varphi_i \in V_h \mid V_i \in [-1, 1], \forall i \in \mathcal{V}\}. \qquad (8)$$

Notice that, owing to (7), functions in $V_h^+$ take values in $[-1, 1]$. For all $v_h := \sum_{i \in \mathcal{V}} V_i \varphi_i \in V_h$, we define the map $P_h^+ : V_h \to V_h^+$ by

$$v_h^+ := P_h^+(v_h) := \sum_{i \in \mathcal{V}} \mathcal{P}_{[-1,1]}(V_i)\varphi_i, \qquad (9)$$

where $\mathcal{P}_{[-1,1]} : \mathbb{R} \to [-1, 1]$ is the $l^2$-orthogonal projection onto the convex set $[-1, 1]$, i.e.,

$$\mathcal{P}_{[-1,1]}(a) := \max(-1, \min(1, a)) = \underset{b \in [-1,1]}{\operatorname{argmin}} |a-b|^2, \qquad (10)$$

for all $a \in \mathbb{R}$. We also define the residual part of $v_h$ as

$$v_h^- := v_h - v_h^+ \in V_h. \qquad (11)$$

## 3.2 Numerical scheme

The main idea of our numerical scheme is to find the discrete solution $u_h \in V_h$ from a nonlinear scheme in which the time evolution is realized on $u_h^+ \in V_h^+$, and the unisolvence of the nonlinear scheme is guaranteed by adding a stabilization term acting on $u_h^-$. More precisely, we apply the implicit Euler method in time and the $H^1$-conforming FEM discretization in space described in §3.1. At each time step $n \in \mathcal{N}$, we find $u_h^n \in V_h$ solving

$$\frac{1}{\tau}(u_h^{n,+} - u_h^{n-1,+}, v_h) + (\nabla u_h^{n,+}, \nabla v_h) + \frac{1}{\epsilon^2}(f(u_h^{n,+}), v_h) + s_h(u_h^{n,-}, v_h) = 0, \qquad (12)$$
for all $v_h \in V_h$, where the stabilization term $s_h$ is inspired from [4] and reads as follows:

$$s_h(v_h, w_h) := \alpha \sum_{i \in \mathcal{V}} \left( \frac{h_i^d}{\tau} + h_i^{d-2} + \frac{h_i^d}{\epsilon^2} \right) V_i W_i, \quad (13)$$

for all $v_h := \sum_{i \in \mathcal{V}} V_i \varphi_i$ and $w_h := \sum_{i \in \mathcal{V}} W_i \varphi_i \in V_h$. The positive parameter $\alpha = O(1)$ will be precised in Lemma 3.3. The initial data $u_h^0$ can be evaluated by using any nodally bound-preserving operator with optimal convergence rate in the $H^1$-norm, like the Lagrange interpolation for smooth initial data (i.e., for $u_0 \in C^0(\bar{D})$) or the nonlinear Ritz map defined in (26) for non-smooth initial data (i.e., for $u_0 \notin C^0(\bar{D})$).

## 3.3 Fundamental properties

This section collects some fundamental properties for our scheme. Some proofs, in particular those of Lemma 3.1 and 3.2 which are simple adaptations from arguments in [4], are outlined for completeness.

The first lemma is devoted to the basic properties of the stabilization term.

**Lemma 3.1 (Sign of stabilization).** The bilinear form $s_h(\cdot, \cdot)$ satisfies the following inequalities:

$$s_h(v_h^{-} - w_h^{-}, v_h^{+} - w_h^{+}) \geq 0, \quad \forall v_h, w_h \in V_h, \quad (14a)$$

$$s_h(v_h^{-}, w_h - v_h^{+}) \le 0, \quad \forall v_h \in V_h, \forall w_h \in V_h^{+}. \quad (14b)$$

*Proof.* We consider arbitrary $v_h := \sum_{i \in \mathcal{V}} V_i \varphi_i$ and $w_h := \sum_{i \in \mathcal{V}} W_i \varphi_i \in V_h$. We write $v_h^{\pm} := \sum_{i \in \mathcal{V}} V_i^{\pm} \varphi_i$ and $w_h^{\pm} := \sum_{i \in \mathcal{V}} W_i^{\pm} \varphi_i \in V_h$. For all $i \in \mathcal{V}$, we have $(V_i^{-} - W_i^{-})(V_i^{+} - W_i^{+}) \ge 0$. Indeed, for all $Z \in [-1, 1]$, since $V_i^{+}$ and $W_i^{+}$ are the projections of $V_i$ and $W_i$ into the convex set $[-1, 1]$, we have

$$(V_i - V_i^{+})(Z - V_i^{+}) \le 0 \quad \text{and} \quad (W_i - W_i^{+})(Z - W_i^{+}) \le 0.$$

Then picking $Z = W_i^+$ and $Z = V_i^+$ in the first and second inequalities, respectively, gives $(V_i^{-} - W_i^{-})(V_i^{+} - W_i^{+}) \ge 0$. Combining this bound with the fact that $\frac{h_i^d}{\tau} + h_i^{d-2} + \frac{h_i^d}{\epsilon^2} > 0$, we infer (14a). The second inequality is a direct application of (14a), by picking $w_h \in V_h^+$ so that $w_h^+ = w_h$ and $w_h^- = 0$. $\square$

The next lemma will be applied for proving the well-posedness of (12).

**Lemma 3.2 (Bound on stabilization).** There exists two constants $C_1, C_2$, independent of $\epsilon, h$ and $\tau$, such that, for all $v_h, w_h \in V_h$,

$$C_1 \alpha \left( \left( \frac{1}{\tau} + \frac{1}{\epsilon^2} \right) \|v_h\|^2 + \|\nabla v_h\|^2 \right) \le s_h(v_h, v_h), \quad (15a)$$

$$s_h(v_h, w_h) \le C_2 \alpha \left( \frac{1}{\tau} + \frac{1}{h_{\min}^2} + \frac{1}{\epsilon^2} \right) \|v_h\| \|w_h\|. \quad (15b)$$

*Proof.* Let $v_h := \sum_{i \in \mathcal{V}} V_i \varphi_i$ and $w_h := \sum_{i \in \mathcal{V}} W_i \varphi_i$. The inequality (15b) is a direct application of the equivalence between the norms $\|v_h\|^2$ and $\sum_{i \in \mathcal{V}} h_i^d V_i^2$ (see, for example, [15, Remark 66.5] and [26, Chapter 15]). For the inequality (15a), the norm equivalence between $\|v_h\|^2$ and $\sum_{i \in \mathcal{V}} h_i^d V_i^2$ implies that

$$\left(\frac{1}{\tau} + \frac{1}{\epsilon^2}\right) \|v_h\|^2 \le C \sum_{i \in \mathcal{V}} \left(\frac{h_i^d}{\tau} + \frac{h_i^d}{\epsilon^2}\right) V_i^2. \quad (16)$$
For the gradient term, we invoke the following property:

$$ \|\nabla v_h\|^2 \le C \sum_{i \in \mathcal{V}} \sum_{j \in \mathcal{I}(i)} |\mathbf{x}_i - \mathbf{x}_j|^{d-2} (V_i - V_j)^2, $$

with $\mathcal{I}(i) := \{j \in \mathcal{V} | \int_D |\varphi_i \varphi_j| > 0\}$ whose cardinal number is uniformly bounded owing to the mesh shape-regularity. This implies that

$$ \|\nabla v_h\|^2 \le C \sum_{i \in \mathcal{V}} h_i^{d-2} V_i^2. \qquad (17) $$

Combining the bounds (16)–(17) completes the proof. □

With the previous lemma, we can prove the existence and uniqueness of the numerical solution $u_h^n \in V_h$ to (12) for all $n \in \mathcal{N}$.

**Lemma 3.3** (Discrete well-posedness). Assume that $\alpha > \frac{1}{C_1}$, where $C_1$ stems from (15a), and that the time step is small enough so that $\tau < \frac{\epsilon^2}{4}$. Then (12) has a unique solution $u_h^n \in V_h$ for all $n \in \mathcal{N}$.

**Proof.** We define the map $F_h : V_h \to V_h$ by

$$ (F_h(v_h), w_h) := (v_h^+, w_h) + \tau(\nabla v_h^+, \nabla w_h) + \tau\epsilon^{-2}(f(v_h^+), w_h) + \tau s_h(v_h^-, w_h), \quad (18) $$

for all $w_h \in V_h$. We are going to show that $F_h$ is continuous, strictly monotone and coercive, so that well-posedness follows from the theory of finite-dimensional monotone operators (see [24, Corollary 10.42]). For all $v_h, w_h \in V_h$, we denote

$$ \theta_h := v_h - w_h, \quad \theta_h^{\oplus} := v_h^+ - w_h^+, \quad \theta_h^{\ominus} := v_h^- - w_h^- $$

i) The continuity of $F_h$ is a direct consequence of the Cauchy–Schwarz inequality and the upper bound of $s_h$ on (15b). Indeed, for all $z_h \in V_h$, we have

$$ \begin{aligned} |(F_h(v_h) - F_h(w_h), z_h)| &\le \|\theta_h^{\oplus}\|\|z_h\| + \tau\|\nabla\theta_h^{\oplus}\|\|\nabla z_h\| \\ &\quad + \tau\epsilon^{-2}\|f(v_h^+) - f(w_h^+)\|\|z_h\| \\ &\quad + C_2\alpha\left(1 + \frac{\tau}{h_{\min}^2} + \frac{\tau}{\epsilon^2}\right)\|\theta_h^{\ominus}\|\|z_h\|. \end{aligned} $$

Then we invoke the Lipschitz continuity of $f$, the boundedness of functions in $V_h^+$, the inverse inequality $\|\nabla v_h\| \le C h_{\min}^{-1} \|v_h\|$, and the fact that $z_h \in V_h$ is arbitrary, to infer that

$$ \|F_h(v_h) - F_h(w_h)\| \le C \max(1, \alpha) \left(1 + \frac{\tau}{h_{\min}^2} + \frac{\tau}{\epsilon^2}\right) \left(\|\theta_h^{\oplus}\| + \|\theta_h^{\ominus}\|\right). $$

We observe that $\|\theta_h^{\oplus}\| + \|\theta_h^{\ominus}\| \le C\|\theta_h\|$, which follows from the equivalence of $\|v_h\|^2$ and $\sum_{i \in \mathcal{V}} h_i^d V_i^2$, $|V_i^+ - W_i^+| \le |V_i - W_i|$ and $|V_i^- - W_i^-| \le |V_i^+ - W_i^+| + |V_i - W_i|$, for all $i \in \mathcal{V}$. Hence, as $\|\theta_h\| \to 0$, $\|F_h(v_h) - F_h(w_h)\| \to 0$. This proves the continuity of $F_h$.

ii) The strict monotonicity of $F_h$ is proved as follows: We start with the identity

$$ (F_h(v_h) - F_h(w_h), \theta_h) = (\theta_h^{\oplus}, \theta_h) + \tau(\nabla\theta_h^{\oplus}, \nabla\theta_h) + \frac{\tau}{\epsilon^2}(f(v_h^+) - f(w_h^+), \theta_h) + \tau s_h(\theta_h^{\ominus}, \theta_h). \quad (19) $$
$$ (\theta_h^{\oplus}, \theta_h) + \tau(\nabla\theta_h^{\oplus}, \nabla\theta_h) \ge \frac{1}{2}\|\theta_h^{\oplus}\|^2 - \frac{1}{2}\|\theta_h^{\ominus}\|^2 + \frac{\tau}{2}\|\nabla\theta_h^{\oplus}\|^2 - \frac{\tau}{2}\|\nabla\theta_h^{\ominus}\|^2. $$

$$ \begin{aligned} \frac{\tau}{\epsilon^2} (f(v_h^+) - f(w_h^+), \theta_h) &= \frac{\tau}{\epsilon^2} (f(v_h^+) - f(w_h^+), \theta_h^{\oplus} + \theta_h^{\ominus}) \\ &= \frac{\tau}{\epsilon^2} \left( (f(v_h^+) - f(w_h^+), \theta_h^{\oplus}) + (f(v_h^+) - f(w_h^+), \theta_h^{\ominus}) \right) \\ &\ge \frac{\tau}{\epsilon^2} \left( \inf_{s \in [-1,1]} f'(s) \|\theta_h^{\oplus}\|^2 - \|f'\|_{L^\infty([-1,1])} \|\theta_h^{\oplus}\| \|\theta_h^{\ominus}\| \right) \\ &\ge \frac{\tau}{\epsilon^2} \left( \inf_{s \in [-1,1]} f'(s) \|\theta_h^{\oplus}\|^2 - \frac{\|f'\|_{L^\infty([-1,1])}}{2} (\|\theta_h^{\oplus}\|^2 + \|\theta_h^{\ominus}\|^2) \right) \\ &\ge -\frac{\tau}{\epsilon^2} (2\|\theta_h^{\oplus}\|^2 + \|\theta_h^{\ominus}\|^2). \end{aligned} $$

$$ \tau s_h(\theta_h^{\ominus}, \theta_h) = \tau s_h(\theta_h^{\ominus}, \theta_h^{\oplus}) + \tau s_h(\theta_h^{\ominus}, \theta_h^{\ominus}) \ge \tau s_h(\theta_h^{\ominus}, \theta_h^{\ominus}) \ge C_1 \alpha \left( \left(1 + \frac{\tau}{\epsilon^2}\right) \|\theta_h^{\ominus}\|^2 + \tau \|\nabla\theta_h^{\ominus}\|^2 \right). $$

$$ \begin{aligned} (F_h(v_h) - F_h(w_h), \theta_h) &\ge \left(\frac{1}{2} - \frac{2\tau}{\epsilon^2}\right) \|\theta_h^{\oplus}\|^2 + \frac{\tau}{2} \|\nabla\theta_h^{\oplus}\|^2 \\ &\quad + \left(\frac{\tau}{\epsilon^2}(C_1\alpha - 1) + (C_1\alpha - \frac{1}{2})\right) \|\theta_h^{\ominus}\|^2 + \tau \left(C_1\alpha - \frac{1}{2}\right) \|\nabla\theta_h^{\ominus}\|^2, \end{aligned} \quad (20) $$

$$ \begin{aligned} (F_h(v_h), v_h) &= (F_h(v_h) - F_h(w_h), \theta_h) \\ &\ge \left(\frac{1}{2} - \frac{2\tau}{\epsilon^2}\right) \|\theta_h^{\oplus}\|^2 + \left(\frac{\tau}{\epsilon^2}(C_1\alpha - 1) + (C_1\alpha - \frac{1}{2})\right) \|\theta_h^{\ominus}\|^2 \\ &= \left(\frac{1}{2} - \frac{2\tau}{\epsilon^2}\right) \|v_h^+\|^2 + \left(\frac{\tau}{\epsilon^2}(C_1\alpha - 1) + (C_1\alpha - \frac{1}{2})\right) \|v_h^-\|^2 \\ &\ge \frac{1}{2} \min \left(\frac{1}{2} - \frac{2\tau}{\epsilon^2}, \frac{\tau}{\epsilon^2}(C_1\alpha - 1) + (C_1\alpha - \frac{1}{2})\right) \|v_h\|^2. \end{aligned} $$

$$ \frac{1}{\tau}(u_h^{n,+} - u_h^{n-1,+}, v_h - u_h^{n,+}) + (\nabla u_h^{n,+}, \nabla(v_h - u_h^{n,+})) + \frac{1}{\epsilon^2}(f(u_h^{n,+}), v_h - u_h^{n,+}) \ge 0. \quad (21) $$

<footnote n="178">178 For the first two terms on the right-hand side, we invoke the identity $(a_h, a_h+b_h) \ge \frac{1}{2} \|a_h\|^2 - \frac{1}{2} \|b_h\|^2$</footnote>

<footnote n="179">179 for all $a_h, b_h \in V_h$, which gives</footnote>

<footnote n="180">180 In order to bound the third term on the right-hand side of (19), we observe that $\inf_{s \in [-1,1]} f'(s) = -1$ and $\|f'\|_{L^\infty([-1,1])} = 2$. Then the Cauchy–Schwarz and Young’s inequalities lead to</footnote>

<footnote n="182">182 Finally, the stabilization term on the right-hand side of (19) is bounded by using (14a) and (15a):</footnote>

<footnote n="183">183 Rearranging the terms gives</footnote>

<footnote n="184">184 which is positive for $\alpha > \frac{1}{C_1}$, $\tau < \frac{\epsilon^2}{4}$ and $v_h \ne w_h$. This proves the strict monotonicity.</footnote>

<footnote n="185">185 iii) For the coercivity of $F_h$, we pick $w_h = 0 \in V_h^+$ in (20) and use the inequality $\|v_h\|^2 \le 2\|v_h^+\|^2 + 2\|v_h^-\|^2$ to infer that</footnote>

<footnote n="187">187 Therefore, $(F_h(v_h), v_h)/\|v_h\| \to \infty$ as $\|v_h\| \to \infty$ for $\alpha > \frac{1}{C_1}$, $\tau < \frac{\epsilon^2}{4}$. This completes the proof. $\square$</footnote>

<footnote n="188">188 Owing to the following lemma, our nonlinear FEM scheme can also be considered as a solver for a constrained variational inequality.</footnote>

<footnote n="189">189 Lemma 3.4 (Constrained variational inequality). Let the assumptions in Lemma 3.3 hold true. Then for all $n \in \mathcal{N}$, $u_h^{n,+} \in V_h^+$ satisfies the following variational inequality for all $v_h \in V_h^+$:</footnote>
In addition, $u_h^{n,+}$ is the unique minimizer in $V_h^+$ of the functional $\tilde{\mathcal{J}}_\epsilon^n: V_h^+ \to \mathbb{R}^+$ defined as

$$ \tilde{\mathcal{J}}_{\epsilon}^{n}(v_h) := \frac{1}{2\tau} \|v_h - u_h^{n-1,+}\|^2 + \mathcal{J}_{\epsilon}(v_h). \quad (22) $$

Proof. Let $v_h \in V_h^+$ and consider $v_h - u_h^{n,+}$ as the test function in (12). This leads to

$$ \frac{1}{\tau}(u_h^{n,+} - u_h^{n-1,+}, v_h - u_h^{n,+}) + (\nabla u_h^{n,+}, \nabla(v_h - u_h^{n,+})) + \frac{1}{\epsilon^2}(f(u_h^{n,+}), v_h - u_h^{n,+}) = s_h(u_h^{n,-}, u_h^{n,+} - v_h). $$

Owing to (14b), $s_h(u_h^{n,-}, u_h^{n,+} - v_h) \ge 0$. This proves (21). Now we show that $u_h^{n,+}$ is the unique minimizer of $\tilde{\mathcal{J}}_\epsilon^n$. For all $v_h, \tilde{v}_h \in V_h^+$ and $\delta \in \mathbb{R}$, setting $\delta w_h := v_h - \tilde{v}_h$ and observing that $V_h^+$ is a convex set, we have

$$ \tilde{\mathcal{J}}_{\epsilon}^{n}(v_h + \delta w_h) = \tilde{\mathcal{J}}_{\epsilon}^{n}(v_h) + \delta \left( \frac{1}{\tau}(v_h - u_h^{n-1,+}, w_h) + (\nabla v_h, \nabla w_h) + \frac{1}{\epsilon^2}(f(v_h), w_h) \right) \\ + \delta^2 \left( \frac{1}{2\tau} \|w_h\|^2 + \frac{1}{2} \|\nabla w_h\|^2 + \frac{1}{2\epsilon^2} \int_D f'(v_h)(w_h)^2 \right) + O(\delta^3), $$

where the terms of order $O(\delta)$ and $O(\delta^2)$ correspond to the first- and second-order variations of the discrete energy functional $\tilde{\mathcal{J}}_\epsilon^n$, respectively. Noticing that $\|f'(v_h)\|_{L^\infty} = \|3(v_h)^2 - 1\|_{L^\infty} \le 2$ for all $v_h \in V_h^+$ and since $\tau < \frac{\epsilon^2}{4}$, we have

$$ \frac{1}{2\tau} \|w_h\|^2 + \frac{1}{2} \|\nabla w_h\|^2 + \frac{1}{2\epsilon^2} \int_D f'(v_h)(w_h)^2 \ge \frac{1}{2\tau} \|w_h\|^2 - \frac{1}{\epsilon^2} \|w_h\|^2 > 0, $$

for all $v_h \in V_h^+$ and $w_h \in V_h^+ \setminus \{0\}$. This implies that $\tilde{\mathcal{J}}_\epsilon^n$ is strictly convex. Hence, $u_h^{n,+}$ is the unique minimizer in $V_h^+$ of $\tilde{\mathcal{J}}_\epsilon^n$. □

Our last result shows that the scheme (12) inherits the energy stability from the continuous level.

**Corollary 3.5 (Stability).** Let the assumptions in Lemma 3.3 be satisfied. Then the following stability property holds true:

$$ \max_{n \in \mathcal{N}} \|u_h^{n,+}\|^2 + \sum_{n \in \mathcal{N}} \|u_h^{n,+} - u_h^{n-1,+}\|^2 + 2 \sum_{n \in \mathcal{N}} \tau \|\nabla u_h^{n,+}\|^2 \le \|u_h^0\|^2 + C\epsilon^{-2}. \quad (23) $$

In addition, for the total energy $\mathcal{J}_\epsilon$ defined in (3c), we have

$$ \max_{n \in \mathcal{N}} \mathcal{J}_{\epsilon}(u_h^{n,+}) + \frac{7}{8} \sum_{n \in \mathcal{N}} \tau \left\| \frac{u_h^{n,+} - u_h^{n-1,+}}{\tau} \right\|^2 \le \mathcal{J}_{\epsilon}(u_h^0). \quad (24) $$

Proof. In order to prove (23), we choose $v_h = u_h^{n,+}$ in (12) to find that

$$ \frac{1}{2} (\|u_h^{n,+}\|^2 - \|u_h^{n-1,+}\|^2 + \|u_h^{n,+} - u_h^{n-1,+}\|^2) + \tau \|\nabla u_h^{n,+}\|^2 \\ = -\tau \epsilon^{-2} (f(u_h^{n,+}), u_h^{n,+}) - \tau s_h(u_h^{n,-}, u_h^{n,+}), \quad \forall n \in \mathcal{N}. $$

The first term on the right-hand side is bounded by $C\tau\epsilon^{-2}$ since $\|u_h^{n,+}\|_{L^\infty} \le 1$. The second term is non-positive owing to (14b) and $0 \in V_h^+$. Summing over $n \in \mathcal{N}$ and since $\sum_{n \in \mathcal{N}} \tau = T = O(1)$
by assumption, we get the expected bound (23). To prove (24), we choose $v_h := \frac{u_h^{n,+} - u_h^{n-1,+}}{\tau}$ in (12) to infer that

$$ \frac{1}{2\tau} \|\nabla u_h^{n,+}\|^2 - \frac{1}{2\tau} \|\nabla u_h^{n-1,+}\|^2 + \frac{1}{2\tau} \|\nabla(u_h^{n,+} - u_h^{n-1,+})\|^2 + \left\| \frac{u_h^{n,+} - u_h^{n-1,+}}{\tau} \right\|^2 \\ + \frac{1}{\epsilon^2} \left( f(u_h^{n,+}), \frac{u_h^{n,+} - u_h^{n-1,+}}{\tau} \right) + \frac{1}{\tau} s_h(u_h^{n,-}, u_h^{n,+} - u_h^{n-1,+}) = 0, \quad \forall n \in \mathcal{N}. \quad (25) $$

The stabilization term in the second line is non-negative owing to (14b). Now we control the first term in the second line. For all $a, b \in [-1, 1]$, the second-order Taylor expansion of $F$ and $\inf_{s \in [-1, 1]} f'(s) = -1$ imply that

$$ F(b) - F(a) \le f(b)(b-a) + \frac{1}{2}(b-a)^2. $$

Since $u_h^{n,+}$ and $u_h^{n-1,+}$ take values in $[-1, 1]$, we infer that

$$ \left(f(u_h^{n,+}), \frac{u_h^{n,+} - u_h^{n-1,+}}{\tau}\right) \ge \int_D \frac{F(u_h^{n,+}) - F(u_h^{n-1,+})}{\tau} - \frac{\tau}{2} \left\| \frac{u_h^{n,+} - u_h^{n-1,+}}{\tau} \right\|^2. $$

Inserting this estimate in (25), and summing over $n \in \mathcal{N}$ proves (24), since $1 - \frac{\tau}{2\epsilon^2} \ge \frac{7}{8}$. $\square$

**Remark 3.6 (Convergence with exponential $\epsilon$-dependence).** Under the assumption $\tau < \frac{\epsilon^2}{4}$, the convergence of the scheme can be proved by using the standard Grönwall inequality, but the final error bound will exhibit exponential dependence on $\epsilon$. For the sake of brevity, we omit this convergence result and its proof and refer to [7, Theorem 6.3].

# 4 Error analysis

In this section, we establish an error estimate for our scheme with only polynomial dependence on $\epsilon^{-1}$. The main results are Theorem 4.4 and Corollary 4.5.

## 4.1 Preliminary tools

In the classical parabolic framework, the error is decomposed as $u - u_h = (u - \Pi_h^E(u)) + (\Pi_h^E(u) - u_h)$, where $\Pi_h^E(u)$ is the elliptic or Ritz projection. This projection is introduced in order to obtain an optimal convergence rate in the $L^\infty(L^2)$-norm. However, this simple choice does not fit well with our scheme owing to the presence of the stabilization term. To solve this issue, we introduce a nonlinear bound-preserving map which is the counterpart of the Ritz projection. For all $v \in H^1(D)$, we consider the following problem: Find $\tilde{v}_h \in V_h$ by solving the equation

$$ (\tilde{v}_h^+ - v, w_h) + (\nabla(\tilde{v}_h^+ - v), \nabla w_h) + \tilde{s}_h(\tilde{v}_h - \tilde{v}_h^+, w_h) = 0, \quad \forall w_h \in V_h, \quad (26) $$

where $\tilde{v}_h^+ := P_h^+(\tilde{v}_h) \in V_h^+$ and $\tilde{s}_h(v_h, w_h) := \tilde{\alpha} \sum_{i \in \mathcal{V}} (h_i^d + h_i^{d-2})V_i W_i$ with a sufficiently large constant $\tilde{\alpha} = O(1)$. We notice that the properties stated in Lemma 3.1 also hold true for the bilinear form $\tilde{s}_h$. Using $\tilde{v}_h^+$, we can define the following nonlinear Ritz map and derive its fundamental properties.
**Lemma 4.1** (Nonlinear Ritz map). We define the map $\Pi_h^+ : H^1(D) \to V_h^+$ by $\Pi_h^+(v) := \tilde{v}_h^+$, the bound-preserving solution to (26), for all $v \in H^1(D)$. Equivalently, we have

236

$$ \Pi_h^+(v) = \arg\min_{v_h \in V_h^+} \|v_h - v\|_{H^1}. \qquad (27) $$

In addition, if $v \in H^2(D)$, $\Pi_h^+(v)$ has optimal approximation properties in the $H^1$-norm.

237

*Proof.* The well-posedness to (26) follows from [4, Theorem 3.2], and the minimization property (27) is proved in [4, Theorem 3.7]. Finally, if $v \in H^2(D)$, we can pick the Lagrange interpolation of $v$ as a candidate in (27), and conclude that $\Pi_h^+(v)$ has optimal approximation properties in the $H^1$-norm. $\square$

238

239

240

241

242

243

244

245

**Remark 4.2** (Optimality in the $L^2$-norm). Although this nonlinear Ritz map $\Pi_h^+$ is optimal in the $H^1$-norm, it is still an open question whether $\Pi_h^+$ has optimal approximation properties in the $L^2$-norm, owing to the non-symmetry in its construction. However, it is observed numerically in [4, 3] that the error in the $L^2$-norm is also optimal.

246

247

248

The classical discrete Grönwall inequality delivers error estimates with exponential dependence on $\epsilon^{-1}$, which is not acceptable in practice. Therefore, the following sharper Grönwall-type inequality [6, Lemma 6.2] is needed.

249

250

**Lemma 4.3** (Shaper discrete Grönwall inequality). Suppose that the non-negative real sequences $(\beta_n)_{n \in \mathcal{N}}$, $(w_j^n)_{n \in \mathcal{N}}$ for all $j = 1, 2, 3$, and the positive numbers $A, B, \chi$ satisfy

$$ w_1^n + \sum_{k \in \{1:n\}} \tau w_2^k \le A + \sum_{k \in \{1:n\}} \tau \beta_k w_1^k + \sum_{k \in \{1:n-1\}} \tau w_3^k, \qquad (28a) $$

$$ \sum_{k \in \{1:n-1\}} \tau w_3^k \le B \left\{ \max_{k \in \{1:n-1\}} (w_1^k) \chi \right\} \sum_{k \in \{1:n-1\}} \tau (w_1^k + w_2^k), \qquad (28b) $$

for all $n \in \mathcal{N}$. Assume that

251

$$ \max_{n \in \mathcal{N}} \tau \beta_n \le \frac{1}{2}, \qquad (29a) $$

$$ 8A \exp \left\{ \frac{2+2\chi}{\chi} \sum_{n \in \mathcal{N}} \tau \beta_n \right\} \le (8B(1+T))^{-\frac{1}{\chi}}. \qquad (29b) $$

Then we have

252

$$ \max_{n \in \mathcal{N}} w_1^n + \sum_{n \in \mathcal{N}} \tau w_2^n \le 8A \exp \left\{ 2 \sum_{n \in \mathcal{N}} \tau \beta_n \right\}. \qquad (30) $$

## 4.2 Error estimates

In this section, we derive an error estimate for our scheme with the help of the tools introduced above. We define $u^n(\cdot) := u(\cdot, t^n)$ for all $n \in \{0:N\}$, and the global piecewise constant in time solution $u_h^+$ by $u_h^+(\cdot, t)|_{J^n} := u_h^{n,+}(\cdot)$ for all $n \in \mathcal{N}$. We also introduce the error decomposition

254

255

256

$$ e_h^n := u^n - u_h^{n,+} = (u^n - \Pi_h^+(u^n)) + (\Pi_h^+(u^n) - u_h^{n,+}) := \eta_h^n + \xi_h^n, \quad \forall n \in \{0:N\}. \qquad (31) $$

Our error bound is devised under the assumption that the following quantity is finite:

257

$$ A_{h,\tau,\epsilon} := \|\tilde{h}\partial_t \nabla^2 u\|_{L^2(L^2)}^2 + \|\tau \partial_{tt} u\|_{L^2(L^2)}^2 + \epsilon^{-4} \|\tilde{h} \nabla^2 u\|_{L^2(L^2)}^2, \qquad (32) $$

where the piecewise constant function $\tilde{h}$ is defined by $\tilde{h}|_K = h_K$ for all $K \in \mathcal{T}_h$. A sufficient condition for $A_{h,\tau,\epsilon}$ to be finite is given in Lemma 2.2. Now we present the error estimate.

258

259
**Theorem 4.4 (Error bound).** Let the assumptions in Lemma 2.2, 2.4 and 3.3 hold true. Assume that $u_h^0 := \Pi_h^+(u_0)$ and that $\tau$ and $\tilde{h}$ are small enough so that the following conditions hold true:

$$A_{h,\tau,\epsilon} \le C\epsilon^{8+24\kappa} \exp(-24C_\lambda), \quad (33a)$$

$$\tau \le \min \left\{ \frac{\epsilon^2}{384}, \frac{1}{18 + 8 \max_{n \in \mathcal{N}} \lambda(t^n)}, C\epsilon^2 \|\partial_t u\|_{L^2(H^2)} \right\}, \quad (33b)$$

with the constants $\kappa, C_\lambda$ from Assumption 2.4. Then we have the following error bound:

$$\max_{n \in \mathcal{N}} \|\xi_h^n\|^2 + \epsilon^2 \sum_{n \in \mathcal{N}} \tau \|\nabla \xi_h^n\|^2 \le C \frac{\exp(8C_\lambda)}{\epsilon^{8\kappa}} A_{h,\tau,\epsilon}. \quad (34)$$

**Proof.** The proof is divided into three steps. In Step 1, we derive the error equation and perform some simple inequalities. In Step 2, we control the nonlinear terms. In Step 3, we apply the shaper discrete Grönwall inequality from Lemma 4.3 to conclude the proof.

**Step 1:** For all $n \in \mathcal{N}$, we define the consistency error $\mathcal{R}^n$ as the linear form on $V_h$ such that

$$\langle \mathcal{R}^n, v_h \rangle := \left( \frac{u^n - u^{n-1}}{\tau}, v_h \right) + (\nabla u^n, \nabla v_h) + \frac{1}{\epsilon^2} (f(u^n), v_h), \quad \forall v_h \in V_h. \quad (35)$$

Subtracting (12) from (35) gives

$$\frac{1}{\tau}(\xi_h^n - \xi_h^{n-1}, v_h) + (\nabla \xi_h^n, \nabla v_h) + \frac{1}{\epsilon^2}(f(u^n) - f(u_h^{n,+}), v_h) = -\frac{1}{\tau}(\eta_h^n - \eta_h^{n-1}, v_h) - (\nabla \eta_h^n, \nabla v_h) + s_h(u_h^{n,-}, v_h) + \langle \mathcal{R}^n, v_h \rangle,$$

for all $v_h \in V_h$. We invoke (26) to infer that

$$-(\nabla \eta_h^n, v_h) = (\eta_h^n, v_h) - \tilde{s}_h(\tilde{u}_h^n - \Pi_h^+(u^n), v_h),$$

where $\tilde{u}_h^n \in V_h$ is the solution to (26) in the case where $v := u^n$. Then we pick $v_h := 2\xi_h^n$ in the above equalities to obtain the error equation

$$\begin{aligned} \|\xi_h^n\|^2 - \|\xi_h^{n-1}\|^2 + \|\xi_h^n - \xi_h^{n-1}\|^2 + 2\tau \|\nabla \xi_h^n\|^2 + \frac{2\tau}{\epsilon^2}(f(u^n) - f(u_h^{n,+}), e_h^n) \\ = 2(\eta_h^n - \eta_h^{n-1}, \xi_h^n) + 2\tau(\eta_h^n, \xi_h^n) \\ - 2\tau \tilde{s}_h(\tilde{u}_h^n - \Pi_h^+(u^n), \xi_h^n) + 2\tau s_h(u_h^{n,-}, \xi_h^n) \\ + 2\tau \langle \mathcal{R}^n, \xi_h^n \rangle + \frac{2\tau}{\epsilon^2}(f(u^n) - f(u_h^{n,+}), \eta_h^n). \end{aligned} \quad (36)$$

The terms in the second line are bounded by Young's inequality:

$$2(\eta_h^n - \eta_h^{n-1}, \xi_h^n) + 2\tau(\eta_h^n, \xi_h^n) \le \int_{J^n} \|\partial_t \eta_h^n\|^2 + \tau\|\eta_h^n\|^2 + 2\tau\|\xi_h^n\|^2.$$

The stabilization terms on the third line are non-positive owing to (14b). Indeed, since $u_h^{n,+}$ and $\Pi_h^+(u^n)$ (equal to $\tilde{u}_h^{n,+}$ by definition) belong to $V_h^+$ and setting $\tilde{u}_h^{n,-} := \tilde{u}_h^n - \tilde{u}_h^{n,+}$, we find that

$$-\tilde{s}_h(\tilde{u}_h^n - \Pi_h^+(u^n), \xi_h^n) = \tilde{s}_h(\tilde{u}_h^{n,-}, u_h^{n,+} - \tilde{u}_h^{n,+}) \le 0, \quad s_h(u_h^{n,-}, \xi_h^n) = s_h(u_h^{n,-}, \Pi_h^+(u^n) - u_h^{n,+}) \le 0.$$

Now we control the consistency term. The consistency error is such that

$$\langle \mathcal{R}^n, \xi_h^n \rangle = \left( \frac{u^n - u^{n-1}}{\tau} - \partial_t u^n, \xi_h^n \right).$$
The Cauchy–Schwarz and Young’s inequalities lead to

$$
\begin{align*}
2\tau |\langle \mathcal{R}^n, \xi_h^n \rangle| &\le \tau \left\| \partial_t u^n - \frac{u^n - u^{n-1}}{\tau} \right\|^2 + \tau \|\xi_h^n\|^2 \\
&\le \frac{1}{\tau} \left\| \int_{J_n} (t - t^{n-1}) \partial_{tt} u \right\|^2 + \tau \|\xi_h^n\|^2 \\
&\le \tau^2 \int_{J_n} \|\partial_{tt} u\|^2 + \tau \|\xi_h^n\|^2.
\end{align*}
$$

Finally, to bound the last term on the right-hand side of the error equation (36), we invoke Young’s

inequality, and the fact that $u^n$ and $u_h^{n,+}$ both take values in $[-1, 1]$, and $\|f'\|_{L^\infty([-1,1])} = 2$ to

infer that

$$
\begin{align*}
\frac{2\tau}{\epsilon^2} (f(u^n) - f(u_h^{n,+}), \eta_h^n) &\le \frac{2\tau}{\epsilon^2} \|f'\|_{L^\infty([-1,1])} |(\eta_h^n + \xi_h^n, \eta_h^n)| \\
&\le \frac{4\tau}{\epsilon^2} \left( \|\eta_h^n\|^2 + \frac{1}{2} (\epsilon^2 \|\xi_h^n\|^2 + \frac{1}{\epsilon^2} \|\eta_h^n\|^2) \right) \\
&\le \frac{2\tau}{\epsilon^2} \left( 2 + \frac{1}{\epsilon^2} \right) \|\eta_h^n\|^2 + 2\tau \|\xi_h^n\|^2.
\end{align*}
$$

Step 2: The goal of this step is to bound from below the nonlinear term on the left-hand side

of (36). For all $n \in \mathcal{N}$, we consider the decomposition

$$
\frac{2\tau}{\epsilon^2} (f(u^n) - f(u_h^{n,+}), e_h^n) = \frac{2\tau(1-\epsilon^2)}{\epsilon^2} (f(u^n) - f(u_h^{n,+}), e_h^n) + 2\tau(f(u^n) - f(u_h^{n,+}), e_h^n).
$$

The second term on the right-hand side is bounded from below using $\inf_{s \in [-1,1]} f'(s) = -1$ as

follows:

$$
2\tau(f(u^n) - f(u_h^{n,+}), e_h^n) \geq -2\tau\|e_h^n\|^2 \geq -4\tau(\|\eta_h^n\|^2 + \|\xi_h^n\|^2).
$$

For the first term on the right-hand side, we invoke the principal eigenvalue defined in (4) and the

inequality $(f(a) - f(b), a - b) \ge f'(a)(a - b)^2 - 3|a - b|^3$ for all $a, b \in [-1, 1]$ to infer that

$$
\begin{align*}
\frac{1}{\epsilon^2} (f(u^n) - f(u_h^{n,+}), e_h^n) & \\
\ge & \frac{1}{\epsilon^2} (f'(u^n) e_h^n, e_h^n) - \frac{3}{\epsilon^2} \|e_h^n\|_{L^3}^3 \\
\ge & -\lambda(t^n) \|e_h^n\|^2 - \|\nabla e_h^n\|^2 - \frac{12}{\epsilon^2} (\|\eta_h^n\|_{L^3}^3 + \|\xi_h^n\|_{L^3}^3) \\
\ge & -2\lambda(t^n) (\|\eta_h^n\|^2 + \|\xi_h^n\|^2) - \left(1 + \frac{1}{\gamma}\right) \|\nabla \eta_h^n\|^2 - (1 + \gamma) \|\nabla \xi_h^n\|^2 - \frac{12}{\epsilon^2} (\|\eta_h^n\|_{L^3}^3 + \|\xi_h^n\|_{L^3}^3),
\end{align*}
$$

with a constant $\gamma > 0$ yet to be determined. We pick $\gamma := \frac{\epsilon^2}{2(1-\epsilon^2)}$ and multiply the above inequality

by $2\tau(1 - \epsilon^2)$ to get

$$
\begin{align*}
\frac{2\tau(1 - \epsilon^2)}{\epsilon^2} (f(u^n) - f(u_h^{n,+}), e_h^n) + 2\tau \|\nabla \xi_h^n\|^2 & \\
&\geq -4(1 - \epsilon^2)\tau\lambda(t^n)(\|\eta_h^n\|^2 + \|\xi_h^n\|^2) - \frac{2(2 - \epsilon^2)(1 - \epsilon^2)}{\epsilon^2}\tau \|\nabla \eta_h^n\|^2 \\
&\quad + \tau\epsilon^2 \|\nabla \xi_h^n\|^2 - \frac{24\tau(1 - \epsilon^2)}{\epsilon^2} (\|\eta_h^n\|_{L^3}^3 + \|\xi_h^n\|_{L^3}^3).
\end{align*}
$$
287 Gathering all the estimates and using that $1 - \epsilon^2 \le 1$, $2 - \epsilon^2 \le 2$, $2 + \frac{1}{\epsilon^2} \le \frac{3}{\epsilon^2}$, $\|\eta_h^n\|_{L^\infty} \le 2$,

288 $\frac{1-\epsilon^2}{\epsilon^2} \le \frac{1}{\epsilon^4}$ to simplify some factors, we arrive at

$$
\begin{aligned}
\|\xi_h^n\|^2 + \|\xi_h^n - \xi_h^{n-1}\|^2 + \tau\epsilon^2 \|\nabla\xi_h^n\|^2 &\le \|\xi_h^{n-1}\|^2 + \int_{J^n} \|\partial_t\eta_h^n\|^2 + \tau^2 \int_{J^n} \|\partial_{tt}u\|^2 \\
&\quad + \left(5 + \frac{54}{\epsilon^4} + 4\lambda(t^n)\right) \tau \|\eta_h^n\|^2 + \frac{4\tau}{\epsilon^2} \|\nabla\eta_h^n\|^2 \\
&\quad + \left(9 + 4\lambda(t^n)\right) \tau \|\xi_h^n\|^2 + \frac{24\tau}{\epsilon^2} \|\xi_h^n\|_{L^3}^3.
\end{aligned}
$$

289 Then, since $\|\xi_h^n - \xi_h^{n-1}\|_{L^\infty} \le \|\xi_h^n\|_{L^\infty} + \|\xi_h^{n-1}\|_{L^\infty} \le 4$, we have

$$
\|\xi_h^n\|_{L^3}^3 \le 4\|\xi_h^n - \xi_h^{n-1}\|_{L^3}^3 + 4\|\xi_h^{n-1}\|_{L^3}^3 \le 16\|\xi_h^n - \xi_h^{n-1}\|^2 + 4\|\xi_h^{n-1}\|_{L^3}^3.
$$

290 This leads to

$$
\begin{aligned}
& \|\xi_h^n\|^2 + \left(1 - \frac{384\tau}{\epsilon^2}\right) \|\xi_h^n - \xi_h^{n-1}\|^2 + \tau\epsilon^2 \|\nabla\xi_h^n\|^2 \\
\le & \|\xi_h^{n-1}\|^2 + \int_{J^n} \|\partial_t\eta_h\|^2 + \tau^2 \int_{J^n} \|\partial_{tt}u\|^2 + \left(5 + \frac{54}{\epsilon^4} + 4\lambda(t^n)\right) \tau \|\eta_h^n\|^2 + \frac{4\tau}{\epsilon^2} \|\nabla\eta_h^n\|^2 \\
& + \left(9 + 4\lambda(t^n)\right) \tau \|\xi_h^n\|^2 + \frac{96\tau}{\epsilon^2} \|\xi_h^{n-1}\|_{L^3}^3.
\end{aligned}
$$

291 **Step 3:** For all $n \in \mathcal{N}$, summing over $k \in \{1:n\}$, observing that $\xi_h^0 = 0$ owing to the choice of initial condition, and invoking the condition $1 - \frac{384\tau}{\epsilon^2} \ge 0$ (see (33b)), we infer that

$$
\|\xi_h^n\|^2 + \epsilon^2 \sum_{k \in \{1:n\}} \tau \|\nabla \xi_h^k\|^2 \le \tilde{A}_{h,\tau,\epsilon}^n + \sum_{k \in \{1:n\}} \left( \beta_k \tau \|\xi_h^k\|^2 + \frac{96\tau}{\epsilon^2} \|\xi_h^{k-1}\|_{L^3}^3 \right),
$$

293 with $\beta_k := 9 + 4\lambda(t^k)$ and

$$
\begin{aligned}
\tilde{A}_{h,\tau,\epsilon}^n &:= \sum_{k \in \{1:n\}} \left\{ \int_{J^k} \|\partial_t \eta_h^k\|^2 + \tau^2 \int_{J^k} \|\partial_{tt} u\|^2 + \left(5 + \frac{54}{\epsilon^4} + 4\lambda(t^k)\right) \tau \|\eta_h^k\|^2 + \frac{4\tau}{\epsilon^2} \|\nabla \eta_h^k\|^2 \right\} \\
&:= \sum_{i \in \{1:4\}} T_i^n,
\end{aligned}
$$

294 with obvious notation. We bound $T_i^n$ for each $i \in \{1:4\}$. Owing to the regularity Assumption 2.2

295 and the approximation property (27), we have, for all $n \in \mathcal{N}$,

$$
\sum_{k \in \{1:n\}} \int_{J^k} \|\partial_t \eta_h^k\|^2 \le C \|\tilde{h} \partial_t \nabla^2 u\|_{L^2(L^2)}^2, \quad (37a)
$$

$$
\sum_{k \in \{1:n\}} \tau \left( \|\eta_h^k\|^2 + \|\nabla \eta_h^k\|^2 \right) \le C \|\tilde{h} \nabla^2 u\|_{L^2(L^2)}^2. \quad (37b)
$$

296 Hence, we find

$$
T_1^n + T_4^n \le C \left( \|\tilde{h} \partial_t \nabla^2 u\|_{L^2(L^2)}^2 + \frac{1}{\epsilon^2} \|\tilde{h} \nabla^2 u\|_{L^2(L^2)}^2 \right).
$$

297 For $T_2^n$, we have

$$
T_2^n = \int_0^{t^n} \|\tau \partial_{tt} u\|^2 \le \|\tau \partial_{tt} u\|_{L^2(L^2)}^2.
$$
For $T_3^n$, we invoke (4) and use the regularity assumption 2.2 to find that, for all $t \in J$,

$$ 
\begin{aligned}
\lambda(t) &= \max \left\{ 0, - \inf_{v \in H^1(D), \|v\| = 1} \left( \epsilon^{-2}(f'(u(\cdot, t)), v^2) + \|\nabla v\|^2 \right) \right\} \\
&\le \max \left\{ 0, - \inf_{v \in H^1(D), \|v\| = 1} \|\nabla v\|^2 + \sup_{v \in H^1(D), \|v\| = 1} \epsilon^{-2}(f'(u(\cdot, t)), v^2) \right\} \le \frac{C}{\epsilon^2},
\end{aligned}
$$

since $\inf_{v \in H^1(D), \|v\| = 1} \|\nabla v\|^2 \ge 0$. This, with (37b), implies that

$$ T_3^n \le C \left( 1 + \frac{1}{\epsilon^2} + \frac{1}{\epsilon^4} \right) \|\tilde{h}\nabla^2 u\|_{L^2(L^2)}^2 \le \frac{C}{\epsilon^4} \|\tilde{h}\nabla^2 u\|_{L^2(L^2)}^2, $$

recall that $\epsilon \le 1$ by assumption. Gathering the bounds on $T_i^n$, for all $i \in \{1:4\}$, we infer that for all $n \in \mathcal{N}$,

$$ \tilde{A}_{h, \tau, \epsilon}^n \le C_3 A_{h, \tau, \epsilon}, $$

where the constant $C_3$ is independent of $\epsilon$, $h$ and $\tau$. From this, we infer that

$$ \|\xi_h^n\|^2 + \epsilon^2 \sum_{k \in \{1:n\}} \tau \|\nabla \xi_h^n\|^2 \le C_3 A_{h, \tau, \epsilon} + \sum_{k \in \{1:n\}} \left( \beta_k \tau \|\xi_h^k\|^2 + \frac{96\tau}{\epsilon^2} \|\xi_h^{k-1}\|_{L^3}^3 \right), \quad \forall n \in \mathcal{N}. $$

This is (28a) in Lemma 4.3 with $w_1^n := \|\xi_h^n\|^2$, $w_2^n := \epsilon^2 \|\nabla \xi_h^n\|^2$, $w_3^n := 96\epsilon^{-2} \|\xi_h^n\|_{L^3}^3$, $A := C_3 A_{h, \tau, \epsilon}$.

Then, we let $\chi := \frac{1}{2}$ and $B := 96C_S\epsilon^{-4}$, with $C_S$ from the following Sobolev embedding:

$$ \|v\|_{L^3}^3 \le \|v\| \|v\|_{L^4}^2 \le C_S \|v\| \|v\|_{H^1}^2, \quad \forall v \in H^1(D). $$

With the help of the above inequality, we infer (28b) in Lemma 4.3. It remains to verify (29a) and (29b). First, (29a) is a direct consequence of assumption (33b). With the help of the spectral control from [9, Lemma 4.2 (c)] and under the assumption $\tau \le C\epsilon^2 \|\partial_t u\|_{L^2(H^2)}$ (see (33b)), we have

$$ \left| \sum_{n \in \mathcal{N}} \tau \lambda(t^n) - \int_J \lambda \right| \le T. $$

Thus, we find

$$ \exp\left\{6 \sum_{n \in \mathcal{N}} \tau \beta_n\right\} = \exp\left\{6 \sum_{n \in \mathcal{N}} \tau(9 + 4\lambda(t^n))\right\} \le \exp\left\{78T + 24 \int_J \lambda\right\} \le C_4 \epsilon^{-24\kappa} \exp\{24C_\lambda\}, $$

with $C_4 := \exp(78T)$, and where we invoked (5). This implies that

$$ 8C_3 A_{h, \tau, \epsilon} \exp\left\{6 \sum_{n \in \mathcal{N}} \tau \beta_n\right\} \le 8C_3 C_4 A_{h, \tau, \epsilon} \epsilon^{-24\kappa} \exp\{24C_\lambda\}. $$

Hence, owing to (33b), we infer that

$$ 8C_3 C_4 A_{h, \tau, \epsilon} \epsilon^{-24\kappa} \exp\{24C_\lambda\} \le (768C_S(1+T)\epsilon^{-4})^{-2}, $$

which shows that (29b) is satisfied. Therefore, we can invoke the shaper discrete Grönwall inequality, which gives

$$ \max_{n \in \mathcal{N}} \|\xi_h^n\|^2 + \epsilon^2 \sum_{n \in \mathcal{N}} \tau \|\nabla \xi_h^n\|^2 \le 8C_3 \exp(2 \sum_{n \in \mathcal{N}} \tau \beta_n) A_{h, \tau, \epsilon} \le C \frac{\exp(8C_\lambda)}{\epsilon^{8\kappa}} A_{h, \tau, \epsilon}. $$

This completes the proof. □
Finally, the triangle inequality leads to the following error estimates for $e_h$.

**Corollary 4.5.** Let the assumptions in Theorem 4.4 be satisfied. Then we have the following error bounds:

$$ \max_{n \in \mathcal{N}} \|u^n - u_h^{n,+}\|^2 \le C \left( \frac{\exp(8C_\lambda)}{\epsilon^{8\kappa}} A_{h,\tau,\epsilon} + \|\tilde{h}\nabla^2 u\|_{L^\infty(L^2)}^2 \right), \quad (38a) $$

$$ \sum_{n \in \mathcal{N}} \tau \|\nabla(u^n - u_h^{n,+})\|^2 \le C \left( \frac{\exp(8C_\lambda)}{\epsilon^{8\kappa+2}} A_{h,\tau,\epsilon} + \|\tilde{h}\nabla^2 u\|_{L^2(L^2)}^2 \right). \quad (38b) $$

**Remark 4.6 (Improved convergence).** By assuming that $\|\eta_h\| \le C\|\tilde{h}^2\nabla^2 u\|$, i.e., that the non-linear Ritz map has optimal convergence properties in the $L^2$-norm, an optimal convergence rate in the $L^\infty(L^2)$-norm and a better dependence on $\epsilon$ can be obtained with more restrictive assumptions on $h$. This can be shown by adapting the arguments from [18]. Details are omitted for brevity.

**Remark 4.7 (Initial condition).** It is possible to choose any other optimally converging approximation of the initial condition. This will only lead to the additional contribution $\|\tilde{h}\nabla^2 u_0\|^2$ in the definition of $A_{h,\tau,\epsilon}$.

**Remark 4.8 (Explicit $\epsilon$-dependence).** Under the assumption $\kappa = 0$, i.e., if there are no topological changes during the evolution, it is proved in [13] that $\|\nabla^m \partial_t^k u\|_{L^\infty(L^2)} \le C\epsilon^{-k-m-\frac{1}{2}}$ for arbitrary $k, m \in \mathbb{N}$. This leads to a more precise bound on $A_{h,\tau,\epsilon}$. To the best of our knowledge, such bounds are not available when topological changes can occur. Moreover, assuming $\|\nabla^2 u\|_{L^2(L^2)} \le C\epsilon^{-2}$, $\|\nabla^2 u\|_{L^\infty(L^2)} \le C\epsilon^{-\frac{5}{2}}$, $\|\partial_t \nabla^2 u\|_{L^2(L^2)} \le C\epsilon^{-3}$ and $\|\partial_{tt} u\|_{L^2(L^2)} \le C\epsilon^{-2}$, a straightforward consequence of Corollary 4.5 are the following error bounds:

$$ \max_{n \in \mathcal{N}} \|u^n - u_h^{n,+}\|^2 \le C \left( \frac{\exp(8C_\lambda)}{\epsilon^{8\kappa}} \left( \frac{h^2}{\epsilon^8} + \frac{\tau^2}{\epsilon^4} \right) + \frac{h^2}{\epsilon^5} \right), \quad (39a) $$

$$ \sum_{n \in \mathcal{N}} \tau \|\nabla(u^n - u_h^{n,+})\|^2 \le C \left( \frac{\exp(8C_\lambda)}{\epsilon^{8\kappa+2}} \left( \frac{h^2}{\epsilon^8} + \frac{\tau^2}{\epsilon^4} \right) + \frac{h^2}{\epsilon^4} \right). \quad (39b) $$

# 5 Numerical experiments

In this section, we numerically verify our scheme on two test cases from the literature. The first test case involves topological changes whereas the second one does not. In addition, we compare our scheme to another bound-preserving scheme based on an operator-splitting technique. The numerical experiments are conducted with the Gridap.jl library in [27] in the Julia programming language.

## 5.1 Presentation of the test cases

Here we present two test cases with the discretization parameters considered in our simulations.

**Test case 1 ([20, Section 6, test 3]):** We set the final time $T := 0.025$, the domain $D := (-1, 1)^2$, and parameter $\epsilon \in \{0.05, 0.1, 1\}$. We enforce a homogeneous Neumann boundary condition and the initial data is

$$ u_0(x, y) := u_0^{(1)}(x, y)u_0^{(2)}(x, y), \quad (40) $$
342 with

$$
\begin{align*}
u_0^{(1)}(x, y) &:= \begin{cases}
\tanh\left(60\left((x-0.5)^2 + (y-0.4)^2 - (0.25)^2\right)\right), & \text{if } x > 0.3, y \in (-1, 1), \\
\tanh\left(60\left((y-0.4)^2 - (0.15)^2\right)\right), & \text{if } x \in [-0.3, 0.3], y \in (-1, 1), \\
\tanh\left(60\left((x+0.5)^2 + (y-0.4)^2 - (0.25)^2\right)\right), & \text{if } x < -0.3, y \in (-1, 1),
\end{cases} \\
u_0^{(2)}(x, y) &:= \begin{cases}
\tanh\left(60(x^2 + (y-0.6)^2 - (0.25)^2)\right), & \text{if } y > 0.4, x \in (-1, 1), \\
\tanh\left(60(x^2 - (0.15)^2)\right), & \text{if } y \in [-0.4, 0.4], x \in (-1, 1), \\
\tanh\left(60(x^2 + (y+0.6)^2 - (0.25)^2)\right), & \text{if } y < -0.4, x \in (-1, 1).
\end{cases}
\end{align*}
$$

343 Noting that there are topological changes during the evolution so that $\kappa \neq 0$.

344 Test case 2 ([17, test 1]): We set the final time $T := 0.032$, the domain $D := (-1, 1)^2$, and
345 parameter $\epsilon \in \{0.05, 0.1, 1\}$. We enforce a homogeneous Neumann boundary condition and the
346 initial data is

$$
u_0(x, y) := \begin{cases} \tanh \frac{d(x, y)}{\sqrt{2\epsilon}}, & \text{if } \frac{x^2}{0.36} + \frac{y^2}{0.04} \ge 1, \\ \tanh \frac{-d(x, y)}{\sqrt{2\epsilon}}, & \text{otherwise,} \end{cases} \quad (41)
$$

347 with $d(x, y) := \frac{x^2}{0.36} + \frac{y^2}{0.04} - 1$. There is no topological change during the evolution so that $\kappa = 0$.

348 Since the above test cases do not have an explicit expression for the exact solution, we construct
349 reference solutions by considering a very fine mesh and a very small time step. In practice,
350 we employ the classical implicit-Euler scheme with $P_1$-FEM as in [18] with $\tau := \frac{h}{10}$ and $h :=
351 1.56 \times 10^{-3}$, i.e., the number of cells is roughly $100 \times 2^{14} \approx 1.6$ million. Instead, the numerical
352 solutions are computed with $\tau := \frac{h}{2}$ and the number of cells is $100 \times 2^l$ with $l \in \{4, 6, 8, 10\}$.

353 5.2 Implementation details

354 In our numerical experiments, we set the weight of the stabilization to α = 1. This turns out to
355 be sufficient in practice. The same choice is considered in [4, 3, 5].

356 The nonlinear scheme (12) is implemented with a Richardson-like iteration coupled with a
357 linearized Newton iteration. At each time step $n \in \mathcal{N}$, we iteratively find $u_h^n := \lim_{k \to \infty} u_h^{n,k}$ such
358 that for all $k \in \mathbb{N}$ and all $v_h \in V_h$,

$$
a(u_h^{n,k}, v_h) = a(u_h^{n,k-1}, v_h) + \omega \left( (u_h^{n-1,+}, v_h) - a(u_h^{n,k-1,+}, v_h) - \tau s_h (u_h^{n,k-1,-}, v_h) \right), \quad (42)
$$

359 with $\omega \in \mathbb{R}^+$ a user-defined parameter, the initial guess $u_h^{n,0} := u_h^{n-1,+}$ and

$$
a(v_h, w_h) := (v_h, w_h) + \tau(\nabla v_h, \nabla w_h) + \frac{\tau}{\epsilon^2}(f(v_h), w_h) \quad \forall v_h, w_h \in V_h. \quad (43)
$$

360 Since $a(\cdot, \cdot)$ is nonlinear, we call a linearized Newton solver to replace the nonlinear function $f$ in
361 (42). This leads to an additional iteration loop where we find $u_h^{n,k} := \lim_{l \to \infty} u_h^{n,k,l}$ by solving

$$
\tilde{a}(u_h^{n,k,l-1}; u_h^{n,k,l}, w_h) = a(u_h^{n,k-1}, w_h) + \omega \left( (u_h^{n-1,+}, w_h) - a(u_h^{n,k-1,+}, w_h) - \tau s_h (u_h^{n,k-1,-}, w_h) \right), \quad (44)
$$

362 with

$$
\tilde{a}(u_h^{n,k,l-1}; v_h, w_h) := (v_h, w_h) + \tau(\nabla v_h, \nabla w_h) \\
\qquad + \frac{\tau}{\epsilon^2} \left( f(u_h^{n,k,l-1}) + f'(u_h^{n,k,l-1})(v_h - u_h^{n,k,l-1}), w_h \right), \quad (45)
$$
363 for all $l \in \mathbb{N}$ and all $v_h, w_h \in V_h$, with the initial guess $u_h^{n,k,0} := u_h^{n,k-1,+}$. The corresponding pseudo-code is presented in Algorithm 1.
364

<caption>**Algorithm 1 Numerical scheme**</caption>

**Require:** initial data $u_h^0$, mesh $\mathcal{T}_h$, time step indices $\mathcal{N}$, relaxation parameter $\omega$, tolerance $\delta$.
1: **for** $n \in \mathcal{N}$ **do**
2: &nbsp;&nbsp;&nbsp;&nbsp;$k \leftarrow 1$, $e_k \leftarrow \infty$
3: &nbsp;&nbsp;&nbsp;&nbsp;$u_h^{n,0,+} \leftarrow u_h^{n-1,+}$
4: &nbsp;&nbsp;&nbsp;&nbsp;**while** $e_k > \delta$ **do**
5: &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;$\triangleright$ Richardson-like iteration
6: &nbsp;&nbsp;&nbsp;&nbsp;**if** $l = 1$ **then**
7: &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;$u_h^{n,k,0} \leftarrow u_h^{n-1,+}$
8: &nbsp;&nbsp;&nbsp;&nbsp;**else**
9: &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;$u_h^{n,k,0} \leftarrow u_h^{n,k-1,+}$
10: **end if**
11: **while** $e_l > \delta$ **do**
12: &nbsp;&nbsp;&nbsp;&nbsp;$\tilde{a}(u_h^{n,k,l-1}; u_h^{n,k,l}, w_h) = a(u_h^{n,k-1}, w_h)$
13: &nbsp;&nbsp;&nbsp;&nbsp;$+\omega((u_h^{n-1,+}, w_h) - a(u_h^{n,k-1,+}, w_h) - \tau s_h(u_h^{n,k-1,-}, w_h))$, $\forall w_h \in V_h$
14: &nbsp;&nbsp;&nbsp;&nbsp;$e_l \leftarrow \|u_h^{n,k,l} - u_h^{n,k,l-1}\|_{H^1}$
15: **end while**
16: $u_h^{n,k} \leftarrow u_h^{n,k,l}$
17: $e_k \leftarrow \|u_h^{n,k} - u_h^{n,k-1}\|_{H^1}$
18: **end while**
19: $u_h^n \leftarrow u_h^{n,k}$
20: **end for**
21: **return** $u_h^+$

365 To illustrate the numerical efficiency of our scheme (12), we compare it with another bound-preserving scheme derived using an operator-splitting technique as follows: At each time step $n \in \mathcal{N}$, find $u_h^n \in V_h$ such that

367

$$ \left( \frac{u_h^{n*} - u_h^{n-1}}{\tau}, v_h \right) + (\nabla u_h^{n*}, \nabla v_h) = 0, \quad \forall v_h \in V_h, \qquad (46a) $$

$$ U_i^n = \frac{U_i^{n*} \exp \left\{ \frac{\tau f(U_i^{n*})}{\epsilon^2 (U_i^{n*}+1)(1-U_i^{n*})} \right\}}{1 + \left( \exp \left\{ \frac{\tau f(U_i^{n*})}{\epsilon^2 (U_i^{n*}+1)(1-U_i^{n*})} \right\} - 1 \right)}, \quad \forall i \in \mathcal{V}, \qquad (46b) $$

368 where $u_h^{n*} := \sum_{i \in \mathcal{V}} U_i^{n*} \varphi_i$. The first step (i.e., the parabolic evolution) is bound-preserving if the mesh size is sufficiently small and satisfies the acute angle condition, and the second step (i.e., the stiff ODE evolution) is a linearized bound-preserving ODE solver proposed in [16, Section 3.2].

369

370

371

## 5.3 Numerical results

372 Under the tolerance $\delta := 10^{-12}$, the stabilization parameter $\alpha = 1$ and the relaxation parameter $\omega = 1$, our scheme converges with 3 Newton iterations and 2 Richardson-like iterations (so 6 iterations in total) at each time step, independently of $\epsilon$. In Tables 1 and 2, we present the relative errors in the $L^2$-norm and $H^1$-seminorm at the final time (i.e., $\frac{\|u(T)-u_h^+(T)\|}{\|u(T)\|}$ and $\frac{\|\nabla(u(T)-u_h^+(T))\|}{\|\nabla u(T)\|}$),

373

374

375
376

377

378

379

380

381

382

383

384

385

386

387

388

389

390

391

392

393

394

395

396

397

398

399

400

401

402

# 6 Conclusion

403

404

405

406

407

408

409

410

In this work, we developed a novel bound-preserving numerical scheme that inherits the energy decay property from the continuous level under reasonable time step constraints. Our analysis establishes an optimal error bound in the $L^2(H^1)$-seminorm with only polynomial $\epsilon^{-1}$-dependence. Numerical experiments comparing the proposed scheme with a bound-preserving operator-splitting approach showcased its superior accuracy and computational efficiency in under-resolved regimes. Future work will focus on two key extensions: (1) The development of higher-order generalizations of the current scheme; (2) Its application to a broader class of phase-field models.

# References

411

412

413

[1] G. Akrivis and B. Li. Error estimates for fully discrete BDF finite element approximations of the Allen-Cahn equation. *IMA J. Numer. Anal.*, 42(1):363–391, 2022.

414

415

[2] S. M. Allen and J. W. Cahn. A microscopic theory for antiphase boundary motion and its application to antiphase domain coarsening. *Acta Metallurgica*, 27(6):1085–1095, 1979.
<table><caption>Table 1: Relative errors and relative cpu times for Test case 1 with topological changes. (For each value of $\epsilon$, cpu times are normalized by the cpu time of our scheme on the coarsest mesh.)</caption><thead><tr><th colspan="4">ε = 0.05</th><th colspan="4">Our scheme</th></tr><tr><th>cells</th><th>L<sup>2</sup></th><th>H<sup>1</sup></th><th>cpu</th><th>L<sup>2</sup></th><th>H<sup>1</sup></th><th>cpu</th></tr></thead><tbody><tr><td>*40<sup>2</sup></td><td>2.25e-1</td><td>6.84e-1</td><td>0.22</td><td>4.25e-1</td><td>2.15e0</td><td>1.00</td></tr><tr><td>80<sup>2</sup></td><td>1.93e-1(0.23)</td><td>6.12e-1(0.16)</td><td>0.78</td><td>3.16e-1(0.43)</td><td>7.76e-1(1.49)</td><td>4.89</td></tr><tr><td>160<sup>2</sup></td><td>1.49e-1(0.37)</td><td>5.14e-1(0.25)</td><td>5.32</td><td>5.48e-1(-0.80)</td><td>1.23e0(-0.67)</td><td>39.43</td></tr><tr><td>320<sup>2</sup></td><td>1.00e-1(0.58)</td><td>4.00e-1(0.36)</td><td>38.41</td><td>2.18e-2(4.67)</td><td>1.05e-1(3.56)</td><td>267.35</td></tr></tbody><thead><tr><th colspan="4">ε = 0.1</th><th colspan="4">Our scheme</th></tr><tr><th>cells</th><th>L<sup>2</sup></th><th>H<sup>1</sup></th><th>cpu</th><th>L<sup>2</sup></th><th>H<sup>1</sup></th><th>cpu</th></tr></thead><tbody><tr><td>*40<sup>2</sup></td><td>1.52e-1</td><td>4.99e-1</td><td>0.10</td><td>3.26e-1</td><td>6.78e-1</td><td>1.00</td></tr><tr><td>80<sup>2</sup></td><td>1.14e-1(0.43)</td><td>4.08e-1(0.30)</td><td>0.41</td><td>4.43e-2(2.93)</td><td>1.72e-1(2.01)</td><td>2.42</td></tr><tr><td>160<sup>2</sup></td><td>7.61e-2(0.59)</td><td>2.85e-1(0.52)</td><td>3.01</td><td>2.05e-2(1.12)</td><td>8.61e-2(1.01)</td><td>21.72</td></tr><tr><td>**320<sup>2</sup></td><td>4.52e-2(0.76)</td><td>1.71e-1(0.73)</td><td>27.58</td><td>1.02e-2(1.01)</td><td>4.18e-2(1.05)</td><td>211.67</td></tr></tbody><thead><tr><th colspan="4">ε = 1</th><th colspan="4">Our scheme</th></tr><tr><th>cells</th><th>L<sup>2</sup></th><th>H<sup>1</sup></th><th>cpu</th><th>L<sup>2</sup></th><th>H<sup>1</sup></th><th>cpu</th></tr></thead><tbody><tr><td>**40<sup>2</sup></td><td>1.05e-1</td><td>5.95e-1</td><td>0.40</td><td>1.05e-1</td><td>5.93e-1</td><td>1.00</td></tr><tr><td>80<sup>2</sup></td><td>6.24e-2(0.77)</td><td>3.39e-1(0.83)</td><td>2.00</td><td>6.17e-2(0.78)</td><td>3.33e-1(0.84)</td><td>7.59</td></tr><tr><td>160<sup>2</sup></td><td>3.35e-2(0.91)</td><td>1.76e-1(0.95)</td><td>9.24</td><td>3.26e-2(0.93)</td><td>1.71e-1(0.97)</td><td>79.40</td></tr><tr><td>320<sup>2</sup></td><td>1.69e-2(0.99)</td><td>8.78e-2(1.01)</td><td>87.36</td><td>1.03e-2(1.02)</td><td>8.26e-2(1.05)</td><td>702.38</td></tr></tbody></table>

416

417

418

[3] A. Amiri, G. R. Barrenechea, and T. Pryer. A nodally bound-preserving finite element method for reaction-convection-diffusion equations. *Math. Models Methods Appl. Sci.*, 34(8):1533–1565, 2024.

419

420

[4] G. R. Barrenechea, E. Georgoulis, T. Pryer, and A. Veeser. A nodally bound-preserving finite element method. *IMA J. Numer. Anal.*, 44(4):2198–2219, 2024.

421

422

[5] G. R. Barrenechea, T. Pryer, and A. Trenam. A nodally bound-preserving discontinuous galerkin method for the drift-diffusion equation. *arXiv preprint arXiv:2410.05040*, 2024.

423

424

[6] S. Bartels. Robustness of error estimates for phase field models at a class of topological changes. *Comput. Methods Appl. Mech. Engrg.*, 288:75–82, 2015.

425

426

[7] S. Bartels. *Numerical Methods for Nonlinear Partial Differential Equations*, volume 47 of Springer Series in Computational Mathematics. Springer, Cham, 2015.

427

428

429

[8] S. Bartels and R. Müller. Quasi-optimal and robust a posteriori error estimates in $L^\infty(L^2)$ for the approximation of Allen-Cahn equations past singularities. *Math. Comp.*, 80(274):761–780, 2011.

430

431

432

[9] S. Bartels, R. Müller, and C. Ortner. Robust a priori and a posteriori error analysis for the approximation of Allen-Cahn and Ginzburg-Landau equations past topological changes. *SIAM J. Numer. Anal.*, 49(1):110–134, 2011.

433

434

435

[10] J. H. Brandts, S. Korotov, and M. Křížek. The discrete maximum principle for linear simplicial finite element approximations of a reaction-diffusion problem. *Linear Algebra Appl.*, 429(10): 2344–2357, 2008.
<table><thead><tr><th>ε = 0.05</th><th colspan="3">Operator-splitting</th><th colspan="3">Our scheme</th></tr><tr><th>cells</th><th>L<sup>2</sup></th><th>H<sup>1</sup></th><th>cpu time</th><th>L<sup>2</sup></th><th>H<sup>1</sup></th><th>cpu time</th></tr></thead><tbody><tr><td>*40<sup>2</sup></td><td>7.01e-2</td><td>5.78e-1</td><td>0.06</td><td>1.78e-1</td><td>9.73e-1</td><td>1.00</td></tr><tr><td>80<sup>2</sup></td><td>6.80e-2 (0.05)</td><td>5.64e-1(0.04)</td><td>0.36</td><td>7.66e-2(1.24)</td><td>6.06e-1(0.69)</td><td>2.28</td></tr><tr><td>160<sup>2</sup></td><td>5.93e-2 (0.20)</td><td>6.05e-1(0.16)</td><td>2.46</td><td>1.83e-2(2.08)</td><td>1.66e-1(1.88)</td><td>22.34</td></tr><tr><td>320<sup>2</sup></td><td>4.86e-2(0.29)</td><td>4.29e-1(0.24)</td><td>17.85</td><td>5.54e-3(1.72)</td><td>5.14e-2(1.70)</td><td>104.39</td></tr></tbody></table>

<table><thead><tr><th>ε = 0.1</th><th colspan="3">Operator-splitting</th><th colspan="3">Our scheme</th></tr><tr><th>cells</th><th>L<sup>2</sup></th><th>H<sup>1</sup></th><th>cpu time</th><th>L<sup>2</sup></th><th>H<sup>1</sup></th><th>cpu time</th></tr></thead><tbody><tr><td>*40<sup>2</sup></td><td>5.82e-2</td><td>4.85e-1</td><td>0.13</td><td>1.07e-1</td><td>7.64e-1</td><td>1.00</td></tr><tr><td>80<sup>2</sup></td><td>4.92e-2 (0.24)</td><td>3.99e-1(0.29)</td><td>0.59</td><td>1.52e-2(2.86)</td><td>1.25e-1(2.65)</td><td>2.67</td></tr><tr><td>160<sup>2</sup></td><td>3.54e-2 (0.48)</td><td>2.87e-1(0.48)</td><td>4.18</td><td>7.00e-3(1.13)</td><td>6.49e-2(0.96)</td><td>23.72</td></tr><tr><td>**320<sup>2</sup></td><td>2.13e-2(0.74)</td><td>1.74e-1(0.73)</td><td>30.20</td><td>2.71e-3(1.37)</td><td>2.75e-2(1.25)</td><td>201.58</td></tr></tbody></table>

<table><caption>Table 2: Relative errors and relative cpu times for Test case 2 without topological changes. (For each value of $\epsilon$, cpu times are normalized by the cpu time of our scheme on the coarsest mesh.)</caption><thead><tr><th>ε = 1</th><th colspan="3">Operator-splitting</th><th colspan="3">Our scheme</th></tr><tr><th>cells</th><th>L<sup>2</sup></th><th>H<sup>1</sup></th><th>cpu time</th><th>L<sup>2</sup></th><th>H<sup>1</sup></th><th>cpu time</th></tr></thead><tbody><tr><td>**40<sup>2</sup></td><td>2.14e-2</td><td>2.45e-1</td><td>0.44</td><td>2.10e-2</td><td>2.43e-1</td><td>1.00</td></tr><tr><td>80<sup>2</sup></td><td>1.31e-2 (0.72)</td><td>1.50e-1(0.73)</td><td>1.71</td><td>1.28e-2(0.72)</td><td>1.48e-1(0.73)</td><td>5.90</td></tr><tr><td>160<sup>2</sup></td><td>7.04e-3 (0.90)</td><td>7.97e-2(0.91)</td><td>11.78</td><td>6.89e-3(0.90)</td><td>7.87e-2(0.92)</td><td>51.67</td></tr><tr><td>320<sup>2</sup></td><td>3.16e-3(1.16)</td><td>3.54e-2(1.18)</td><td>85.12</td><td>3.07e-3(1.17)</td><td>3.48e-2(1.18)</td><td>457.78</td></tr></tbody></table>

436

437

[11] X. Chen. Spectrum for the Allen-Cahn, Cahn-Hilliard, and phase-field equations for generic interfaces. *Comm. Partial Differential Equations*, 19(7-8):1371-1395, 1994.

438

439

440

[12] K. Chrysafinos. Stability analysis and best approximation error estimates of discontinuous time-stepping schemes for the Allen-Cahn equation. *ESAIM Math. Model. Numer. Anal.*, 53 (2):551-583, 2019.

441

442

[13] P. de Mottoni and M. Schatzman. Geometrical evolution of developed interfaces. *Trans. Amer. Math. Soc.*, 347(5):1533-1589, 1995.

443

444

[14] Q. Du, L. Ju, X. Li, and Z. Qiao. Maximum bound principles for a class of semilinear parabolic equations and exponential time-differencing schemes. *SIAM Rev.*, 63(2):317-359, 2021.

445

446

[15] A. Ern and J.-L. Guermond. *Finite Elements III. First-Order and Time-Dependent PDEs*, volume 74 of *Texts in Applied Mathematics*. Springer, Cham, 2021.

447

448

449

[16] A. Ern, J.-L. Guermond, and Z. Wang. Asymptotic and invariant-domain preserving schemes for scalar conservation equations with stiff source terms and multiple equilibrium points. *J. Sci. Comput.*, 100(3):Paper No. 83, 30, 2024.

450

451

452

[17] X. Feng and Y. Li. Analysis of symmetric interior penalty discontinuous Galerkin methods for the Allen-Cahn equation and the mean curvature flow. *IMA J. Numer. Anal.*, 35(4): 1622-1651, 2015.

453

454

[18] X. Feng and A. Prohl. Numerical analysis of the Allen-Cahn equation and approximation for mean curvature flows. *Numer. Math.*, 94(1):33-65, 2003.

455

456

[19] X. Feng and A. Prohl. Analysis of a fully discrete finite element method for the phase field model and approximation of its sharp interface limits. *Math. Comp.*, 73(246):541-567, 2004.
457

458

[20] X. Feng and H.-J. Wu. A posteriori error estimates and an adaptive finite element method for the Allen-Cahn equation and the mean curvature flow. *J. Sci. Comput.*, 24(2):121-146, 2005.

459

460

[21] P. Grisvard. *Elliptic Problems in Nonsmooth Domains*, volume 69 of *Classics in Applied Mathematics*. SIAM, Philadelphia, PA, 2011.

461

462

[22] B. Li, J. Yang, and Z. Zhou. Arbitrarily high-order exponential cut-off methods for preserving maximum principle of parabolic equations. *SIAM J. Sci. Comput.*, 42(6):A3957–A3978, 2020.

463

464

465

[23] C. Liu, B. Riviere, J. Shen, and X. Zhang. A simple and efficient convex optimization based bound-preserving high order accurate limiter for cahn-hilliard-navier-stokes system. *SIAM J. Sci. Comput.*, 46(3):A1923–A1948, 2024.

466

467

[24] M. Renardy and R. C. Rogers. *An Introduction to Partial Differential Equations*, volume 13 of *Texts in Applied Mathematics*. Springer-Verlag, New York, second edition, 2004.

468

469

[25] J. Shen and X. Zhang. Discrete maximum principle of a high order finite difference scheme for a generalized Allen-Cahn equation. *Commun. Math. Sci.*, 20(5):1409-1436, 2022.

470

471

[26] V. Thomée. *Galerkin Finite Element Methods for Parabolic Problems*, volume 25 of Springer Series in Computational Mathematics. Springer-Verlag, Berlin, second edition, 2006.

472

473

[27] F. Verdugo and S. Badia. The software design of Gridap: a finite element package based on the Julia JIT compiler. *Comput. Phys. Commun.*, 276:Paper No. 108341, 24, 2022.

474

475

476

[28] J. Yang, Z. Yuan, and Z. Zhou. Arbitrarily high-order maximum bound preserving schemes with cut-off postprocessing for Allen-Cahn equations. *J. Sci. Comput.*, 90(2):Paper No. 76, 36, 2022.
<figure><img src="image_3.png" /><figcaption>Figure 1: Numerical illustration for Test case 1 at the final time with 100 × 2¹⁰ cells.</figcaption></figure>
<figure><img src="image_4.png" /><figcaption>Figure 2: Numerical illustration for Test case 2 at the final time with 100 × 2¹⁰ cells.</figcaption></figure>