A Uniquely Solvable, Positivity-Preserving and Unconditionally Energy Stable Numerical Scheme for the Functionalized Cahn-Hilliard Equation with Logarithmic Potential

We consider the following functionalized Cahn-Hilliard (FCH) equation:

Here, we assume that $T>0$ is a given final time and $\Omega=(0,1)^{d}$ with the spatial dimension $1\leq d\leq 3$. For the sake of simplicity, below we analyze the system (1.1) in a periodic setting on $\Omega$ such that the variables $\phi$, $\mu$, $\omega$ satisfy periodic boundary conditions. With some minor modifications, our results can be extended to the case with homogeneous Neumann boundary conditions, or boundary conditions of mixed periodic-homogeneous Neumann type.

In the system (1.1), the scalar function $\phi:\Omega\times(0,T)\to\mathbb{R}$ denotes local concentrations of the two components in a binary mixture (e.g., an amphiphilic material in a solvent). Here we define $\phi$ as the difference of volume fractions. In view of its physical interpretation, only the values $\phi\in[-1,1]$ are admissible. The pure states correspond to the values $-1$ (the pure solvent) and $1$ (the pure amphiphile), whereas $\phi\in(-1,1)$ indicates the presence of a mixing state. The scalar function $\mu:\Omega\times(0,T)\to\mathbb{R}$ is referred to as the chemical potential, which is the first variational derivative of the following free energy functional:

where $\epsilon>0$, $\eta>0$ and $p\in\{1,2\}$ are given parameters. The energy (1.2) extends the classical Cahn-Hilliard energy in the literature (see [1])

where $F$ stands for the (Helmholtz) free energy density of mixing and its derivative is denoted by $f=F^{\prime}$. Then the scalar function $\omega:\Omega\times(0,T)\to\mathbb{R}$ in (1.1) corresponds to the first variational derivative of $E_{\mathrm{CH}}$, which is the chemical potential in the classical Cahn-Hilliard theory for binary mixtures.

In (1.2), the positive parameter $\epsilon\ll 1$ characterizes the width of inner structures, e.g., the (diffuse) interfaces between two components. Next, let us pay some more attention on the parameter $\eta\in\mathbb{R}$, which plays an important role in the modelling and analysis. Indeed, the sign of $\eta$ can change fundamentally the structure of the problem as well as the physically motivating examples [2]. When $\eta=0$, the energy $E$ reduces to the well-known phase-field formulation of the Willmore functional (PFW) that approximates the Canham-Helfrich bending energy of surfaces. The corresponding phase-field model has been used to study deformations of elastic vesicles subject to possible volume/surface constraints (see e.g., [3, 4]). When $\eta<0$, the energy $E$ is related to the so-called Cahn-Hilliard-Willmore (CHW) energy, which was introduced in [5, 6, 7] to investigate strong anisotropy effects arising in the growth and coarsening of thin films. In this paper, our interest will focus on the case $\eta>0$ such that $E$ is referred to as the functionalized Cahn-Hilliard (FCH) energy in the literature. It was first proposed in [8] to model phase separation of mixtures with an amphiphilic structure. Later on, the extended FCH model was used to describe nanoscale morphology changes in functionalized polymer chains [9], amphiphilic bilayer interfaces [2, 10, 11], and membranes undergoing pearling bifurcations [12, 13]. With a positive parameter $\eta$, the FCH energy (1.2) reflects the balance between the square of the variational derivative of $E_{\mathrm{CH}}$ against itself: the first term stabilizes bilayers, pores and micelle structures, while the second term promotes the growth of interfaces and competition between these geometries (see e.g., [14]). Minimizers of the FCH energy are approximate critical points of $E_{\mathrm{CH}}$, which however favor more surface area. This fact dramatically changes the nature of the energy landscape and presents rich morphological structures that are quite different from the phase separation process described by the classical Cahn-Hilliard energy. The FCH energy can naturally incorporate the propensity of the amphiphilic surfactant phase to drive the creation of stable bilayers, or homoclinic interfaces with an intrinsic width [2]. In particular, the bilayer structure has the distinct feature that it separate two identical phases by a thin region of a second phase. For recent works on the minimization problem of the FCH energy, bilayer-structures, pearled patterns, nanostructure morphology changes and network bifurcations, we refer to [13, 14, 15, 16, 17] and the references cited therein. Finally, we note that the exponent $p$ in (1.2) distinguishes some natural limits of the FCH energy, when $p=1$, it represents the strong functionalization, while the choice $p=2$ corresponds to the weak functionalization (see [10, 12, 15] for detailed discussions).

In this study, we work with the following thermodynamically relevant potential of logarithmic type:

where $\lambda\in\mathbb{R}$ is a given parameter. The potential function $F$ was proposed in the original work [1] on the phase separation of binary alloys and it also arises in the Flory–Huggins theory for polymer solutions [18]. The first order derivative of $F$ is denoted by $f$ such that

When $\lambda>2$, $F$ has a double-well structure with two minima $\pm\,r_{*}\in(-1,1)$, where $r_{*}$ is the positive root of the equation $f(r)=0$. Besides, we observe that $\lim_{r\to\pm 1}f(r)=\pm\infty$. This singular nature of $f$ near the pure states $\pm 1$ brings great challenges in the study of related phase-field equations, including the classical Cahn-Hilliard equation driven by $E_{\mathrm{CH}}$. On the other hand, to avoid possible difficulties caused by the singularity of $f$, a widely used approach is to work with an approximate polynomial like

which is sometimes called the “shallow quenching” approximation of the logarithmic potential (1.4) (see e.g., [19, 20]). Concerning the Cahn-Hilliard equation $\partial_{t}\phi=\nabla\cdot[\mathcal{M}\nabla(-\epsilon^{2}\Delta\phi+f(\phi))]$ with logarithmic potential (1.4), we refer to [19, 21, 22, 23, 24, 25, 26] for extensive theoretic analysis on well-posedness as well as long-time behavior of global solutions (see also the review articles [20, 27]), while for contributions in the numerical studies, we recall [28, 29, 30, 31, 32, 33, 34, 35, 36] and the references therein.

Under suitable boundary conditions (e.g., the periodic or Neumann type), the FCH equation (1.1) can be viewed as an $H^{-1}$ gradient flow of the FCH energy (1.2), which not only dissipates the energy but also preserves the mass along time evolution. As long as a solution $\phi$ exists on $(0,T)$ and is sufficiently regular, it holds

and

Concerning the theoretic analysis of equation (1.1), when $\eta>0$ and the periodic boundary conditions are imposed, Dai et al. [37] proved the existence of global weak solutions in the case with a general regular potential including (1.6) and a degenerate mobility $\mathcal{M}=\mathcal{M}(\phi)$ (see also [38]). For the case with the regular potential (1.6) and a constant mobility, existence and uniqueness of global solutions in the Gevrey class were established in [39] for $\eta\in\mathbb{R}$, again in the periodic setting. We also mention [40, 41] in which some theoretic results for the Cahn-Hilliard-Willmore equation (with $\eta<0$) can be found. The case associated with a logarithmic potential (1.4) is more difficult. In the recent work [42], Schimperna and Wu investigated the initial boundary value problem of (1.1) subject to homogeneous Neumann boundary conditions with $\eta\in\mathbb{R}$. After overcoming difficulties from the singular potential and its interaction with higher-order spatial derivatives, they proved existence, uniqueness and regularity of a global weak solution and characterized its long-time behavior. For a higher-order parabolic equation, in general we do not have maximum principle for its solution. Nevertheless, the singular character of the nonlinearity $f$ can help to keep the weak solutions of (1.1) staying away from the pure states at a point-wise level such that $-1<\phi(x,t)<1$ almost everywhere in $\Omega\times(0,T)$. This is sometimes referred to as a positivity-preserving property, since both $1+\phi>0$ and $1-\phi>0$ hold. When the spatial dimension is less than or equal to two, an improved conclusion has been drawn in [42]: if the initial datum $\phi_{0}$ is strictly separated from $\pm 1$, then there exists a uniform distance $\delta\in(0,1)$ such that $\|\phi(t)\|_{C(\overline{\Omega})}\leq 1-\delta$ for all $t\in[0,T]$. This fact is crucial for obtaining higher order spatial regularity of solutions, since the singular potential $f$ is now confined on a compact subset of $(-1,1)$ and thus is smooth. We note that similar properties have been proven for the classical Cahn-Hilliard equation with a logarithmic potential and for some extended systems, see for instance, [24, 25, 26] and the references therein.

Based on the analytic results obtained in [42], we are going to provide a detailed numerical analysis of the FCH equation (1.1) with logarithmic potential (1.4), subject to periodic boundary conditions. Observe that (1.1) is a sixth order parabolic equation for $\phi$ (see (2.3)–(2.4) below). Those nonlinear terms in the chemical potential $\mu$ such as $\beta(\phi)\beta^{\prime}(\phi)$, $\beta^{\prime\prime}(\phi)|\nabla\phi|^{2}$, $\Delta\beta(\phi)$ with $\beta(\phi)\triangleq f(\phi)+\lambda\phi$ (see (2.4)) present a highly nonlinear and singular nature of the problem, which turns out to be more involved than the fourth order Cahn-Hilliard equation. This feature yields great challenges in the numerical study. One obvious difficulty is to overcome the step restriction in the time-space discretization such that an explicit numerical scheme will require a rather restrictive Courant-Friedrichs-Lewy (CFL) condition like $\Delta t\leq Ch^{6}$, with $\Delta t$ and $h$ being the time and spatial step sizes, respectively. On the other hand, a fully implicit scheme like backward Euler may be difficult to prove its unique solvability and only conditionally energy stable. Our goal is to present a numerical scheme that is uniquely solvable, preserving the mass conservation, energy dissipation and the positivity property, under some mild refinement requirement. Besides, we shall apply some efficient numerical solvers to solve the numerical scheme in the long-time scale simulation.

There have been some progresses on the numerical study of the FCH equation (1.1) with $\eta>0$ and a polynomial potential like (1.6). In [14], Jones proposed a semi-implicit numerical scheme and proved its energy stability but the unique solvability was not given. He also made comparisons among semi-implicit, fully-implicit and explicit exponential time differencing (ETD) numerical schemes. A fully implicit scheme with pseudo-spectral approximation in space was introduced in Christlieb et al. [43], where the authors performed numerical tests to show its accuracy and efficiency but did not prove energy stability and solvability. In [44], Guo et al. presented a local discontinuous Galerkin method to overcome the difficulty associated with higher order spatial derivatives. Energy stability was shown for the semi-discrete (time-continuous) scheme therein. Later on, a convex-concave splitting scheme was proposed in Feng et al. [45], where the unique solvability and unconditional energy stability were theoretically verified. Besides, the authors performed the convergence analysis and applied an efficient preconditioned steepest descent algorithm to solve their scheme. Recently, theoretical analysis on the energy stability of a stabilized semi-implicit scheme for the FCH equation was carried out in [46]. We also refer to [47] for an unconditionally energy stable second-order numerical scheme based on the scalar auxiliary variable (SAV) approach, to [48] for a highly accurate, linear and unconditionally energy stable large time-stepping scheme, and to [49] for numerical comparison between those schemes based on the stable SAV approach and the classical BDF method.

It is worth mentioning that all the previous numerical studies mentioned above are associated with a regular (polynomial) potential. To the best of our knowledge, there is no result on the numerical analysis of the FCH equation (1.1) with the logarithmic potential (1.4) in the existing literature. Here we aim to make the first attempt in this direction. The main contributions of this paper are summarized below.

(1) We propose a first-order semi-implicit finite difference numerical scheme (see (2.18)–(2.19)), based on a suitable convex-concave decomposition of the highly nonlinear FCH energy (see Proposition 2.3).

One main difficulty comes from the troublesome term $-\int_{\Omega}\beta(\phi)\Delta\phi\,\mathrm{d}x$ in the FCH energy, which is neither convex nor concave (cf. (2.8)). Using integration by parts and the periodic boundary conditions, we can rewrite it as $\int_{\Omega}\beta^{\prime}(\phi)|\nabla\phi|^{2}\,\mathrm{d}x$ instead. Thanks to the strict separation property of $\phi$ (see Proposition 2.2), we are able to show that the new functional (and its discrete analogue) is strictly convex on a suitable closed set (cf. [42, 50] for conclusions in some weaker settings). This observation enables us to derive the desired convex splitting scheme that treats the convex part of $E$ implicitly and the concave part explicitly. Our argument is different from the case with the regular potential (1.6) studied in [45] such that no auxiliary term is required here to guarantee the necessary convexity.

(2) We prove that the proposed numerical scheme (2.18)–(2.19) is uniquely solvable, without any restriction on the size of the time step and the grid (see Theorem 3.1). Moreover, solutions to this fully discrete scheme preserve the mass at each time step (see Proposition 2.4) and satisfy the positivity-preserving property in the point-wise sense (see (3.1)).

The positivity-preserving property at the discrete level is crucial, since it guarantees the numerical scheme to be unconditionally well-defined. Concerning the Cahn-Hilliard equation with logarithmic potential (1.4), this issue was first tackled in [30] where Copetti and Elliott worked with a backward Euler scheme. Under certain constraint on the time step, they proved the unique solvability and positivity-preserving property of numerical solutions by applying a regularization of the singular potential and then passing to the limit. An alternative approach was recently introduced in Chen et al. [29], based on the convex-concave decomposition technique (see e.g., [51, 52]). They proposed a first order semi-implicit convex splitting scheme as well as a second order accurate scheme involving the implicit backward differential formulation (BDF). Then they overcame the time step restriction and theoretically justified the unique solvability and energy stability of those schemes. In particular, the logarithmic part was treated in an implicit way and its convexity as well as singularity helped to prevent the numerical solution from approaching the singular values $\pm 1$. Later in Chen et al. [28], a similar idea was applied to show the positivity-preserving property of a second order accurate scheme with a modified Crank-Nicolson approximation. Further applications can be found in recent works on the Cahn-Hilliard equation with a Flory-Huggins-de Gennes energy [53, 54], the Poisson-Nernst-Planck system [55] and the reaction-diffusion system with detailed balance [56] etc.

In this study, we successfully generalize the argument in [29] to the sixth order FCH equation (1.1). Roughly speaking, for any $n\in\mathbb{N}$, we look for a discrete solution $\phi^{n+1}$ to the numerical scheme (2.18)–(2.19) by solving a minimization problem form certain discrete energy $\mathcal{J}^{n}$ derived from the convex-concave energy decomposition (see (3.2)). If the numerical solution in the previous time step (denoted by $\phi^{n}$) takes all its grid values in $(-1,1)$ and its discrete average also belongs to $(-1,1)$, we can use the specific form of the logarithmic term $\beta(\phi)$ and in particular, its singular nature near $\pm 1$ to show that the solution $\phi^{n+1}$ that is a (global) minimum of $\mathcal{J}^{n}$, must exist and can only possibly occur at an interior point in the numerical solution variable domain $\mathcal{A}_{h}$ (see (3.4)). Besides, uniqueness of the numerical solution is a direct consequence of the strict convexity of the discrete energy function.

(3) We verify the unconditional energy stability of the proposed semi-implicit numerical scheme (2.18)–(2.19), see Theorem 4.2. The proof is based on the convex-concave decomposition of the discrete FCH energy given in Corollary 2.1. Moreover, the unconditional energy stability enables us to derive a uniform in time $H^{2}_{h}$-bound for the numerical solutions (see Corollary 4.1), which plays an important role in the subsequent convergence analysis.

(4) We perform an optimal rate convergence analysis for the numerical scheme (2.18)–(2.19) and establish an error estimate with $O(\Delta t+h^{2})$ accuracy in $l^{\infty}(0,T;L_{h}^{2})\cap l^{2}(0,T;H^{3}_{h})$, under a linear refinement requirement on the time step size such that $\Delta t\leq C_{1}h$ (see Theorem 5.3).

To obtain the expected error estimate under the above linear refinement constraint, a crucial step is to construct a higher order asymptotic expansion of the exact solution up to the second order temporal and spatial accuracy (see (5.11)). A similar idea was used in [55] for the optimal rate convergence analysis of the Poisson-Nernst-Planck system. Here in our case, based on the lower order error estimate in the discrete space $l^{\infty}(0,T;H_{h}^{-1})\cap l^{2}(0,T;H^{2}_{h})$ (see (5.35)) we can first derive a rough estimate on the modified numerical error function in $l^{\infty}(0,T;L_{h}^{2})$ (see (5.29)), which combined with the inverse inequality further leads to an error estimate in $l^{\infty}(0,T;L_{h}^{\infty})$ (see (5.30)), provided that $\Delta t$ and $h$ are suitably small and satisfy the linear constraint $\Delta t\leq C_{1}h$. Furthermore, this $L_{h}^{\infty}$ estimate allows us to show that all the numerical solutions can be kept away from the pure states $\pm 1$ with a uniform positive distance (see (5.31)). Hence, the singular term $\beta$ (and its derivatives) is confined on a compact subset of $(-1,1)$, becoming regular and uniformly bounded. With this crucial observation, we can proceed to derive the required higher order error estimate in $l^{\infty}(0,T;L_{h}^{2})\cap l^{2}(0,T;H^{3}_{h})$ by using the energy method.

The remaining part of this paper is organized as follows. In Section 2 we present some analytic results on the continuous problem and propose the numerical scheme. In Section 3, we provide detailed proofs on unique solvability of the scheme and the positivity-preserving property of numerical solutions. The unconditional energy stability is then established in Section 4. In Section 5, we perform an optimal rate convergence analysis and derive the error estimates. Numerical experiments are presented in Section 6 to verify the accuracy and efficiency of our scheme. The mass conservation, energy dissipation and positivity-preserving properties are verified numerically. In the final Section 7, we make some concluding remarks.