Numerical analysis of a hybridized discontinuous Galerkin method for the Cahn-Hilliard problem

TL;DR

This paper presents a comprehensive numerical analysis of a hybridized discontinuous Galerkin method applied to the Cahn-Hilliard problem, establishing stability, convergence, and new discrete inequalities for the method.

Contribution

It introduces a stable, convergent hybridizable discontinuous Galerkin scheme for the Cahn-Hilliard equations with new discrete inequalities and error estimates.

Findings

Proved unconditional stability of the scheme.

Derived a priori error estimates for the method.

Established discrete Agmon and Gagliardo-Nirenberg inequalities.

Abstract

The mixed form of the Cahn-Hilliard equations is discretized by the hybridizable discontinuous Galerkin method. For any chemical energy density, existence and uniqueness of the numerical solution is obtained. The scheme is proved to be unconditionally stable. Convergence of the method is obtained by deriving a priori error estimates that are valid for the Ginzburg-Lindau chemical energy density and for convex domains. The paper also contains discrete functional tools, namely discrete Agmon and Gagliardo-Nirenberg inequalities, which are proved to be valid in the hybridizable discontinuous Galerkin spaces.

Tables

Table 1. Table 1: L 2 superscript 𝐿 2 L^{2} error ∥ c h j − c h 6 ∥ L 2 ​ ( Ω ) subscript delimited-∥∥ subscript 𝑐 subscript ℎ 𝑗 subscript 𝑐 subscript ℎ 6 superscript 𝐿 2 Ω \mathinner{\lVert c_{h_{j}}-c_{h_{6}}\rVert}_{L^{2}(\Omega)} between the discrete approximations for the order parameter at the j th superscript 𝑗 th j^{\text{th}} refinement ( j = 3 , 4 , 5 𝑗 3 4 5 j=3,4,5 ) and the solution on the fine grid ( j = 6 ) 𝑗 6 (j=6) for various fixed interface parameters κ 𝜅 \kappa computed at the final time T = 0.1 𝑇 0.1 T=0.1 and the corresponding estimated rates of convergence.

			$\kappa =1/2^8$		$\kappa =1/2^{10}$		$\kappa =1/2^{12}$
$j$	$h_j$	${\tau }_j$	Error	Rate	Error	Rate	Error	Rate
3	$1/2^3$	$0.1/2^6$	5.652$\times {10}^{-2}$	—	1.422$\times {10}^{-1}$	—	1.441$\times {10}^{-1}$	—
4	$1/2^4$	$0.1/2^8$	1.206$\times {10}^{-2}$	2.229	3.078$\times {10}^{-2}$	2.208	7.822$\times {10}^{-2}$	0.881
5	$1/2^5$	$0.1/2^{10}$	2.403$\times {10}^{-3}$	2.327	6.340$\times {10}^{-3}$	2.279	1.939$\times {10}^{-2}$	2.012

Equations

$${\partial }_{t}c-\Delta \mu$$

$$\mu$$

$$\nabla c\cdot \mathbf{n}$$

$$\nabla \mu \cdot \mathbf{n}$$

$$c$$

$$\Phi (c)={\Phi }_{+}(c)+{\Phi }_{-}(c),$$

$$S_h$$

$$M_h$$

$${\hat{S}}_h$$

$${\mathbf{S}}_h=S_h\times {\hat{S}}_h,{\mathbf{M}}_h=M_h\times {\hat{S}}_h.$$

$$H^s({\mathcal{E}}_h)=\left\{v\in L^2(\Omega )\text{mathchar}58\forall E\in {\mathcal{E}}_h,v{|}_E\in H^s(E)\right\}.$$

$$j_0({\mathbf{u}}_h,{\mathbf{v}}_h)$$

$$j_1({\mathbf{u}}_h,{\mathbf{v}}_h)$$

$$({\mathbf{u}}_h,{\mathbf{v}}_h{)}_{0,h}$$

$$({\mathbf{u}}_h,{\mathbf{v}}_h{)}_{1,h}$$

$${\Vert {\mathbf{v}}_h\Vert }_{0,h}$$

$${\Vert {\mathbf{v}}_h\Vert }_{1,h}$$

$$j_0(\mathbf{v},\mathbf{v})\le h^2\sum _{E\in {\mathcal{E}}_h}\frac{1}{h_E}{\Vert v-\hat{v}\Vert }_{L^2(\partial E)}^2\le h^2\Vert \mathbf{v}{\Vert }_{1,h}^2,\forall \mathbf{v}=(v,\hat{v})\in H^1({\mathcal{E}}_h)\times L^2({\Gamma }_h).$$

$${\Vert \mathbf{v}\Vert }_{0,h}^2={\Vert v\Vert }_{L^2(\Omega )}^2+j_0(\mathbf{v},\mathbf{v})$$

$$|j_0(\mathbf{u},\mathbf{v})|\le j_0(\mathbf{u},\mathbf{u}{)}^{1/2}{j}_0(\mathbf{v},\mathbf{v}{)}^{1/2}$$

$$\mathinner{\!\left\lVert\boldsymbol{v}\right\rVert}_{1,h,\star}=\mathinner{\biggl{(}\mathinner{\!\left\lVert\boldsymbol{v}\right\rVert}_{1,h}^{2}+\sum_{E\in\mathcal{E}_{h}}h_{E}\mathinner{\!\left\lVert\nabla v\cdot\mathbf{n}_{E}\right\rVert}_{L^{2}(\partial E)}^{2}\biggr{)}}^{1/2}.$$

$$\Vert {\mathbf{v}}_h{\Vert }_{1,h}\le {\Vert {\mathbf{v}}_h\Vert }_{1,h,\star }\lesssim \Vert {\mathbf{v}}_h{\Vert }_{1,h},\forall {\mathbf{v}}_h\in {\mathbf{S}}_h.$$

$$\mathinner{\!\left\lVert v\right\rVert}_{\mathrm{DG}}=\bigg{(}\sum_{E\in\mathcal{E}_{h}}\mathinner{\!\left\lVert\nabla v\right\rVert}_{L^{2}(\Omega)}^{2}+\sum_{e\in\Gamma_{h}^{0}}\frac{1}{h_{e}}\mathinner{\!\left\lVert[v]\right\rVert}_{L^{2}(e)}^{2}\bigg{)}^{1/2}.$$

$${\Vert v\Vert }_{\mathrm{DG}}\lesssim {\Vert \mathbf{v}\Vert }_{1,h},\forall \mathbf{v}\in H^1({\mathcal{E}}_h)\times L^2({\Gamma }_h).$$

$${\Vert [v]\Vert }_{L^2(e)}^2={\Vert v{|}_{E_{+}}-v{|}_{E_{-}}\Vert }_{L^2(e)}^2\le {\Vert v{|}_{E_{+}}-\hat{v}\Vert }_{L^2(e)}^2+{\Vert v{|}_{E_{-}}-\hat{v}\Vert }_{L^2(e)}^2.$$

$$\begin{array}{rl}\sum _{e\in {\Gamma }_h^0}\frac{1}{h_e}{\Vert [v]\Vert }_{L^2(e)}^2& \le \sum _{e\in {\Gamma }_h^0}\frac{1}{h_e}\left({\Vert v{|}_{E_{+}}-\hat{v}\Vert }_{L^2(e)}^2+{\Vert v{|}_{E_{-}}-\hat{v}\Vert }_{L^2(e)}^2\right)\\ & \le \sum _{E\in {\mathcal{E}}_h}\sum _{e\in {\mathcal{F}}_E\cap {\Gamma }_h^0}\frac{1}{h_e}\left({\Vert v{|}_{E_{+}}-\hat{v}\Vert }_{L^2(e)}^2+{\Vert v{|}_{E_{-}}-\hat{v}\Vert }_{L^2(e)}^2\right)\\ & \lesssim \sum _{E\in {\mathcal{E}}_h}\frac{1}{h_E}\left({\Vert v-\hat{v}\Vert }_{L^2(\partial E_{+})}^2+{\Vert v-\hat{v}\Vert }_{L^2(\partial E_{-})}^2\right)\\ & \lesssim \sum _{E\in {\mathcal{E}}_h}\frac{1}{h_E}{\Vert v-\hat{v}\Vert }_{L^2(\partial E)}^2,\end{array}$$

$$\overline{v}=\frac{1}{|\Omega |}{\int }_{\Omega }v,\forall v\in L^2(\Omega ).$$

$$\frac{1}{p^{\star }}=\frac{1}{2}-\frac{1}{d}.$$

$${\Vert v\Vert }_{L^p(\Omega )}\le C_P{\left({\Vert \mathbf{v}\Vert }_{1,h}^2+|\overline{v}{|}^2\right)}^{1/2},\forall \mathbf{v}=(v,\hat{v})\in H^1({\mathcal{E}}_h)\times L^2({\Gamma }_h).$$

$$\Vert v{\Vert }_{L^p(\Omega )}^2\lesssim \Vert v{\Vert }_{\mathrm{DG}}^2+|\Psi (v){|}^2,$$

$${\Vert \mathbf{v}\Vert }_{0,h}\le (C_P^2+h^2{)}^{1/2}\Vert \mathbf{v}{\Vert }_{1,h}\lesssim {\Vert \mathbf{v}\Vert }_{1,h},\forall \mathbf{v}\in (H^1({\mathcal{E}}_h)\cap L_0^2(\Omega ))\times L^2({\Gamma }_h).$$

$${\Vert v_h\Vert }_{L^q(\Omega )}$$

$${\Vert \nabla v_h\Vert }_{L^2(E)}$$

$$\Vert {\mathbf{v}}_h{\Vert }_{1,h}$$

$$\begin{gathered} a_{\mathcal{D}}(\mathbf{u},\mathbf{v})=\sum _{E\in {\mathcal{E}}_h}{\int }_E\nabla u\cdot \nabla vdx-\sum _{E\in {\mathcal{E}}_h}{\int }_{\partial E}\left(\nabla u\cdot {\mathbf{n}}_E\left(v-\hat{v}\right)+\left(u-\hat{u}\right)\nabla v\cdot {\mathbf{n}}_E\right)ds \\ +\sum _{E\in {\mathcal{E}}_h}\frac{\sigma }{h_E}{\int }_{\partial E}\left(u-\hat{u}\right)\left(v-\hat{v}\right)ds. \end{gathered}$$

Full text

1 Introduction

The Cahn–Hilliard equation was originally proposed in [1] as a model for phase separation in binary alloys. Since then, it has become fundamental to the phase field theory for moving interface problems. Some notable applications include tumor growth [2, 3] and multi-phase flows [4, 5, 6]. In its primal form, the Cahn–Hilliard equation is a fourth order nonlinear parabolic equation; thus, its numerical approximation presents a significant computational challenge. The conforming finite element approximation of fourth order elliptic operators is difficult, as the natural functional setting for weak solutions demands $H^{2}$ regularity. The construction and implementation of $H^{2}$ conforming elements are challenging tasks, especially in three dimensions. Alternatives to conforming methods for the Cahn–Hilliard problem in primal form are non-conforming approximations of $H^{2}$ such as the Morley element [7], or $\mathcal{C}^{0}$ interior-penalty methods [8, 9, 10]. The latter approach is particularly attractive for three dimensional simulations due to ease in which a basis can be constructed in higher dimensions.

Due to the difficulty involved with the numerical treatment of higher order derivatives, the mixed form of the Cahn–Hilliard equation is often preferred as it involves instead the solution of a coupled second order system. To our knowledge, this approach to the Cahn–Hilliard problem was first considered in [11] with classical $\mathcal{C}^{0}$ elements. However, it is well known that such methods violate the local mass balance satisfied at the continuous level. To this end, fully non-conforming discontinuous Galerkin (DG) methods are a suitable choice. For the use of DG methods in the numerical solution of the mixed Cahn–Hilliard system, we refer to [12, 13, 14] for the LDG method and [15, 16, 17] for the IPDG method.

Despite their advantages, DG methods come with a significant increase in the number of globally coupled degrees of freedom over conforming $\mathcal{C}^{0}$ elements. To address this additional computational burden, the hybridizable discontinuous Galerkin (HDG) methods were introduced in [18]. Key to the HDG methods is the introduction of additional unknowns on the mesh skeleton which act as Lagrange multipliers enforcing the continuity of the normal component of the numerical flux. As a consequence, the element unknowns can be eliminated locally through static condensation leading to a reduction in the total number of globally coupled degrees of freedom. The HDG method has seen success across a wide variety of elliptic and parabolic problems, and has recently been applied to the Cahn–Hilliard [2, 19] and Cahn–Hilliard–Navier–Stokes systems [4]. Closely related to the HDG method is the hybrid high order (HHO) method introduced in [20], which has been applied to the Cahn–Hilliard problem in [21, 22]. Similar to the HDG method, the HHO method introduces additional degrees of freedom on the mesh skeleton in order to leverage static condensation to reduce the size of the global system. The HHO method differs from the use of local reconstruction operators and face-based stabilizations. For more information on the ties between the HHO and HDG methods, we refer the reader to [23].

Regarding the theoretical analysis of HDG methods for the Cahn–Hilliard system, we mention that optimal error bounds for the hybridizable LDG method are proven in [19]. However, to the best of our knowledge, the theoretical analysis of the hybridized IPDG methods in [2, 4] is missing from the literature. We remark that the additional facet unknowns in our HDG scheme precludes the use of discrete functional analysis tools used previously in [15, 16] to analyze the IPDG method. Fortunately, appropriate analogues of these tools have been extended to the HHO setting in [21], and using similar techniques we will show that they hold also in the HDG setting.

The main contributions of this paper are (i) the unconditional unique solvability (Theorem 1) of HDG established by using the Minty–Browder and Brouwer fixed point theorems, (ii) the unconditional stability for any $\mathcal{C}^{2}$ potential function (Theorem 2), (iii) the $L^{\infty}$ stability of the order parameter (Theorem 3) established for convex domains and for the Ginzburg–Landau potential without any regularization, truncation, or extension, and (iv) optimal a priori error estimates in the broken $H^{1}$ norm (Theorem 4) for the Ginzburg–Landau potential and for convex domains. It is worth noting that a major challenge in the analysis is to avoid the use of any modifications or assumptions on the potential function such as the ones made in [17, 24]. Following similar strategies to [15, 25], we successfully avoid such assumptions by proving discrete Agmon and Gagliardo–Nirenberg inequalities in the HDG setting, which are useful stand-alone results that can be applied to other problems.

2 Model problem, notation and preliminaries

Let $\Omega\subset\mathbb{R}^{d}$, $d=2,3$ be a bounded, open, polygonal ($d=2$) or polyhedral ($d=3$) domain with outer unit normal $\boldsymbol{n}$. We consider the fourth order Cahn–Hilliard equation, rewritten as a second order system: find a pair $(c,\mu)$ satisfying

[TABLE]

We assume that the scalar potential function $\Phi$ admits a concave-convex decomposition; it suffices to assume $\Phi\in\mathcal{C}^{2}$. In other words, we can write

[TABLE]

with $\Phi_{+}$ convex and $\Phi_{-}$ concave.

2.1 Basic results on broken Sobolev and polynomial spaces

Let $\mathcal{E}_{h}$ be a conforming shape-regular mesh of $\Omega$, made of simplices $E$ with boundary $\partial E$ and diameter $h_{E}$. The mesh size is $h=\max_{E\in\mathcal{E}_{h}}h_{E}$. Let $\Gamma_{h}^{0}$ (resp $\Gamma_{h}^{b}$) be the set of interior (resp. boundary) faces and let $\Gamma_{h}=\Gamma_{h}^{0}\cup\Gamma_{h}^{b}$. Further, let $\mathcal{F}_{E}$ denote the set of all the faces of an element $E\in\mathcal{E}_{h}$. We assume that the family $\mathinner{\left\{\mathcal{E}_{h}\right\}}_{h>0}$ is quasi-uniform.

Let $k\geq 1$ be a fixed integer. We introduce a pair of broken polynomial spaces on $\mathcal{E}_{h}$:

[TABLE]

Moreover, the HDG method requires the following broken polynomial space defined on $\Gamma_{h}$:

[TABLE]

Here, $L_{0}^{2}(\Omega)$ is the zero mean value subspace of $L^{2}(\Omega)$ and $\mathbb{P}_{k}(\mathcal{O})$ denotes the space of polynomials of degree less than or equal to $k$ defined on the open set $\mathcal{O}$.

As the spaces $S_{h}$ and $M_{h}$ are non-conforming, we introduce the broken gradient operator $\nabla_{h}v_{h}$ by the restriction $(\nabla_{h}v_{h})|_{E}=\nabla(v_{h}|_{E})$. Moreover, the trace of a function $v_{h}\in S_{h}$ may be double-valued on interior facets. To each interior facet $e\in\Gamma_{h}^{0}$, we associate a unique normal vector $\mathbf{n}_{e}$ and denote by $E_{+}$ and $E_{-}$ the neighboring elements of $e$ such that $\mathbf{n}_{e}$ points from $E_{-}$ to $E_{+}$. We introduce the jump $[\cdot]$ and average $\{\cdot\}$ of $v_{h}\in S_{h}$ across an interior facet $e\in\Gamma_{h}^{0}$ as follows: let $[v_{h}]=v_{h}|_{E_{+}}-v_{h}|_{E_{-}}$ and $\{v_{h}\}=(v_{h}|_{E_{+}}+v_{h}|_{E_{-}})/2$. On boundary faces $e\in\Gamma_{h}^{b}$, we set $[v_{h}]=\{v_{h}\}=v_{h}|_{E}\text{ on }e$, where $E$ is the element such that $e\subset\partial E\cap\partial\Omega$.

We adopt the following notation for various product spaces of interest in this work:

[TABLE]

Pairs in these product spaces will be denoted using boldface; for example, $\boldsymbol{s}_{h}=(s_{h},\hat{s}_{h})\in\boldsymbol{S}_{h}$. Throughout we use the notation $a\lesssim b$ to denote $a\leq Cb$ where $C$ is a generic constant independent of the mesh parameters $h$ and $\tau$, but possibly dependent on the polynomial degree $k$, the spatial dimension $d$, and the domain $\Omega$.

Given an integer $s\geq 1$, we define the broken Sobolev space:

[TABLE]

Define $j_{0},j_{1}\mathrel{\mathop{\mathchar 58\relax}}\boldsymbol{S}_{h}\times\boldsymbol{S}_{h}\rightarrow\mathbb{R}$ as follows

[TABLE]

We equip $\boldsymbol{S}_{h}$ with the following (semi-)inner-products:

[TABLE]

as well as their induced (semi-)norms:

[TABLE]

The space $\boldsymbol{M}_{h}$ equipped with $(\cdot,\cdot)_{1,h}$ is an inner-product space. We note that $\mathinner{\!\left\lVert\cdot\right\rVert}_{1,h}$ is not a norm on $\boldsymbol{S}_{h}$. From the definition of $j_{0}(\cdot,\cdot)$ and the fact that $h_{E}\leq h$ , we have that

[TABLE]

Therefore we have for any $\boldsymbol{v}=(v,\hat{v})\in H^{1}(\mathcal{E}_{h})\times L^{2}(\Gamma_{h})$

[TABLE]

Moreover, on the broken Sobolev space $H^{2}(\mathcal{E}_{h})\times L^{2}(\Gamma_{h})$ we introduce another norm (with $\mathbf{n}_{E}$ denoting the unit outward normal vector to $\partial E$):

[TABLE]

It is equivalent to the norm $\|\cdot\|_{1,h}$ on $\boldsymbol{S}_{h}$ using trace and inverse inequalities:

[TABLE]

We recall the standard DG semi-norm for $H^{1}(\mathcal{E}_{h})$:

[TABLE]

The following bound holds:

[TABLE]

Indeed, the triangle inequality yields for any $\hat{v}\in L^{2}(\Gamma_{h}^{0})$

[TABLE]

We sum over all interior faces and use that $h_{e}\approx h_{E}$ which results from the shape and contact regularity of the mesh (see Lemma 1.41 and Lemma 1.42 in [26]). Thus, we obtain

[TABLE]

which is sufficient to conclude. Throughout the paper, we use the notation $\overline{v}$ to denote the mean value operator of any function $v$:

[TABLE]

With this notation, we recall Poincaré’s inequality valid on broken Sobolev spaces.

Lemma 1** (Poincaré inequality in $H^{1}(\mathcal{E}_{h})$)**

Let $p^{\star}$ be the exponent of the Sobolev embedding of $H^{1}(\Omega)$ into $L^{p}(\Omega)$ defined by

[TABLE]

Then, for each $p\leq p^{\star}$, there is a constant $C_{P}>0$ independent of $h$ such that

[TABLE]

Proof 1

From [27, 28], we have for $v\in H^{1}(\mathcal{E}_{h})$

[TABLE]

with $\Psi\in H^{1}(\mathcal{E}_{h})^{\prime}$ any bounded linear functional satisfying $\Psi(1)=1$. Choosing $\Psi\mathrel{\mathop{\mathchar 58\relax}}v\mapsto\overline{v}$, the result then follows with (20). \qed

As a consequence, with (23), we have

[TABLE]

From (24), it is evident that $\|\cdot\|_{1,h}$ defines a norm on $\boldsymbol{M}_{h}$. We now recall important inequalities which hold in the discrete spaces, thanks to the quasi-uniformity assumption on the mesh.

[TABLE]

The proofs of (25) and (26) can be found in [29]. To obtain (27), we use (26) and the definition of $j_{1}$.

We consider the symmetric discretization of the Laplace operator $-\Delta u$: $a_{\mathcal{D}}\mathrel{\mathop{\mathchar 58\relax}}(H^{2}(\mathcal{E}_{h}),L^{2}(\Gamma_{h}))\times(H^{2}(\mathcal{E}_{h}),L^{2}(\Gamma_{h}))\to\mathbb{R}$:

[TABLE]

The parameter $\sigma>0$ is a user-specified penalty parameter. We recall the following basic results concerning the bilinear form $a_{\mathcal{D}}$ [30, 31]:

Lemma 2** (Coercivity and continuity)**

Provided the penalty parameter $\sigma>0$ is chosen sufficiently large, the bilinear form is coercive on $\boldsymbol{S}_{h}$: there exists a constant $C_{\mathrm{coer}}>0$ such that

[TABLE]

Moreover, the bilinear form $a_{\mathcal{D}}$ is continuous on $\boldsymbol{S}_{h}\times\boldsymbol{S}_{h}$: there exists a constant $C_{\mathrm{cont}}>0$ such that

[TABLE]

and furthermore, there exists a constant $C_{b}^{\star}>0$ such that

[TABLE]

Henceforth, we will always assume $\sigma>0$ is sufficiently large to ensure the coercivity of the bilinear form $a_{\mathcal{D}}$.

3 Discrete functional analysis tools

The main goal of this section is to show discrete counter-parts to the Agmon (Lemma 6) and Gagliardo–Nirenberg (Lemma 7) inequalities in the HDG setting following ideas in [21, 32]. Such inequalities are important to establish stability estimates (see Theorem 3) which allow us to show convergence of the HDG scheme without any regularization of the potential function $\Phi$. We start with recalling and defining continuous and discrete Green and Laplace operators.

3.1 Continuous Green operator

We recall the continuous Green operator $G\mathrel{\mathop{\mathchar 58\relax}}V^{\prime}\rightarrow V$, where $V=H^{1}(\Omega)\cap L_{0}^{2}(\Omega)$. For any $f\in V^{\prime}$, $G(f)$ is the unique function in $V$ such that

[TABLE]

The embedding $V\subset L_{0}^{2}(\Omega)$ is dense and therefore, if $f\in L_{0}^{2}(\Omega)$ and $(\cdot,\cdot)_{\Omega}$ denote the $L^{2}$ inner-product, we find:

[TABLE]

In addition, if $\Omega$ is convex, we have

[TABLE]

3.2 Discrete Green operator

We introduce a discrete analogue of the Green operator, $G$, on the space $\boldsymbol{S}_{h}$. Consider $\boldsymbol{G}_{h}\mathrel{\mathop{\mathchar 58\relax}}\boldsymbol{S}_{h}\to\boldsymbol{M}_{h}$ satisfying:

[TABLE]

Note that the right-hand side of (36) defines, for fixed $\boldsymbol{w}_{h}\in\boldsymbol{S}_{h}$, a bounded linear functional on $\boldsymbol{M}_{h}$ by the Cauchy-Schwarz’s inequality and equivalence of norms on finite dimensional spaces. This fact, combined with the fact that $a_{\mathcal{D}}(\cdot,\cdot)$ is coercive on $\boldsymbol{M}_{h}$, shows that the operator $\boldsymbol{G}_{h}$ is well defined by the Lax-Milgram theorem.

3.3 Discrete Laplace operator

We introduce a discrete Laplace operator $\boldsymbol{\Delta}_{h}\mathrel{\mathop{\mathchar 58\relax}}\boldsymbol{S}_{h}\to\boldsymbol{M}_{h}$ as the unique solution to

[TABLE]

That $\boldsymbol{\Delta}_{h}$ is well-defined follows from the Riesz representation theorem, as the right-hand side defines a bounded linear functional on $\boldsymbol{M}_{h}$ while $(\cdot,\cdot)_{0,h}$ defines an inner-product on $\boldsymbol{S}_{h}$. We now show that $-\boldsymbol{G}_{h}$ is the inverse of the discrete Laplacian $\boldsymbol{\Delta}_{h}$ restricted to $\boldsymbol{M}_{h}$. To this end, the definitions of $\boldsymbol{G}_{h}$ (36) and $\boldsymbol{\Delta}_{h}$ (37) yield

[TABLE]

so that

[TABLE]

Choosing $\boldsymbol{v}_{h}=\boldsymbol{G}_{h}\boldsymbol{\Delta}_{h}\boldsymbol{w}_{h}+\boldsymbol{w}_{h}\in\boldsymbol{M}_{h}$, using the coercivity of the bilinear form $a_{\mathcal{D}}$ (29), and noting that $\mathinner{\!\left\lVert\cdot\right\rVert}_{1,h}$ is a norm on $\boldsymbol{M}_{h}$, we find that

[TABLE]

3.4 Properties of $\boldsymbol{\Delta}_{h}$ and $\boldsymbol{G}_{h}$

To set notation, for $\boldsymbol{w}_{h}\in\boldsymbol{S}_{h}$, we write

[TABLE]

In other words, $\Delta_{h}\boldsymbol{w}_{h}$ (or $G_{h}\boldsymbol{w}_{h}$) and $\hat{\Delta}_{h}\boldsymbol{w}_{h}$ (or $\hat{G}_{h}\boldsymbol{w}_{h}$) refer to, respectively, the element and face degrees of freedom obtained from $\boldsymbol{\Delta}_{h}\boldsymbol{w}_{h}$ (or $\boldsymbol{G}_{h}\boldsymbol{w}_{h}$).

We now proceed to show properties for $\boldsymbol{\Delta}_{h}$ and $\boldsymbol{G}_{h}$. To this end, we use the local $L^{2}$ projections. Let $E\in\mathcal{E}_{h}$, $e\in\Gamma_{h}$ and denote by $\pi_{h}$ and $\hat{\pi}_{h}$ the orthogonal $L^{2}$-projections satisfying

[TABLE]

Lemma 3** (Properties of $\pi_{h}$ and $\hat{\pi}_{h}$)**

Let $E\in\mathcal{E}_{h}$ and $1\leq p\leq\infty$. We have that

[TABLE]

Moreover, the following approximation results hold: let $s\in\mathbb{N}$ such that $1\leq s\leq k+1$. Then, it holds for all $v\in W^{s,p}(E)$ that

[TABLE]

where $e\in\mathcal{F}_{E}$. Moreover, if $v\in H^{s}(E)$, then it holds that

[TABLE]

Proof 2

Proofs of (43)-(46) can be found in [25]. The estimate (47) follows from (46) by observing that the best approximation property of the $L^{2}$-projection $\hat{\pi}_{h}$ ensures that

[TABLE]

\qed

With the above properties, we show the following inverse estimates for the discrete Laplacian.

Lemma 4

The following estimates hold:

[TABLE]

Proof 3

We begin by showing (49). Choosing $\boldsymbol{v}_{h}=\boldsymbol{\Delta}_{h}\boldsymbol{w}_{h}$ in (37) and using the continuity of the bilinear form $a_{\mathcal{D}}$ (30) and (27), we find

[TABLE]

This shows (49). Next, we show (50). As $\Delta_{h}\boldsymbol{w}_{h}\in M_{h}\subset L_{0}^{2}(\Omega)$, we have

[TABLE]

By the definition of the discrete Laplace operator (37) and the inner-product (9),

[TABLE]

*Observe that Lemma 3 implies that *

[TABLE]

Returning to (52), using (53), the continuity of the bilinear form $a_{\mathcal{D}}$ (30), (16), (54) and (49), we obtain

[TABLE]

Therefore, for any $\varphi\in H^{1}(\Omega)\cap L_{0}^{2}(\Omega)$ satisfying $\mathinner{\!\left\lVert\varphi\right\rVert}_{H^{1}(\Omega)}=1$, it holds that

[TABLE]

The result follows. \qed

We proceed to show approximation properties of the discrete Green operator.

Lemma 5

Assume that $\Omega$ is convex. Then, for any $\boldsymbol{w}_{h}\in\boldsymbol{M}_{h}$, the following estimates hold:

[TABLE]

where we define $\boldsymbol{\pi}_{h}(Gw_{h})=(\pi_{h}Gw_{h},\hat{\pi}_{h}Gw_{h})\in\boldsymbol{M}_{h}$.

Proof 4

By consistency and definition of $G(w_{h})$, we obtain that

[TABLE]

where $\boldsymbol{G}w_{h}=(Gw_{h},Gw_{h}|_{\Gamma_{h}})$. With (36), the following orthogonality relation easily follows:

[TABLE]

Thus, $\boldsymbol{G}_{h}\boldsymbol{w}_{h}$ is a modified HDG elliptic projection of $\boldsymbol{G}w_{h}$. Hence, the rest of the proof naturally modifies standard convergence proofs for HDG applied to elliptic problems. The details are provided in 0.A.\qed

3.5 Discrete Agmon inequality

Recall Agmon’s inequality (see e.g. [33, Lemma 4.10]):

[TABLE]

We now establish the following discrete counter-part in the HDG setting following closely the arguments in [21, 32].

Lemma 6** (Discrete Agmon Inequality)**

*Assume that $\Omega$ is convex. Then, the following inequality holds: *

[TABLE]

Proof 5

Considering (40) component-by-component, we observe that $w_{h}=-G_{h}\boldsymbol{\Delta}_{h}\boldsymbol{w}_{h}$. This fact, combined with the triangle inequality, yields (recall that $\boldsymbol{\Delta}_{h}\boldsymbol{w}_{h}=(\Delta_{h}\boldsymbol{w}_{h},\hat{\Delta}_{h}\boldsymbol{w}_{h})$)

[TABLE]

To bound $T_{1}$, we note that the $L^{\infty}$-stability of the $L^{2}$-projection operator (43) ensures that

[TABLE]

and thus the continuous Agmon inequality (60) yields (since $G(\Delta_{h}\boldsymbol{w}_{h})\in H^{2}(\Omega)$)

[TABLE]

By (34), (35), and (50), we find

[TABLE]

Next, we bound the term $T_{2}$. By the quasi-uniformity assumption and (25), we find that

[TABLE]

and thus by Lemma 5, (35), and (49),

[TABLE]

We conclude by noting that $d\leq 3$ and thus $h^{(3-d)/2}\lesssim 1$. \qed

Corollary 1

*Assume that $\Omega$ is convex. Then, the following inequality holds: *

[TABLE]

Proof 6

Let $\boldsymbol{w}_{h}\in\boldsymbol{S}_{h}$ and set $\overline{\boldsymbol{w}}_{h}=(\overline{w}_{h},\overline{w}_{h}|_{\Gamma_{h}})\in\boldsymbol{S}_{h}$. Then, $\boldsymbol{w}_{h}-\overline{\boldsymbol{w}}_{h}\in\boldsymbol{M}_{h}$, and thus by Lemma 6,

[TABLE]

Observe that $\mathinner{\!\left\lVert\overline{\boldsymbol{w}}_{h}\right\rVert}_{1,h}=0$, and

[TABLE]

and thus testing (70) with $\boldsymbol{v}_{h}=\boldsymbol{\Delta}_{h}\overline{\boldsymbol{w}}_{h}$,

[TABLE]

The result follows. \qed

3.6 Discrete Gagliardo–Nirenberg inequality

We recall the Gagliardo–Nirenberg inequality for bounded domains ([34], Theorem 1.2 in [35]). for $2\leq p\leq p^{\star}$ with $p^{\star}$ defined as in Lemma 1, it holds that

[TABLE]

where

[TABLE]

We show the following discrete counter-part.

Lemma 7** (Discrete Gagliardo–Nirenberg inequality)**

Suppose that $2\leq p\leq p^{\star}$, with $p^{\star}$ defined in Lemma 1 and let $\alpha$ be defined by (73). Then, it holds that

[TABLE]

As a consequence, we have

[TABLE]

Proof 7

Fix $\boldsymbol{w}_{h}\in\boldsymbol{M}_{h}$. To simplify notation, denote by $\boldsymbol{z}_{h}=(z_{h},\hat{z}_{h})=\boldsymbol{\Delta}_{h}\boldsymbol{w}_{h}$. Recall from (40) that $w_{h}=-G_{h}\boldsymbol{z}_{h}$. By the triangle inequality, we have that

[TABLE]

To bound the first term in the right-hand side of (76), we use the $W^{1,p}$-stability of the $L^{2}$-projection (44), the continuous Gagliardo–Nirenberg inequality (72), (34), (35), and (50):

[TABLE]

To bound the second term in the right-hand side of (76), we use (25), the definition of the norm $\mathinner{\!\left\lVert\cdot\right\rVert}_{1,h}$, and Lemma 5 since $\boldsymbol{z}_{h}\in\boldsymbol{M}_{h}$:

[TABLE]

Finally, by (49), we have

[TABLE]

We conclude by noting that $h^{d(\frac{1}{p}-\frac{1}{2})+\alpha}\lesssim 1$ since for $d\leq 3$ and $2\leq p\leq p^{\star}$ we have

[TABLE]

As in the proof of Corollary 1, we fix $\boldsymbol{w}_{h}\in\boldsymbol{S}_{h}$ and consider $\overline{\boldsymbol{w}}_{h}=(\overline{w}_{h},\overline{w}_{h}|_{\Gamma_{h}})$. Then $\boldsymbol{w}_{h}-\overline{\boldsymbol{w}}_{h}\in\boldsymbol{M}_{h}$ and we can apply (74). We then conclude by noting that $\|\overline{\boldsymbol{w}}_{h}\|_{1,h}=0$ and $\boldsymbol{\Delta}_{h}\overline{\boldsymbol{w}}_{h}=\mathbf{0}$. \qed

4 The numerical method

Let $0=t_{0}<t_{1}<\dots<t_{N}=T$ be a partition of the time interval $(0,T)$ into $N$ subintervals of equal length $\tau$. We discretize (1) in time using the first order implicit Euler method, and we employ the Eyre convex-concave splitting scheme [36] for the nonlinear chemical energy density. We treat the convex part implicitly and the concave part explicitly. The fully discrete numerical scheme reads: for any $1\leq n\leq N$, given $c_{h}^{n-1}\in S_{h}$, find $(\boldsymbol{c}_{h}^{n},\boldsymbol{\mu}_{h}^{n})\in\boldsymbol{S}_{h}\times\boldsymbol{S}_{h}$ satisfying

[TABLE]

Here, $\delta_{\tau}$ denotes the backward difference operator:

[TABLE]

We define $\boldsymbol{c}_{h}^{0}\in\boldsymbol{S}_{h}$ as the solution to the variational problem: given $c_{0}\in H^{2}(\Omega)\cap L_{0}^{2}(\Omega)$, set $\boldsymbol{c}_{0}=(c_{0},c_{0}|_{\Gamma_{h}})$ and find $\boldsymbol{c}_{h}^{0}\in\boldsymbol{S}_{h}$ such that

[TABLE]

To see that this variational problem is well-posed, first observe that any solution is unique. Indeed, if $\boldsymbol{w}_{1},\boldsymbol{w}_{2}\in\boldsymbol{S}_{h}$ solve (83), then $\boldsymbol{w}_{1}-\boldsymbol{w}_{2}\in\boldsymbol{M}_{h}$ and

[TABLE]

By testing (84) with $\boldsymbol{v}_{h}=\boldsymbol{w}_{1}-\boldsymbol{w}_{2}$ and using the coercivity of $a_{\mathcal{D}}(\cdot,\cdot)$ (29), we conclude that $\boldsymbol{w}_{1}=\boldsymbol{w}_{2}$. To show existence, it is enough to realize that if $\boldsymbol{w}_{h}\in\boldsymbol{M}_{h}$ is the unique solution of the auxiliary problem

[TABLE]

the existence of which is guaranteed by the Lax–Milgram theorem, then $\boldsymbol{w}_{h}+(\overline{c_{0}},\overline{c_{0}})$ solves (83), since for all $\boldsymbol{v}_{h}\in\boldsymbol{S}_{h}$,

[TABLE]

In the above, we used that for any $\boldsymbol{\theta}\in H^{2}(\mathcal{E}_{h})\times L^{2}(\Gamma_{h})$ and any constant $\alpha\in\mathbb{R}$, we have that

[TABLE]

Proposition 1** (Global mass conservation)**

The HDG scheme (81) satisfies the discrete global mass conservation

[TABLE]

Proof 8

Choose $\boldsymbol{\chi}_{h}=(1,1)$ in (81a) and use (86) to conclude that $(\delta_{\tau}c_{h}^{n},1)=0$ for any $n$, which yields the first equality. The second equality holds from (83). The third equality is a property of problem (1), which follows from integrating (1a) over $\Omega$, using the divergence theorem, and boundary condition (1d). \qed

The remaining part of this section focuses on showing that a unique solution exists for (81).

Theorem 1** (Unconditional unique solvability)**

Given $c_{h}^{n-1}\in S_{h}$, there exists a unique solution $(\boldsymbol{c}_{h}^{n},\boldsymbol{\mu}_{h}^{n})\in\boldsymbol{S}_{h}\times\boldsymbol{S}_{h}$ to (81) for any $\tau,h>0$.

Proof 9

The proof is divided into two steps. We first show that problem (81) is equivalent to problem (88) posed in $\boldsymbol{M}_{h}$, see Lemma 8. Then, we employ Brouwer’s fixed point theorem and the Minty-Browder theorem to show existence of solutions to (88) in Lemma 11 and Lemma 12 respectively. The result then immediately follows.\qed

4.1 An equivalent problem

Following the strategy in [6], we consider an equivalent problem to (81): for any $1\leq n\leq N$, find $(\boldsymbol{y}_{h}^{n},\boldsymbol{w}_{h}^{n})\in\boldsymbol{M}_{h}\times\boldsymbol{S}_{h}$

[TABLE]

Lemma 8

The unique solvability of the problem defined by (81) is equivalent to the unique solvability of the problem defined by (88).

Proof 10

Assume that $(\boldsymbol{c}_{h}^{n},\boldsymbol{\mu}_{h}^{n})\in\boldsymbol{S}_{h}\times\boldsymbol{S}_{h}$ is the unique solution to (81). We now exhibit a solution to (88). In particular, we show that $\boldsymbol{y}_{h}^{n}=(c_{h}^{n}-\overline{c_{0}},\hat{c}_{h}^{n}-\overline{c_{0}})$ and $\boldsymbol{w}_{h}^{n}=(\mu_{h}^{n}-\overline{\mu_{h}^{n}},\hat{\mu}_{h}^{n}-\overline{\mu_{h}^{n}})$ solve (88). Indeed, $\delta_{\tau}y_{h}^{n}=\delta_{\tau}c_{h}^{n}$ and $a_{\mathcal{D}}(\boldsymbol{w}_{h}^{n},\mathring{\boldsymbol{\chi}_{h}})=a_{\mathcal{D}}(\boldsymbol{\mu}_{h}^{n},\mathring{\boldsymbol{\chi}_{h}})$ from (86). Hence, it follows that (88a) holds. Noticing that $(\mu_{h}^{n}-\overline{\mu_{h}^{n}},\mathring{\phi}_{h})_{\Omega}=(\mu_{h}^{n},\mathring{\phi}_{h})_{\Omega}$, it also follows that (88b) holds.

Conversely, assume that $(\boldsymbol{y}_{h}^{n},\boldsymbol{w}_{h}^{n})\in\boldsymbol{M}_{h}\times\boldsymbol{M}_{h}$ solves (88). Let $f(y_{h}^{n})=\Phi_{+}^{\prime}(y_{h}^{n}+\overline{c_{0}})+\Phi_{-}^{\prime}(y_{h}^{n-1}+\overline{c}_{0})$. We show that $\boldsymbol{c}_{h}^{n}=(y_{h}^{n}+\overline{c_{0}},\hat{y}_{h}^{n}+\overline{c_{0}})$ and $\boldsymbol{\mu}_{h}^{n}=(w_{h}^{n}+\overline{f(y_{h}^{n})},\hat{w}_{h}^{n}+\overline{f(y_{h}^{n})})$ solve (81). To this end, let $\boldsymbol{\chi}_{h}\in\boldsymbol{S}_{h}$ and take $\mathring{\boldsymbol{\chi}}_{h}=(\chi_{h}-\overline{\chi_{h}},\hat{\chi}_{h}-\overline{\chi_{h}})\in\boldsymbol{M}_{h}$ in (88a). Since $\delta_{\tau}y_{h}^{n}\in M_{h}$, we have that $(\delta_{\tau}y_{h}^{n},\chi_{h}-\overline{\chi_{h}})=(\delta_{\tau}y_{h}^{n},\chi_{h})=(\delta_{\tau}c_{h}^{n},\chi_{h})$. In addition, from (86), $a_{\mathcal{D}}(\boldsymbol{w}_{h}^{n},\mathring{\boldsymbol{\chi}}_{h})=a_{\mathcal{D}}(\boldsymbol{w}_{h}^{n},\boldsymbol{\chi}_{h})=a_{\mathcal{D}}(\boldsymbol{\mu}_{h}^{n},\boldsymbol{\chi}_{h})$. Thus, (81a) holds. To show that (81b) holds, let $\boldsymbol{\phi}_{h}\in\boldsymbol{S}_{h}$ and take $\mathring{\boldsymbol{\phi}}_{h}=(\phi_{h}-\overline{\phi_{h}},\hat{\phi}_{h}-\overline{\phi_{h}})\in\boldsymbol{M}_{h}$ in (88b). Observe that from (86), $a_{\mathcal{D}}(\boldsymbol{y}_{h}^{n},\mathring{\boldsymbol{\phi}}_{h})=a_{\mathcal{D}}(\boldsymbol{c}_{h}^{n},\mathring{\boldsymbol{\phi}}_{h})=a_{\mathcal{D}}(\boldsymbol{c}_{h}^{n},\boldsymbol{\phi}_{h})$. Further, since $w_{h}^{n}\in M_{h}$, we obtain

[TABLE]

This shows that (88) is satisfied. Since we have shown how to obtain a solution to (81) from (88) and vice-versa, the unique solvability of each problem follows from the unique solvability of the other problem via a standard proof by contradiction. \qed

4.2 Existence and uniqueness of the discrete solution to

(88)

To prove the existence and uniqueness of the discrete solution, we will rely on the following Lemma 9 and Lemma 10 [37].

Lemma 9** (Corollary to Brouwer’s fixed point theorem)**

Let $(V,(\cdot,\cdot)_{V})$ be a finite-dimensional Hilbert space and let $f\mathrel{\mathop{\mathchar 58\relax}}V\to V$ be a continuous mapping with the following property: there exists an $M>0$ such that

[TABLE]

Then, there exists $v_{0}\in V$ such that

[TABLE]

Lemma 10** (Minty–Browder theorem)**

Let $(V,\mathinner{\!\left\lVert\cdot\right\rVert}_{V})$ be a separable and reflexive Banach space and let $A\mathrel{\mathop{\mathchar 58\relax}}V\to V^{\prime}$ be a coercive and hemicontinuous monotone operator. Then, given any $f\in V^{\prime}$, there exists a $u\in V$ such that

[TABLE]

Moreover, if $A$ is strictly monotone, then the solution $u$ is unique.

We now employ Lemma 9 to show that given $w_{h}^{n}\in M_{h}$, there exists a unique $\boldsymbol{y}_{h}^{n}\in\boldsymbol{M}_{h}$ solving (88b).

Lemma 11

Given $w_{h}^{n},y_{h}^{n-1}\in M_{h}$, there exists a unique $\boldsymbol{y}_{h}^{n}\in\boldsymbol{M}_{h}$ such that

[TABLE]

Proof 11

For fixed $\boldsymbol{v}_{h}\in\boldsymbol{M}_{h}$, we define a mapping $\boldsymbol{\mathcal{F}}\mathrel{\mathop{\mathchar 58\relax}}\boldsymbol{M}_{h}\to\boldsymbol{M}_{h}$ as follows: for all $\mathring{\boldsymbol{\phi}_{h}}\in\boldsymbol{M}_{h}$,

[TABLE]

The mapping $\boldsymbol{\mathcal{F}}$ is well-defined thanks to the Riesz representation theorem. Choosing $\mathring{\boldsymbol{\phi}}_{h}=\boldsymbol{v}_{h}$ in (91) and using the coercivity of the bilinear form $a_{\mathcal{D}}$ (29), we find

[TABLE]

Expanding $\Phi_{+}^{\prime}(v_{h}+\overline{c_{0}})$ around $\overline{c_{0}}$, we find there exists $\xi_{h}$ between $\overline{c_{0}}$ and $v_{h}+\overline{c_{0}}$ such that

[TABLE]

and thus by the convexity of $\Phi_{+}$ we have

[TABLE]

Next, the Cauchy-Schwarz’s inequality, discrete Poincaré’s inequality (23), and Young’s inequality yield

[TABLE]

so that, upon choosing $\epsilon=C_{P}^{2}/(C_{\mathrm{coer}}\kappa)$, we find

[TABLE]

Define the sphere $\Xi$ in $\boldsymbol{M}_{h}$ as follows:

[TABLE]

It holds that $\mathinner{\left(\boldsymbol{\mathcal{F}}(\boldsymbol{v}_{h}),\boldsymbol{v}_{h}\right)}_{1,h}\geq 0$ for any $\boldsymbol{v}_{h}\in\Xi$. Consequently, Brouwer’s fixed point theorem guarantees the existence of a function $\boldsymbol{y}_{h}^{n}\in\boldsymbol{M}_{h}$ such that $\boldsymbol{\mathcal{F}}(\boldsymbol{y}_{h}^{n})=0$ and

[TABLE]

Equivalently, $\boldsymbol{y}_{h}^{n}\in\boldsymbol{M}_{h}$ satisfies (90).

Lastly, we show that the solution $\boldsymbol{y}_{h}^{n}\in\boldsymbol{M}_{h}$ solving (90) is unique. To this end, assume $\boldsymbol{y}_{1}\in\boldsymbol{M}_{h}$ and $\boldsymbol{y}_{2}\in\boldsymbol{M}_{h}$ are two solutions of (90). In both cases, we test (90) with $\mathring{\boldsymbol{\phi}}_{h}=\boldsymbol{y}_{1}-\boldsymbol{y}_{2}$, subtract the two resulting equations, and use the coercivity of the bilinear form $a_{\mathcal{D}}$ to find

[TABLE]

As $\Phi_{+}\in\mathcal{C}^{2}$ is convex, $\Phi_{+}^{\prime}$ is nondecreasing, and thus

[TABLE]

Hence $\mathinner{\!\left\lVert\boldsymbol{y}_{1}-\boldsymbol{y}_{2}\right\rVert}_{1,h}^{2}\leq 0$, so that $\boldsymbol{y}_{1}=\boldsymbol{y}_{2}$ as required. \qed

We are now ready to show that (88) is well posed by employing Lemma 10 and Lemma 11.

Lemma 12

The scheme (88) is uniquely solvable for any fixed $\tau$, $h$, and $\kappa$.

Proof 12

For any $\boldsymbol{w}_{h}=(w_{h},\hat{w}_{h})\in\boldsymbol{M}_{h}$, let $\boldsymbol{y}_{h}^{n}(w_{h})\in\boldsymbol{M}_{h}$ be the unique solution to (90). There exists a well-defined operator $\boldsymbol{\mathcal{G}}\mathrel{\mathop{\mathchar 58\relax}}\boldsymbol{M}_{h}\to\boldsymbol{M}_{h}^{\prime}$ satisfying

[TABLE]

We begin by showing $\boldsymbol{\mathcal{G}}$ is coercive. Take $\mathring{\boldsymbol{\chi}}_{h}=\boldsymbol{w}_{h}$ in (101) and note that by the coercivity of the bilinear form $a_{\mathcal{D}}(\cdot,\cdot)$ (29), we have that

[TABLE]

The Cauchy–Schwarz’s, Poincaré’s (23), and Young’s inequalities then yield for any $\epsilon_{1}>0$:

[TABLE]

We now aim to find a lower bound on $(y_{h}^{n},w_{h})_{\Omega}$. To this end, we test (88b) with $\mathring{\boldsymbol{\phi}}_{h}=\boldsymbol{y}_{h}^{n}$ and once again use the coercivity of the bilinear form $a_{\mathcal{D}}$ (29):

[TABLE]

Proceeding as in Lemma 11, see the derivation of (94) (95), we find that

[TABLE]

and for any $\epsilon_{2}>0$ that

[TABLE]

Thus, we obtain

[TABLE]

Combining (103) and (107), we have

[TABLE]

Choosing $\epsilon_{1}=C_{P}^{2}/(C_{\mathrm{coer}}\tau),\epsilon_{2}=C_{P}^{2}/(2C_{\mathrm{coer}}\kappa)$ we find

[TABLE]

and thus

[TABLE]

Next, we show that the operator $\boldsymbol{\mathcal{G}}$ is bounded. For any test function $\mathring{\boldsymbol{\chi}}_{h}\in\boldsymbol{M}_{h}$, the Cauchy–Schwarz’s inequality, Poincaré inequality’s (23), and continuity of the bilinear form $a_{\mathcal{D}}$ (30) yield

[TABLE]

As $\boldsymbol{y}_{h}^{n}(w_{h})\in\boldsymbol{M}_{h}$ is the unique solution of (90), we take $\mathring{\boldsymbol{\phi}}_{h}=\boldsymbol{y}_{h}^{n}$ in (90) and use the coercivity of $a_{\mathcal{D}}$ to find

[TABLE]

Recall from (94) and (106) that there exists some $\epsilon_{3}>0$ such that

[TABLE]

where we have used the Cauchy–Schwarz’s, broken Poincaré’s, and Young’s inequalities to obtain the second inequality. Choosing $\epsilon_{3}=(2C_{P}^{2})/(C_{\mathrm{coer}}\kappa)$ and rearranging yields

[TABLE]

and using the fact that $\sqrt{a+b}\leq\sqrt{a}+\sqrt{b}$ for $a,b\geq 0$, we find

[TABLE]

Combining (111) and (115),

[TABLE]

Consequently, it holds that

[TABLE]

In other words, the operator $\boldsymbol{\mathcal{G}}$ maps bounded sets in $\boldsymbol{M}_{h}$ to bounded sets in $\boldsymbol{M}_{h}^{\prime}$.

Next, we show hemicontinuity of the operator $\boldsymbol{\mathcal{G}}$. It suffices to show for any $\boldsymbol{w}_{h},\boldsymbol{v}_{h},\mathring{\boldsymbol{\chi}}_{h}\in\boldsymbol{M}_{h}$ the (sequential) continuity of the mapping $g\mathrel{\mathop{\mathchar 58\relax}}\mathbb{R}\rightarrow\mathbb{R}$ defined by

[TABLE]

To this end, take an arbitrary $t^{\star}\in\mathbb{R}$ and a sequence $(t_{i})_{i\in\mathbb{N}}\subset\mathbb{R}$ such that $t_{i}\to t^{\star}$ as $i\to\infty$. From the definition of $\boldsymbol{\mathcal{G}}$ (101) and the bilinearity and continuity of the form $a_{\mathcal{D}}$ (30), we find that

[TABLE]

*Given $w_{h}+t^{\star}v_{h}\in M_{h}$ and $w_{h}+t_{i}v_{h}\in M_{h}$, we have respectively the solutions $\boldsymbol{y}_{h}^{n}(w_{h}+t^{\star}v_{h})$ and $\boldsymbol{y}_{h}^{n}(w_{h}+t_{i}v_{h})$ to (90). We can choose $\mathring{\boldsymbol{\phi}}_{h}=\boldsymbol{y}_{h}^{n}(w_{h}+t^{\star}v_{h})-\boldsymbol{y}_{h}^{n}(w_{h}+t_{i}v_{h})$ in (90), subtract the two resulting equations, and use the bilinearity and coercivity of $a_{\mathcal{D}}$ (29) to find: *

[TABLE]

For the first term, we use (100) to have

[TABLE]

Thus, with Cauchy–Schwarz’s inequality and the discrete Poincaré’s inequality (23), we obtain:

[TABLE]

Using this bound in (119), we obtain

[TABLE]

Since the right-hand side above tends to zero as $i$ tends to infinity, this shows that $\boldsymbol{\mathcal{G}}$ is hemicontinuous. Finally, we show the strict monotonicity of $\boldsymbol{\mathcal{G}}$. For any $\boldsymbol{w}_{h}=(w_{h},\hat{w}_{h}),\boldsymbol{s}_{h}=(s_{h},\hat{s}_{h})\in\boldsymbol{M}_{h}$, it holds that

[TABLE]

where we have used the coercivity of the bilinear form $a_{\mathcal{D}}$ (29). To each given $w_{h},s_{h}\in M_{h}$, we have corresponding unique solutions $\boldsymbol{y}_{h}^{n}(w_{h})$ and $\boldsymbol{y}_{h}^{n}(s_{h})$ to (90). Therefore, by testing (90) with $\mathring{\boldsymbol{\phi}}_{h}=\boldsymbol{y}_{h}^{n}(w_{h})-\boldsymbol{y}_{h}^{n}(s_{h})\in\boldsymbol{M}_{h}$ for each solution and subtracting the two resulting equations yields

[TABLE]

*where we have used the coercivity of $a_{\mathcal{D}}$ (29) and (100). Consequently, *

[TABLE]

and thus $\boldsymbol{\mathcal{G}}_{h}$ is monotone. Moreover, as $\mathinner{\lVert\cdot\rVert}_{1,h}$ defines a norm on $\boldsymbol{M}_{h}$, it immediately follows that the inequality in (122) is strict if $\boldsymbol{w}_{h}\neq\boldsymbol{s}_{h}$; hence, $\boldsymbol{\mathcal{G}}_{h}$ is strictly monotone.

We have now verified all of the assumptions of the Minty–Browder theorem, and we conclude that there exists a unique solution $\boldsymbol{w}_{h}^{n}$ satisfying

[TABLE]

implying that $(\boldsymbol{y}_{h}^{n}(w_{h}^{n}),\boldsymbol{w}_{h}^{n})\in\boldsymbol{M}_{h}\times\boldsymbol{M}_{h}$ is the unique solution to (88). \qed

5 Stability analysis for $\mathcal{C}^{2}$ chemical energy density

In this section, we show stability estimates for the discrete solutions. We will make use of the following operator. We define an operator $\boldsymbol{\mathcal{J}}_{h}\mathrel{\mathop{\mathchar 58\relax}}M_{h}\to\boldsymbol{M}_{h}$ such that

[TABLE]

which is well-defined by the Lax–Milgram theorem.

Lemma 13

For all $w_{h}\in M_{h}$ it holds that

[TABLE]

Further, it holds that

[TABLE]

Proof 13

Fix $w_{h}\in M_{h}$ and $\boldsymbol{v}_{h}\in\boldsymbol{S}_{h}$. Set $\overline{\boldsymbol{v}_{h}}=(\overline{v_{h}},\overline{v_{h}}|_{\Gamma_{h}})\in\boldsymbol{S}_{h}$. With (86), we have

[TABLE]

The proof of (126) is in 0.B. \qed

Theorem 2** (Unconditional energy stability)**

The following bound holds for all $1\leq m\leq N$:

[TABLE]

Proof 14

Choosing $\boldsymbol{\chi}_{h}=\tau\boldsymbol{\mu}^{n}_{h}$ in (81a), $\boldsymbol{\phi}_{h}=\boldsymbol{c}_{h}^{n}-\boldsymbol{c}_{h}^{n-1}$ in (81b), and subtracting the two resulting equations, we find

[TABLE]

*With a Taylor expansion of $\Phi_{+}$ and $\Phi_{-}$, there exist $\xi_{h}$ and $\zeta_{h}$ between $c_{h}^{n-1}$ and $c_{h}^{n}$ such that *

[TABLE]

Adding (130) and (131) and integrating over $\Omega$,

[TABLE]

As $\Phi_{+}$ and $-\Phi_{-}$ are convex, the second and third terms on the right-hand side of (132) are non-negative, and hence

[TABLE]

Further, by the symmetry, bilinearity, and coercivity of $a_{\mathcal{D}}$,

[TABLE]

Returning to (129) and using (133), (134), and the coercivity of the bilinear form $a_{\mathcal{D}}$ (29), we find

[TABLE]

Summing from $n=1$ to $n=m$ with $m\leq N$, we have

[TABLE]

By the coercivity (29) and continuity (30) of the bilinear form $a_{\mathcal{D}}$,

[TABLE]

which is the desired bound, save for one term.

To bound the remaining term, take $\boldsymbol{\chi}_{h}=\boldsymbol{\mathcal{J}}_{h}(\delta_{\tau}\boldsymbol{c}_{h}^{n})$ in (81a), and use the definition of $\boldsymbol{\mathcal{J}}_{h}$, (29), and (30). We have, on the one hand,

[TABLE]

and on the other hand,

[TABLE]

Thus, we obtain

[TABLE]

Thus, with (137), we obtain the remaining bound. This completes the proof. \qed

6 Analysis for the Ginzburg–Landau potential

For the remainder of the article, it will be assumed that the chemical energy density $\Phi$ is the Ginzburg–Landau potential:

[TABLE]

6.1 Uniform a priori bounds

Theorem 3

Assume that $\nabla c_{0}\cdot\boldsymbol{n}=0$ on the boundary $\partial\Omega$. There exists a constant $C$ independent of $h$ and $\tau$ but depending linearly on the final time $T$ such that the following holds.

[TABLE]

In addition, if $\Omega$ is convex, there is a constant $C$ independent of $h$ and $\tau$ such that

[TABLE]

Proof 15

To begin, we define $\boldsymbol{c}_{h}^{-1}=\boldsymbol{c}_{h}^{0}$ and $\boldsymbol{\mu}_{h}^{0}\in\boldsymbol{S}_{h}$ as the solution to the following variational problem

[TABLE]

Recall that $\boldsymbol{c}_{h}^{0}\in\boldsymbol{S}_{h}$ is defined by (83).

Therefore by testing (83) with $\boldsymbol{v}_{h}=\boldsymbol{c}_{h}^{0}$, we obtain from the coercivity and continuity of the bilinear form $a_{\mathcal{D}}$ that

[TABLE]

Choosing $\boldsymbol{\phi}_{h}=\boldsymbol{\mu}_{h}^{0}$ in (144) and using (83), we find

[TABLE]

and thus by the Cauchy–Schwarz’s inequality,

[TABLE]

where we have used that $c_{0}$ is a given quantity and that by the definition of the chemical energy density $\Phi$ (141) and the discrete Poincaré inequality Lemma 1,

[TABLE]

We note that since $\overline{c_{0}}=\overline{c_{h}^{0}}$:

[TABLE]

thanks to (145). Further, by Lemma 1, we have

[TABLE]

Next, we show that

[TABLE]

For $1\leq n\leq N$, we subtract (81b) with $n$ replaced by $n-1$ from (81b) and test with $\boldsymbol{\phi}_{h}=\boldsymbol{\mu}_{h}^{n}$ to find

[TABLE]

Testing (81a) with $\boldsymbol{\chi}_{h}^{n}=\kappa\tau\delta_{\tau}\boldsymbol{c}_{h}^{n}$, using the symmetry of $a_{\mathcal{D}}$, and combining with (150), we find

[TABLE]

Using the identity $2x(x-y)=x^{2}+(x-y)^{2}-y^{2}$, we find

[TABLE]

Moreover, it holds that

[TABLE]

Using (126), (20), and Young’s inequality, we have that

[TABLE]

*Thus, by Hölder’s, Cauchy–Schwarz’s, and Young’s inequalities, we find *

[TABLE]

By the triangle inequality,

[TABLE]

But for each $n\geq 1$ and for $2\leq p\leq p^{\star}$, we have from (128) and the fact that $\overline{c_{h}^{n}}=\overline{c_{0}}$,

[TABLE]

It is easy to show (with for instance Lemma 1 and (145)) that

[TABLE]

We conclude that for $2\leq p\leq p^{\star}$,

[TABLE]

This is true for all $n\geq 0$, thus we have

[TABLE]

With Young’s inequality, (155) becomes

[TABLE]

Using (1), we find

[TABLE]

Summing (158) from $n=1$ to $n=m$ (recall that we have defined $\boldsymbol{c}_{h}^{-1}$ such that $\delta_{\tau}\boldsymbol{c}_{h}^{0}=0$), we obtain

[TABLE]

Using (128) and (147), this bound implies

[TABLE]

We can conclude by using a discrete Grönwall inequality, but this will yield a constant that depends exponentially in time. Instead we directly obtain a bound on $\|\mu_{h}^{n}\|_{L^{2}(\Omega)}$. Testing (81b) with $\boldsymbol{\phi}_{h}=\boldsymbol{\mu}_{h}^{n}$, we have

[TABLE]

and thus using the Cauchy–Schwarz’s inequality, the continuity of the bilinear form $a_{\mathcal{D}}$, Young’s inequality, (156) and (128), we have

[TABLE]

and thus multiplying both sides of (162) by $\tau$, summing from $1\leq n\leq N$, and using (128), we find

[TABLE]

with $C$ depending linearly on the final time $T$. Returning to (160), we find for all $1\leq m\leq N$,

[TABLE]

with $C$ depending linearly on the final time $T$. Finally, we prove that

[TABLE]

For readibility, let $\boldsymbol{\chi}_{h}=\boldsymbol{\Delta}_{h}\boldsymbol{c}_{h}^{n}$. Testing the definition of the discrete Laplacian (37) with $\boldsymbol{v}_{h}=\boldsymbol{\chi}_{h}=(\chi_{h},\hat{\chi}_{h})$, and using (81b) and the definitions of $\Phi_{+}$ and $\Phi_{-}$ in (141), we find that

[TABLE]

The Cauchy–Schwarz’s inequality and Young’s inequality then yield

[TABLE]

Applying (156) and (164), we obtain

[TABLE]

Finally, by the discrete Agmon inequality Corollary 1, we have if the domain is convex

[TABLE]

and the right-hand side of (169) is uniformly bounded thanks to (128) and (168). \qed

6.2 Error analysis

The main goal of this section is to prove the following convergence result.

Theorem 4** (Error estimate)**

Let $2\leq s\leq k+1$. Assume that the weak solution has the following regularity:

[TABLE]

Further, assume the domain $\Omega\subset\mathbb{R}^{d}$, $d\in\mathinner{\left\{2,3\right\}}$ is convex. There exists $\tau_{0}>0$ such that for any $\tau\leq\tau_{0}$, the following error estimate holds. For any $1\leq m\leq N$,

[TABLE]

The proof of this estimate requires several intermediate results; thus, it is presented in subsection 6.4. We begin the analysis by introducing useful operators and interpolants.

6.3 Interpolation and intermediate results

Recall the definitions of the $L^{2}$ projections (42) and their approximation properties Lemma 3. For a given $v\in H^{s}(\Omega)$, set $\boldsymbol{v}=(v,v|_{\Gamma_{h}})$ and denote by $\boldsymbol{\pi}_{h}v=(\pi_{h}v,\hat{\pi}_{h}v)$, defined by (42). From  Lemma 3, it follows that

[TABLE]

Further, we consider the elliptic projection $\boldsymbol{\Pi}_{h}\mathrel{\mathop{\mathchar 58\relax}}H^{2}(\Omega)\rightarrow\boldsymbol{S}_{h}$ defined as follows. Fix $w\in H^{2}(\Omega)$ and set $\boldsymbol{w}=(w,w|_{\Gamma_{h}})$. The function $\boldsymbol{\Pi}_{h}w=(\Pi_{h}w,\hat{\Pi}_{h}w)$ is the unique function in $\boldsymbol{S}_{h}$ satisfying:

[TABLE]

Lemma 14** (Approximation properties of $\boldsymbol{\Pi}_{h}$)**

The following approximation properties hold. Fix $s$ such that $2\leq s\leq k+1$ and let $w\in H^{s}(\Omega)$. Define $\boldsymbol{w}=(w,\hat{w}|_{\Gamma_{h}})$.

[TABLE]

Proof 16

Observe by the definition of the $L^{2}$-projection $\pi_{h}$ (42), it holds that

[TABLE]

so that $\pi_{h}w-\Pi_{h}w\in M_{h}$. Consequently, by (18), the coercivity and continuity of $a_{\mathcal{D}}(\cdot,\cdot)$ (29) and (32), we have

[TABLE]

and therefore by the triangle inequality and the continuity of $a_{\mathcal{D}}(\cdot,\cdot)$ (31),

[TABLE]

The approximation properties of the projection $\boldsymbol{\pi}_{h}$ Lemma 3 yield

[TABLE]

Since the domain $\Omega$ is convex, we can prove optimal $L^{2}$-estimates for $\Pi_{h}$ using a standard duality argument. To this end, consider the auxiliary Neumann problem: find $z\in H^{2}(\Omega)\cap L_{0}^{2}(\Omega)$ such that

[TABLE]

We find that, on the one hand,

[TABLE]

while on the other hand, the single-valuedness of $z$ across element interfaces yields,

[TABLE]

where $\boldsymbol{z}=(z,z|_{\Gamma_{h}})$. Using the regularity of $w$ and $z$ and the fact that $\nabla z\cdot n=0$ on $\partial\Omega$, we find

[TABLE]

while an element-by-element integration by parts shows that

[TABLE]

Therefore,

[TABLE]

Now, observe that since $z\in L_{0}^{2}(\Omega)$, $\boldsymbol{\pi}_{h}z\in\boldsymbol{M}_{h}$. Consequently, by the definition of the elliptic projection (172),

[TABLE]

and therefore, by the continuity of $a_{\mathcal{D}}(\cdot,\cdot)$ (32), we have

[TABLE]

and the approximation properties of the projection $\boldsymbol{\pi}_{h}$ Lemma 3 yields

[TABLE]

The result now follows from a standard elliptic regularity argument and (179). \qed

For notational brevity, we denote by $\boldsymbol{v}^{n}=\boldsymbol{v}(t_{n})$ for a function $\boldsymbol{v}\in L^{1}(0,T;H^{2}(\mathcal{E}_{h})\times L^{2}(\Gamma_{h}))$ and $1\leq n\leq N$. We define the errors

[TABLE]

To derive the error equations required for the a priori analysis, observe that the following consistency property holds.

Lemma 15** (Consistency)**

Let $(c,\mu)$ be the exact solution to (1). Define $\boldsymbol{c}=(c,c|_{\Gamma_{h}})$ and $\boldsymbol{\mu}=(\mu,\mu|_{\Gamma_{h}})$. Under the same regularity assumptions as Theorem 4, we have

[TABLE]

By the consistency property (15) and the definition of the numerical scheme (81), for any $\boldsymbol{\chi}_{h},\boldsymbol{\phi}_{h}\in\boldsymbol{S}_{h}$, we obtain the following

[TABLE]

By the definition of the elliptic projection (172) and the orthogonal $L^{2}$-projection (42), we equivalently have for any $\boldsymbol{\chi}_{h},\boldsymbol{\phi}_{h}\in\boldsymbol{S}_{h}$

[TABLE]

In the proof of Theorem 4, a major difficulty is in handling the first term in the right-hand side of (193b), see term $T_{4}$ in (210). The following result is critical in controlling that term.

Lemma 16

Fix $1\leq n\leq N$ and set for brevity $Q=(c^{n})^{3}-(c_{h}^{n})^{3}$. Define $\boldsymbol{z}_{h}=(\pi_{h}Q,\{\pi_{h}Q\}|_{\Gamma_{h}})\in\boldsymbol{S}_{h}$. Assume that $c^{n}\in H^{s}(\Omega)$ for $2\leq s\leq k+1$. We have

[TABLE]

Proof 17

From the definition of the $\mathinner{\!\left\lVert\cdot\right\rVert}_{1,h}$-norm (12) and by the $H^{1}$-stability of the $L^{2}$-projection (44), we can write

[TABLE]

To bound the first term in the right-hand side, observe that

[TABLE]

so that on each element,

[TABLE]

*Thus, with Hölder’s inequality, we have *

[TABLE]

Since $c^{n}\in H^{2}(\Omega)$, we have by Sobolev’s embedding:

[TABLE]

By the discrete Gagliardo–Nirenberg inequality (75), (128) and (142), we have

[TABLE]

Therefore, the $H^{2}$ regularity of $c^{n}$ and (143) yield

[TABLE]

Using eq. 23, we write

[TABLE]

With the triangle inequality, we have

[TABLE]

With the approximation bound (173), we finally obtain

[TABLE]

We note that the terms on boundary faces in the right-hand side of (195) vanish. It remains to bound the terms corresponding to interior faces. For $e\in\mathcal{F}_{E}\cap\Gamma_{h}^{0}$, we have by the triangle inequality

[TABLE]

In the last inequality, we used (46). As $[(c^{n})^{3}]=0$, it suffices to estimate $[(c_{h}^{n})^{3}]$ to bound the third term on the right-hand side of (204). To this end, note that

[TABLE]

and therefore,

[TABLE]

where we have used (143). Thus, we have

[TABLE]

We have by (20) and the approximation properties (173)

[TABLE]

We then conclude by combining the bounds above. \qed

6.4 Proof of Theorem 4

Proof 18

From Proposition 1 and (172), we see that $\delta_{\tau}e_{c}^{n}\in M_{h}$. Testing (193a) with $\boldsymbol{\chi}_{h}^{n}=\boldsymbol{\mathcal{J}}_{h}(\delta_{\tau}e_{c}^{n})\in\boldsymbol{M}_{h}$, (193b) with $\boldsymbol{\phi}_{h}^{n}=\delta_{\tau}\boldsymbol{e}_{c}^{n}\in\boldsymbol{M}_{h}$ and subtracting, we find

[TABLE]

Then, it follows from Lemma 13 and the symmetry of the bilinear form $a_{\mathcal{D}}$ that

[TABLE]

Moreover, by the definition of the operator $\boldsymbol{\mathcal{J}}_{h}$, we have

[TABLE]

Therefore, with the coercivity of $a_{\mathcal{D}}$ (29), we obtain

[TABLE]

We proceed to bound each $T_{i}$, $i=1,\dots,5$. Fix $s$ such that $2\leq s\leq k+1$. By the Cauchy–Schwarz’s inequality, a Taylor expansion, the discrete Poincaré inequality (23) since $\boldsymbol{\mathcal{J}}_{h}(\delta_{\tau}e_{c,h}^{n})\in\boldsymbol{M}_{h}$, Young’s inequality, and (173), we have for any $\epsilon>0$

[TABLE]

We bound $T_{2}$ in a similar fashion:

[TABLE]

To bound $T_{3}$, we use the definition of the elliptic projection (172) and the boundedness of the bilinear form $a_{\mathcal{D}}$

[TABLE]

Therefore by the triangle inequality, Young’s inequality, (171), (173), we obtain

[TABLE]

To bound $T_{4}$, we first recall the definitions of $\boldsymbol{z}_{h}\in\boldsymbol{S}_{h}$ and $Q$ given in Lemma 16 and note that by the definition of the $L^{2}$-projection $\pi_{h}$, we have

[TABLE]

By Lemma 13, Lemma 16, the boundedness (30) of the bilinear form $a_{\mathcal{D}}$, and Young’s inequality, we have

[TABLE]

Finally, to bound $T_{5}$, observe that by Lemma 13, Lemma 16, triangle inequality, (20), (173), and Young’s inequality, we have

[TABLE]

Collecting the bounds on $T_{i}$, $i=1\dots 5$, choosing $\epsilon=C_{\mathrm{coer}}/5$, and rearranging, we find there are constants $C$ and $\tilde{C}$ independent of $h$, $\tau$, $c$, $\mu$, and $\kappa$ such that

[TABLE]

Multiplying by $2\tau$, summing from $n=1$ to $n=m$, noting that $\boldsymbol{e}_{c}^{0}=\mathbf{0}$ by construction, and applying a discrete Grönwall inequality yield for any $1\leq m\leq N$,

[TABLE]

under the assumption that $\tau$ is small enough, namely $\tau\leq(\kappa C_{\mathrm{coer}})/(4\tilde{C})$. Hence, the triangle inequality and Lemma 14 yield the bound on $\|\boldsymbol{c}^{n}-\boldsymbol{c}_{h}^{n}\|_{1,h}$. To obtain a bound on $\boldsymbol{e}_{\mu}^{n}$, test (193a) with $\boldsymbol{e}_{\mu}^{n}$. With the coercivity property of $a_{\mathcal{D}}$, we obtain

[TABLE]

The term $W_{1}$ is bounded by (126), (20), and Young’s inequality:

[TABLE]

Choosing $\boldsymbol{\chi}_{h}=(1,1)$ in (193a), we observe that

[TABLE]

since $\delta_{\tau}e_{c}^{n}\in M_{h}$. Thus, we bound $W_{2}$ with the Cauchy-Schwarz’s and Poincaré’s inequalities, the approximation properties of $\Pi_{h}$ and the usual Taylor expansions (similar to the bounds of $T_{1}$ and $T_{2}$) as follows:

[TABLE]

Similar to $T_{3}$, we derive that

[TABLE]

Combining the above bounds, multiplying by $\tau$, summing the resulting inequality from $n=1$ to $n=m$, and using (219) and the triangle inequality yields the error bound for the chemical potential $\mu$. \qed

7 Numerical experiment

To illustrate the rates of convergence experimentally, we consider a scenario taken from [17] involving the merging or separation of two droplets of one fluid surrounded by another. We take as the initial condition for the order parameter $c_{0}$:

[TABLE]

on the two dimensional domain $\Omega=(0,1)^{2}$. The region occupied by the droplets is indicated by $c=1$, while the region occupied by the surrounding fluid is indicated by $c=-1$. This example has been implemented using Netgen/NGSolve [38, 39]. We apply the HDG scheme (81) with $k=1$ for a sequence of mesh and time step sizes $h_{j}=1/2^{j}$ and $\tau_{j}=0.1/2^{2j}$ for $j=3,4,5,6$ and interface parameters $\kappa=1/2^{8},1/2^{10},1/2^{12}$ until the end time $T=0.1$. In the absence of an analytical solution, we instead compute the $L^{2}$ error $\mathinner{\lVert c_{h_{j}}-c_{h_{6}}\rVert}_{L^{2}(\Omega)}$ at the final time $T=0.1$ for $j=3,4,5$. The order of convergence is estimated via

[TABLE]

for $j=3,4,5$. Here, for notational brevity, we have suppressed the time index. The approximate solutions are plotted in Figure 1 and the approximate errors and rates of convergence are listed in Table 1. We observe merging of the two droplets for $\kappa=1/2^{8},1/2^{10}$ and separation of the two droplets for $\kappa=1/2^{12}$. In each case, the order of convergence in the $L^{2}$-norm appears to approach two. As we have taken $\tau_{j}=O(h_{j}^{2})$, we observe first order convergence in time and second order convergence in space as predicted by Theorem 4.

8 Conclusion

In this paper, we have analyzed a hybridized IPDG method for solving the mixed Cahn–Hilliard system combined with a convex-concave splitting of the chemical energy density and a first order implicit Euler method. We proved the unconditional unique solvability of the nonlinear algebraic system arising from our discretization using techniques from the theory of monotone operators. We showed the unconditional stability of the scheme for any $\mathcal{C}^{2}$ potential function, and established the $L^{\infty}$ stability of the order parameter for the Ginzburg–Landau potential on convex domains. Next, we derived optimal a priori error estimates in space and time for the Ginzburg–Landau potential in a mesh-dependent $H^{1}$-like norm. Finally, we observed the expected rate of convergence from our theoretical results through a numerical experiment.

Appendix 0.A Proof details for Lemma 5

Proof 19

We begin by proving (57). To simplify notation, define $\boldsymbol{\xi}=\boldsymbol{G}w_{h}-\boldsymbol{\pi}_{h}G(w_{h})$, $\boldsymbol{\eta}_{h}=\boldsymbol{\pi}_{h}G(w_{h})-\boldsymbol{G}_{h}\boldsymbol{w}_{h}$ so that $\boldsymbol{G}w_{h}-\boldsymbol{G}_{h}\boldsymbol{w}_{h}=\boldsymbol{\xi}+\boldsymbol{\eta}_{h}$. We have from (59) that

[TABLE]

Choosing $\boldsymbol{v}_{h}=\boldsymbol{\eta}_{h}\in\boldsymbol{M}_{h}$ and using the coercivity (29) and continuity (30) of the bilinear form $a_{\mathcal{D}}$ and the bound on $j_{0}$ (16), we have

[TABLE]

Using (171), Lemma 3 and elliptic regularity, we have

[TABLE]

and therefore,

[TABLE]

This shows (57). Next, we show (58) using a standard duality argument. Consider the following boundary value problem: find $z\in H^{2}(\Omega)\cap L_{0}^{2}(\Omega)$ such that

[TABLE]

Since $\Omega$ is convex, $z\in H^{2}(\Omega)$, we have the following identity for $\boldsymbol{z}=(z,z|_{\Gamma_{h}})$.

[TABLE]

As $z$ solves (228), $(\pi_{h}G(w_{h})-G_{h}\boldsymbol{w}_{h},G(w_{h}))_{\Omega}=(\pi_{h}G(w_{h})-G_{h}\boldsymbol{w}_{h},\pi_{h}Gw_{h})_{\Omega}$, and the bilinear form $a_{\mathcal{D}}$ is symmetric, we have that

[TABLE]

Let $\boldsymbol{\pi}_{h}z=(\pi_{h}z,\hat{\pi}_{h}z)\in\boldsymbol{M}_{h}$. By (59),

[TABLE]

Observe that the first term on the right-hand side of (231) is bounded by (32), (226), (227), (18), (171), (46) and elliptic regularity:

[TABLE]

To bound the second term on the right-hand side of (231), we use (16), (171) and elliptic regularity:

[TABLE]

*Therefore, *

[TABLE]

\qed

Appendix 0.B Proof of (126)

Proof 20

To prove (126), fix $w_{h}\in M_{h}$ and $v\in H^{1}(\mathcal{E}_{h})$. Define define $\boldsymbol{z}_{h}=(\pi_{h}v,\{\pi_{h}v\}|_{\Gamma_{h}})\in\boldsymbol{S}_{h}$. By the definition of the orthogonal $L^{2}$-projection $\pi_{h}\mathrel{\mathop{\mathchar 58\relax}}H^{1}(\mathcal{E}_{h})\to S_{h}$, the definition of the operator $\boldsymbol{\mathcal{J}}_{h}$, and the boundedness of the bilinear form $a_{\mathcal{D}}$ we have

[TABLE]

Observe that

[TABLE]

*where $\mathcal{E}_{e}$ is the set of the two neighboring elements of $e$, denoted by $E_{e,1}$ and $E_{e,2}$. Therefore, *

[TABLE]

Thus, by the stability of the $L^{2}$-projection in the DG norm we find

[TABLE]

The result follows. \qed