# A positivity-preserving and convergent numerical scheme for the binary
fluid-surfactant system

Yuzhe Qin  Research Center for Mathematics and Mathematics Education, Beijing
Normal University at Zhuhai, 519087, China  Laboratory of Mathematics and
Complex Systems (Ministry of Education), School of Mathematical Sciences,
Beijing Normal University, Beijing 100875, People’s Republic of China  Cheng
Wang Corresponding author. Mathematics Department, University of
Massachusetts, North Dartmouth, MA 02747, USA  Zhengru Zhang  Laboratory of
Mathematics and Complex Systems (Ministry of Education), School of
Mathematical Sciences, Beijing Normal University, Beijing 100875, People’s
Republic of China

###### Abstract

In this paper, we develop a first order (in time) numerical scheme for the
binary fluid surfactant phase field model. The free energy contains a double-
well potential, a nonlinear coupling entropy and a Flory-Huggins potential.
The resulting coupled system consists of two Cahn-Hilliard type equations.
This system is solved numerically by finite difference spatial approximation,
in combination with convex splitting temporal discretization. We prove the
proposed scheme is unique solvable, positivity-preserving and unconditionally
energy stable. In addition, an optimal rate convergence analysis is provided
for the proposed numerical scheme, which will be the first such result for the
binary fluid-surfactant system. Newton iteration is used to solve the discrete
system. Some numerical experiments are performed to validate the accuracy and
energy stability of the proposed scheme.

††Email addresses: yzqin@bnu.edu.cn(Y. Qin<EMAIL_ADDRESS>(C. Wang),
zrzhang@bnu.edu.cn(Z. Zhang)

Key words: Binary fluid-surfactant system, convex splitting, positivity-
preserving, unconditional energy stability, Newton iteration

##  1 Introduction

Two important characteristics of surfactants in binary fluid is that they can
move towards the fluid interface due to their amphiphilic structure and they
can reduce the interfacial tension and system energy [37]. Therefore,
surfactants have various applications in many fields such as biotechnology and
industry because of their features [2, 35]. In the past two decades, there
have been a number of excellent studies related to models with surfactants
[17, 27, 29, 53, 61, 63]. Often, there are two different ideas to model the
interfacial dynamics with surfactants. One is the sharp interface model which
has a long history dated back to one century ago [19, 48], and this kind of
model has been adopted in [24, 26]. In fact, sharp interface models have made
great progresses in explaining kinetics of diffusional phase transformations
and simulating multiphase systems with surfactants at one time. However, there
are some difficulties stemming from the interface interactions with various
complex processes during the course of phase transformations [35]. The other
approach method is the known phase-field method [16, 28, 18, 54]. This
approach makes use of an appropriate free energy functional to character the
interfacial dynamics, and it has been adopted to investigate the interfacial
dynamics with surfactants in [17, 29, 49, 51, 52]. In particular, phase-field
method was used to in [29] to study the phase transition behaviors of the
monolayer microemulsion system, formed by surfactant molecules. Generally, the
free energy of binary fluid-surfactant model consists of the following two
parts: the first part is the classical Ginzburg-Landau double well potential,
which is used to describe a binary mixture, and the other part, called
nonlinear coupling entropy term, has a historical evolution process, and is
used to account for the influence of the surfactant in boosting the formation
of interfaces. Laradji et al. in the pioneering work [29] introduced two phase
field variables to represent the local densities of the fluids, as well as the
local concentration of the surfactant, respectively. As mentioned in [27], an
extra diffusion term was added to prevent the model from becoming unbounded
and a Ginzburg-Landau type potential for the concentration variable to allow
the coexistence of the two bulk states. In order to restrict the range of the
concentration variable, the authors added the logarithmic Flory-Huggins
potential in [53], based on the nonlinear coupling entropy similar to [27,
29]. In consideration of penalizing the concentration to accumulate along the
fluid interface, the authors changed the nonlinear coupled entropy slightly in
[17]. In addition, a further modified model was considered in [50] by adding
the Flory-Huggins potential for the local concentration variable as well, in
comparison with the model in [17].

In this paper, we focus on constructing unconditionally energy stable
numerical schemes for the binary fluid-surfactant model in [50]. There have
been some works about numerical approximation to multi-phase models [7, 25].
Owing to the stiff nonlinear terms originated from the thin interface
thickness parameter, there are a lot of subtle difficulties to construct
numerical schemes with unconditional energy stability, especially for the
second order accurate (in time) scheme. Lots of efforts have been made to
solve these problems [7, 59, 66], etc. Since a simple fully implicit or
explicit type discretization brings extremely severe time step size constraint
on the interfacial width [1, 15, 47], a semi-implicit method was adopted in
[50]. However, the author mentioned that it suffers from a small CFL
conditional number. Recently, Gu et al. in [20] constructed an energy stable
finite difference scheme for the binary fluid-surfactant system, which is
based on the convex splitting approach [13, 43, 56, 58]: implicit treatment
for the convex part and explicit treatment for the concave part. Meanwhile, it
is observed that, the convexity analysis for one mixed term has not been
theoretically justified in [20], due to the multi variables involved in the
system. In addition, the positivity-preserving property has not been
theoretically proved, so that the well-defined nature of the numerical scheme
is not available. More recently, Yang et al. constructed the linear and stable
schemes for the binary fluid-surfactant system with constant mobility in [63],
using the invariant energy quadratization (IEQ) technique [6, 21, 23, 60, 62,
64, 65, 68]. In this approach, the free energy is transformed into an
equivalent quadratic form by introducing appropriate auxiliary variables, and
all nonlinear terms in this system are treated semi-explicitly [63]. The
energy stability has been derived for the IEQ method, while such a stability
has to be based on an alternate energy functional (involved with auxiliary
variables), not for the original energy functional, as always in the IEQ
approach. Moreover, the positivity-preserving property is not available to the
IEQ-based numerical method, because of the explicit treatment for the
nonlinear logarithmic term. In addition to the IEQ idea, Zhu et al. proposed
the scalar auxiliary variable (SAV) method to the surfactant model in [69],
following similar ideas in [46, 45]. The SAV approach introduces a constant-
coefficient linear equation to solve at each time step, and the energy
stability could be derived for an alternate energy involved with a scalar
variable. The convergence and error estimate for the SAV approach, for the
typical Cahn-Hilliard equation with double-well potential, has also been
established in recent works [33, 44]. However, an application of the SAV
approach to the surfactant model could not overcome the difficulty to
theoretically justify the positivity-preserving property, due to the explicit
treatment of the logarithmic term. In turn, the Flory-Huggins energy potential
has to be re-defined and extended around and beyond the singular phase
variable values. Also see a more recent work [39] of SAV-based numerical
algorithm for the surfactant model.

Among the existing numerical methods, different approaches have different
advantages. Here, we pay attention to the the convex splitting approach,
originated from the pioneering work of Eyre [13]. The idea is that the energy
admits a splitting into purely convex and concave parts, that is,
$E=E_{c}-E_{e}$, where $E_{c}$ and $E_{e}$ are both convex. Such an idea has
had wide applications in various gradient flow models, including either first
or second order accurate schemes. See the related works for the phase field
crystal (PFC) and the modified PFC (MPFC) equation [57, 55], the epitaxial
thin film growth models [4, 14], and the Cahn-Hilliard flow coupled with fluid
motion [3, 8], etc.

Meanwhile, there have been extensive works of linear numerical schemes for the
Cahn-Hilliard and epitaxial thin film equations [30, 31, 32], in which
stabilized implicit-explicit approach has been applied, and energy stability
has been theoretically provided. In the case of a Flory-Huggins energy
potential with singular logarithmic terms, the positivity-preserving property
has been recently established in [12] for the corresponding Allen-Cahn
equation, based on the maximum principle arguments. The advantage of such a
linear scheme (corresponding to an explicit treatment of the nonlinear
logarithmic term) is associated with the computational efficiency, so that a
nonlinear Newton iteration is not required. On the other hand, this approach
works very well for the positivity preserving analysis for the Allen-Cahn
gradient flow, due to the availability of maximum principle, while its direct
extension to the Cahn-Hilliard gradient flow would face a serious theoretical
difficulty. In this paper, we design a uniquely solvable, positivity-
preserving, unconditionally energy stable, and first order in time convergent
scheme for the binary fluid-surfactant system, based on the convex-splitting
idea, combined with the centered difference spatial approximation. For the
theoretical analysis of the positivity-preserving property, we make use of the
singular nature of the logarithmic function, and prove that such a singular
nature prevents the numerical solution approaches the singular limit values,
following similar ideas of in the analysis for the Cahn-Hilliard model [5, 9,
10, 11], as well as the one for the Poisson-Nernst-Planck system [34, 38],
droplet liquid film model [67], etc. In addition, an optimal rate convergence
analysis is provided, which is the first such work for the surfactant model.
The key difficulty in such an analysis is associated with the logarithmic
potential term and the coupled term. In this article, we can make full use of
the convexity of energy associated with the nonlinear terms to directly deal
with all logarithmic terms and coupled terms, because the convexity of energy
indicates the corresponding nonlinear error inner product is always non-
negative.

The rest of the paper is organized as follows. In Section 2, we give a brief
introduction to the binary fluid-surfactant phase field model and state its
energy law. In Section 3, the numerical scheme is proposed and analyzed, and
we prove the unique solvability, positivity-preserving property, as well as
the energy stability. An optimal rate convergence estimate is also provided.
Some numerical experiments are presented in Section 4. Finally, some
conclusions are made in Section 5 .

##  2 The mathematical model: binary fluid-surfactant system

In this paper, we consider the two-dimensional (2-D) binary fluid-surfactant
system. With the domain given by
$\Omega=\left(0,L_{x}\right)\times\left(0,L_{y}\right)$, the binary fluid-
surfactant system is formulated as

$\displaystyle\phi_{t}$ $\displaystyle=M_{1}\Delta\mu_{\phi},$ (2.1a)
$\displaystyle\rho_{t}$
$\displaystyle=M_{2}\nabla\cdot\left(M\left(\rho\right)\nabla\mu_{\rho}\right),$
(2.1b) $\displaystyle\mu_{\phi}$ $\displaystyle=\frac{\delta
G}{\delta\phi}=\frac{f^{\prime}(\phi)}{\varepsilon}-\varepsilon\Delta\phi+\alpha\nabla\cdot\Big{(}\frac{\rho\nabla\phi}{|\nabla\phi|}\Big{)},$
(2.1c) $\displaystyle\mu_{\rho}$ $\displaystyle=\frac{\delta
G}{\delta\rho}=-\alpha|\nabla\phi|+\beta H^{\prime}(\rho),$ (2.1d)

with the periodic boundary condition and
$M\left(\rho\right)=\rho\left(1-\rho\right)$. The PDE system (2.1) corresponds
to the following free energy functional

$G\left(\phi,\rho\right)=\int_{\Omega}\Big{(}\frac{f\left(\phi\right)}{\varepsilon}+\frac{\varepsilon}{2}|\nabla\phi|^{2}+\frac{\alpha}{2}\left(\rho-|\nabla\phi|\right)^{2}+\beta
H\left(\rho\right)\Big{)}\mathrm{d}\bm{x},$ (2.2)

where

$f\left(\phi\right)=\frac{1}{4}\phi^{2}\left(1-\phi\right)^{2},\quad
H\left(\rho\right)=\rho\ln\rho+\left(1-\rho\right)\ln\left(1-\rho\right),$

and $\alpha,\beta,\varepsilon$ are all small positive parameters. In this
paper we assume $M\left(\rho\right)=1$ and $M_{1}=M_{2}=\mathcal{M}$ for
simplicity. Furthermore, to avoid the singularity in calculating the coupled
energy $\left(\rho-|\nabla\phi|\right)^{2}$, we use
$\sqrt{\phi_{x}^{2}+\phi_{y}^{2}+\delta^{2}}$ to approximate $|\nabla\phi|$.
And also, we add diffuse terms $\frac{\eta^{2}}{2}|\Delta\phi|^{2}$ and
$\frac{\xi}{2}|\nabla\rho|^{2}$ in the energy density, so that the new free
energy functional becomes

$G_{new}\left(\phi,\rho\right)=\int_{\Omega}\Big{(}\frac{f\left(\phi\right)}{\varepsilon}+\frac{\varepsilon}{2}|\nabla\phi|^{2}+\frac{\eta^{2}}{2}|\Delta\phi|^{2}+\frac{\xi}{2}|\nabla\rho|^{2}+\frac{\alpha}{2}\left(\rho-|\nabla\phi|\right)^{2}+\beta
H\left(\rho\right)\Big{)}\mathrm{d}\bm{x}.$ (2.3)

For simplicity, we still use $G$ to express $G_{new}$. In turn, the
corresponding chemical potentials become

$\displaystyle\mu_{\phi}$ $\displaystyle=\frac{\delta
G}{\delta\phi}=\frac{f^{\prime}(\phi)}{\varepsilon}-\varepsilon\Delta\phi+\eta^{2}\Delta^{2}\phi+\alpha\nabla\cdot\Big{(}\frac{\rho\nabla\phi}{|\nabla\phi|}\Big{)},$
(2.4a) $\displaystyle\mu_{\rho}$ $\displaystyle=\frac{\delta
G}{\delta\rho}=-\alpha|\nabla\phi|-\xi\Delta\rho+\beta H^{\prime}(\rho).$
(2.4b)

As always in the gradient system, the energy dissipation property is always
valid:

$\frac{\mathrm{d}}{\mathrm{d}t}G\left(\phi\left(t\right),\rho\left(t\right)\right)\leq
0,\quad t>0.$

Besides, the appearance of the Flory-Huggins energy indicates a positivity
property for the density variable, $0<\rho<1$ at a point-wise level. Our
primary aim is to develop a numerical scheme inheriting these properties at a
theoretical level.

######  Remark 2.1.

Typically, in the context of physical models, the Dirichlet energy $\|\nabla
u\|^{2}$ represents surface tension, whereas the higher order term $\|\Delta
u\|^{2}$ represents bending rigidity. In general, it may be assumed that all
orders of the energy density are represented in the expansion of the energy

$E(u)=f(u)+a_{0}u^{2}+a_{2}\|\nabla u\|^{2}+a_{4}\|\Delta u\|^{2}+...,$

where f may be non-quadratic, and the coefficients $a_{k}$ may depend upon
$u$, or derivatives of $u$ but are usually constants. On the other hand, it is
typical to drop all higher order terms beyond those that are
physically/mathematically necessary to make the PDE well posed. One would
assume that the corresponding coefficients are sufficiently small so as to
justify this. In most cases, it may be more reasonable to assume that
$a_{2}=0$ (or is negligible) and only keep terms of order higher than two.

For the binary fluid-surfactant system (2.1), combined with the physical
energy (2.2), we are able to construct a numerical scheme with an energy
stability, while the optimal rate convergence analysis will face essential
theoretical difficulties, due to the highly nonlinear and singular 1-Laplacian
term involved for the variable $\phi$. To overcome this subtle difficulty, we
add an additional bi-harmonic diffusion term for $\phi$, as well as a regular
diffusion term for the $\rho$, in the energy representation (2.3). As a
result, both the energy stability and optimal rate convergence analysis could
be theoretically justified, as will be demonstrated in the later section. In
particular, an optimal rate convergence analysis will be the first such result
for the binary fluid-surfactant system. Such an approach of adding higher
order bi-harmonic diffusion process has been reported in many related
nonlinear physical systems, in which the diffusion terms have played essential
roles in the stability analysis, such as the planetary geostrophic equations
of oceanic geophysical fluid model [36, 40, 41, 42], etc.

##  3 The numerical scheme

In this section, we present a convex-concave decomposition of the energy
$\eqref{energy}$, and propose a convex splitting scheme based on such a
decomposition. The unique solvability, energy stability, positivity-preserving
property will be analyzed afterward.

###  3.1 The convex-concave decomposition of the energy

######  Lemma 3.1.

Suppose that $\Omega=\left(0,L_{x}\right)\times\left(0,L_{y}\right)$ and
$\phi,\rho:\Omega\to\mathbb{R}$ are periodic and sufficiently regular. Define
the following energy functionals

$\displaystyle G_{c}\left(\phi,\rho\right)=$
$\displaystyle\int_{\Omega}\frac{1}{4\varepsilon}(\phi-\frac{1}{2})^{4}+\frac{1}{64\varepsilon}+\frac{\varepsilon}{2}|\nabla\phi|^{2}+\frac{\eta^{2}}{2}|\Delta\phi|^{2}+\beta
H\left(\rho\right)+\frac{\xi}{2}|\nabla\rho|^{2}$
$\displaystyle+\frac{\alpha}{2}\left\\{\left(\rho-|\nabla^{\delta}\phi|\right)^{2}+(\sqrt{2}-1)\rho^{2}+\frac{1}{\delta}|\nabla\phi|^{2}\right\\}\mathrm{d}\bm{x},$
$\displaystyle G_{e}\left(\phi,\rho\right)=$
$\displaystyle\int_{\Omega}\frac{1}{8\varepsilon}(\phi-\frac{1}{2})^{2}+\frac{\alpha}{2}\left\\{(\sqrt{2}-1)\rho^{2}+\frac{1}{\delta}|\nabla\phi|^{2}\right\\}\mathrm{d}\bm{x},$

with $|\nabla^{\delta}\phi|:=(|\nabla\phi|^{2}+\delta^{2})^{\frac{1}{2}}$.
Then $G_{c}\left(\phi,\rho\right)$ and $G_{e}\left(\phi,\rho\right)$ are both
convex with respect to $\phi$ and $\rho$, with
$G\left(\phi,\rho\right)=G_{c}\left(\phi,\rho\right)-G_{e}\left(\phi,\rho\right)$.

###### Proof.

We focus on the convexity analysis of $G_{c}\left(\phi,\rho\right)$ and
$G_{e}\left(\phi,\rho\right)$. Let

$\displaystyle
e_{c}\left(\phi,\phi_{x},\phi_{y},\Delta\phi,\rho,\rho_{x},\rho_{y}\right)=$
$\displaystyle\frac{1}{4\varepsilon}(\phi-\frac{1}{2})^{4}+\frac{1}{64\varepsilon}+\frac{\varepsilon}{2}|\nabla\phi|^{2}+\frac{\eta^{2}}{2}|\Delta\phi|^{2}+\beta
H\left(\rho\right)+\frac{\xi}{2}|\nabla\rho|^{2}$
$\displaystyle+\frac{\alpha}{2}\left\\{\left(\rho-|\nabla\phi|\right)^{2}+(\sqrt{2}-1)\rho^{2}+\frac{1}{\delta}|\nabla\phi|^{2}\right\\},$
$\displaystyle
e_{e}\left(\phi,\phi_{x},\phi_{y},\Delta\phi,\rho,\rho_{x},\rho_{y}\right)=$
$\displaystyle\frac{1}{8\varepsilon}(\phi-\frac{1}{2})^{2}+\frac{\alpha}{2}\left\\{(\sqrt{2}-1)\rho^{2}+\frac{1}{\delta}|\nabla\phi|^{2}\right\\}.$

We also denote

$\displaystyle e_{c_{1}}\left(\bm{v}\right)$ $\displaystyle\triangleq
e_{c1}\left(\phi,\phi_{x},\phi_{y},\Delta\phi,\rho,\rho_{x},\rho_{y}\right)=\frac{1}{4\varepsilon}(\phi-\frac{1}{2})^{4}+\frac{1}{64\varepsilon}+\frac{\varepsilon}{2}|\nabla\phi|^{2}+\frac{\eta^{2}}{2}|\Delta\phi|^{2}+\beta
H\left(\rho\right)+\frac{\xi}{2}|\nabla\rho|^{2},$ $\displaystyle
e_{c_{2}}\left(\bm{v}\right)$ $\displaystyle\triangleq
e_{c2}\left(\phi,\phi_{x},\phi_{y},\Delta\phi,\rho,\rho_{x},\rho_{y}\right)=\frac{\alpha}{2}\left\\{\left(\rho-|\nabla^{\delta}\phi|\right)^{2}+(\sqrt{2}-1)\rho^{2}+\frac{1}{\delta}|\nabla\phi|^{2}\right\\},$
$\displaystyle e_{e}\left(\bm{v}\right)$ $\displaystyle\triangleq
e_{e}\left(\phi,\phi_{x},\phi_{y},\Delta\phi,\rho,\rho_{x},\rho_{y}\right)=\frac{1}{8\varepsilon}(\phi-\frac{1}{2})^{2}+\frac{\alpha}{2}\left\\{(\sqrt{2}-1)\rho^{2}+\frac{1}{\delta}|\nabla\phi|^{2}\right\\},$

where

${\bm{v}}=\left(v_{1},v_{2},v_{3},v_{4},v_{5},v_{6},v_{7}\right)\triangleq\left(\phi,\phi_{x},\phi_{y},\Delta\phi,\rho,\rho_{x},\rho_{y}\right).$

Then we have

$e_{c}\left(\bm{v}\right)=e_{c_{1}}\left(\bm{v}\right)+e_{c_{2}}\left(\bm{v}\right),\quad
G_{c}\left(\phi,\rho\right)=\int_{\Omega}e_{c}\left(\bm{v}\right)\mathrm{d}\bm{x},\quad
G_{e}\left(\phi,\rho\right)=\int_{\Omega}e_{e}\left(\bm{v}\right)\mathrm{d}\bm{x},$

and the following inequalities are derived:

$\displaystyle\partial_{v_{1}}^{2}e_{c_{1}}\left(v_{1},v_{2},v_{3},v_{4},v_{5},v_{6},v_{7}\right)$
$\displaystyle=\frac{3}{\varepsilon}(\phi-\frac{1}{2})^{2}\geq 0,$
$\displaystyle\partial_{v_{2}}^{2}e_{c_{1}}\left(v_{1},v_{2},v_{3},v_{4},v_{5},v_{6},v_{7}\right)$
$\displaystyle=\varepsilon>0,$
$\displaystyle\partial_{v_{3}}^{2}e_{c_{1}}\left(v_{1},v_{2},v_{3},v_{4},v_{5},v_{6},v_{7}\right)$
$\displaystyle=\varepsilon>0,$
$\displaystyle\partial_{v_{4}}^{2}e_{c_{1}}\left(v_{1},v_{2},v_{3},v_{4},v_{5},v_{6},v_{7}\right)$
$\displaystyle=\eta^{2}>0,$
$\displaystyle\partial_{v_{5}}^{2}e_{c_{1}}\left(v_{1},v_{2},v_{3},v_{4},v_{5},v_{6},v_{7}\right)$
$\displaystyle=\frac{\beta}{\rho\left(1-\rho\right)}>0,$
$\displaystyle\partial_{v_{6}}^{2}e_{c_{1}}\left(v_{1},v_{2},v_{3},v_{4},v_{5},v_{6},v_{7}\right)$
$\displaystyle=\xi>0,$
$\displaystyle\partial_{v_{7}}^{2}e_{c_{1}}\left(v_{1},v_{2},v_{3},v_{4},v_{5},v_{6},v_{7}\right)$
$\displaystyle=\xi>0,$
$\displaystyle\partial_{v_{1}}^{2}e_{e}\left(v_{1},v_{2},v_{3},v_{4},v_{5},v_{6},v_{7}\right)$
$\displaystyle=\frac{1}{4\varepsilon}>0,$
$\displaystyle\partial_{v_{2}}^{2}e_{e}\left(v_{1},v_{2},v_{3},v_{4},v_{5},v_{6},v_{7}\right)$
$\displaystyle=\frac{\alpha}{\delta}>0,$
$\displaystyle\partial_{v_{3}}^{2}e_{a}\left(v_{1},v_{2},v_{3},v_{4},v_{5},v_{6},v_{7}\right)$
$\displaystyle=\frac{\alpha}{\delta}>0,$
$\displaystyle\partial_{v_{4}}^{2}e_{c_{1}}\left(v_{1},v_{2},v_{3},v_{4},v_{5},v_{6},v_{7}\right)$
$\displaystyle=0,$
$\displaystyle\partial_{v_{5}}^{2}e_{e}\left(v_{1},v_{2},v_{3},v_{4},v_{5},v_{6},v_{7}\right)$
$\displaystyle=(\sqrt{2}-1)\alpha>0.$

These facts imply that both $e_{c_{1}}\left(\bm{v}\right)$ and
$e_{e}\left(\bm{v}\right)$ are convex. To get the convexity of
$e_{c_{2}}\left(\bm{v}\right)$, we have to analyze the Hessian matrix of
$e_{c_{2}}\left(\hat{\bm{v}}\right)$ as follows, where
$\hat{\bm{v}}=\left(v_{2},v_{3},v_{5}\right)$:

$H\left(\hat{\bm{v}}\right)=\alpha\left[\begin{matrix}\displaystyle\sqrt{2}&-\frac{v_{2}}{\sqrt{v_{2}^{2}+v_{3}^{2}+\delta^{2}}}&-\frac{v_{3}}{\sqrt{v_{2}^{2}+v_{3}^{2}+\delta^{2}}}\\\
-\frac{v_{2}}{\sqrt{v_{2}^{2}+v_{3}^{2}+\delta^{2}}}&-\frac{v_{5}\left(v_{3}^{2}+\delta^{2}\right)}{\left(v_{2}^{2}+v_{3}^{2}+\delta^{2}\right)^{\frac{3}{2}}}+\frac{1}{\delta}+1&\frac{v_{2}v_{3}v_{5}}{\left(v_{2}^{2}+v_{3}^{2}+\delta^{2}\right)^{\frac{3}{2}}}\\\
-\frac{v_{3}}{\sqrt{v_{2}^{2}+v_{3}^{2}+\delta^{2}}}&\frac{v_{2}v_{3}v_{5}}{\left(v_{2}^{2}+v_{3}^{2}+\delta^{2}\right)^{\frac{3}{2}}}&-\frac{v_{5}\left(v_{2}^{2}+\delta^{2}\right)}{\left(v_{2}^{2}+v_{3}^{2}+\delta^{2}\right)^{\frac{3}{2}}}+\frac{1}{\delta}+1\end{matrix}\right].$
(3.3)

Since $0<\rho<1$, i.e. $v_{5}\in\left(0,1\right)$, a careful application of
calculus reveals that $H\left(\hat{\bm{v}}\right)$ is diagonally dominated, so
that it is non-negative definite. This in turn indicates the convexity of
$e_{c_{2}}$. Therefore, we obtain the following inequality, according to the
definition of convex function:

$\displaystyle
e_{c}\left(\lambda\bm{w}+\left(1-\lambda\right)\bm{v}\right)\leq\lambda
e_{c}\left(\bm{w}\right)+\left(1-\lambda\right)e_{c}\left(\bm{v}\right),$
(3.4a) $\displaystyle
e_{e}\left(\lambda\bm{w}+\left(1-\lambda\right)\bm{v}\right)\leq\lambda
e_{e}\left(\bm{w}\right)+\left(1-\lambda\right)e_{e}\left(\bm{v}\right),$
(3.4b)

where $\lambda\in\left(0,1\right),\bm{w},\bm{v}\in\mathbb{R}^{7}.$ Integrating
both sides of $\eqref{convex_c}$ and $\eqref{convex_a}$ leads to

$G_{c}\left(\lambda\phi_{1}+\left(1-\lambda\right)\phi_{2},\lambda\rho_{1}+\left(1-\lambda\right)\rho_{2}\right)\leq\lambda
G_{c}\left(\phi_{1},\rho_{1}\right)+\left(1-\lambda\right)G_{c}\left(\phi_{2},\rho_{2}\right),$

and

$G_{e}\left(\lambda\phi_{1}+\left(1-\lambda\right)\phi_{2},\lambda\rho_{1}+\left(1-\lambda\right)\rho_{2}\right)\leq\lambda
G_{e}\left(\phi_{1},\rho_{1}\right)+\left(1-\lambda\right)G_{e}\left(\phi_{2},\rho_{2}\right),$

which indicates that both $G_{c}\left(\phi,\rho\right)$ and
$G_{e}\left(\phi,\rho\right)$ are convex with respect to $\phi$ and $\rho$. ∎

As a generalization of the theorem presented in [57], the following lemma is
the foundation of energy stability for binary fluid-surfactant, or more
generally, for two variable functional. The proof is similar to [57], so we
skip it for the sake of brevity.

######  Lemma 3.2.

Assume $\phi,\varphi,\rho,\psi:\Omega\to\mathbb{R}$ are periodic and smooth
enough. If $G=G_{c}-G_{e}$ gives a convex-concave decomposition, then we have

$G\left(\phi,\rho\right)-G\left(\varphi,\psi\right)\leq\left(\delta_{\phi}G_{c}\left(\phi,\rho\right)-\delta_{\phi}G_{e}\left(\varphi,\psi\right),\phi-\varphi\right)_{L^{2}}+\left(\delta_{\rho}G_{c}\left(\phi,\rho\right)-\delta_{\rho}G_{e}\left(\varphi,\psi\right),\rho-\psi\right)_{L^{2}},$
(3.5)

where $\delta$ denotes the variational derivative.

Given a time step $\Delta t>0$, we construct the discrete-time, continuous-
space scheme of the binary fluid-surfactant system $\eqref{bfs}$ as follows

$\displaystyle\frac{\phi^{n+1}-\phi^{n}}{\Delta t}$
$\displaystyle=\mathcal{M}\Delta\mu_{\phi}^{n+1},$ (3.6a)
$\displaystyle\frac{\rho^{n+1}-\rho^{n}}{\Delta t}$
$\displaystyle=\mathcal{M}\Delta\mu_{\rho}^{n+1},$ (3.6b)
$\displaystyle\mu_{\phi}^{n+1}$
$\displaystyle=\delta_{\phi}G_{c}\left(\phi^{n+1},\rho^{n+1}\right)-\delta_{\phi}G_{e}\left(\phi^{n},\rho^{n}\right),$
(3.6c) $\displaystyle\mu_{\rho}^{n+1}$
$\displaystyle=\delta_{\rho}G_{c}\left(\phi^{n+1},\rho^{n+1}\right)-\delta_{\rho}G_{e}\left(\phi^{n},\rho^{n}\right).$
(3.6d)

###  3.2 The spatial discretization and the fully discrete numerical scheme

The centered difference approximation is applied to discretize the space. Here
we first recall some basic notations of this methodology, and we use the
similar notations and results for some discrete functions and operators
introduced in [22, 57]. Let $\Omega=(0,L_{x})\times(0,L_{y})$, where for
simplicity, we assume $L_{x}=L_{y}:=L>0$. Let $N\in\mathbb{N}$ be given, and
define the grid mesh size $h:=L/N$. Such a uniform mesh size assumption is
only for simplicity of presentation. We define the following two uniform,
infinite node sets with grid spacing $h>0$:

$E:=\\{p_{i+\frac{1}{2}}~{}|~{}i\in\mathbb{Z}\\},\quad
C:=\\{p_{i}~{}|~{}i\in\mathbb{Z}\\},$ (3.7)

where $p_{i}=p(i):=(i-\frac{1}{2})\cdot h$. Consider the following 2-D
discrete $N^{2}$-peroidic function spaces:

$\displaystyle\mathcal{C}_{per}:=\\{\nu:C\times
C\to\mathbb{R}|\nu_{i,j}=\nu_{i+\alpha N,j+\beta N},~{}\forall
i,j,\alpha,\beta\in\mathbb{Z}\\},$
$\displaystyle\mathcal{E}_{per}^{x}:=\\{\nu:E\times
C\to\mathbb{R}|\nu_{i+\frac{1}{2},j}=\nu_{i+\frac{1}{2}+\alpha N,j+\beta
N},~{}\forall i,j,\alpha,\beta\in\mathbb{Z}\\}.$

The spaces $\mathcal{E}_{per}^{y}$ can be analogously defined. Here we use the
Greek symbols $\nu_{i,j}=\nu(p_{i},p_{j})$, et cetera. The functions of
$\mathcal{C}_{per}$ are called $cell-centered$ functions, while the functions
of $\mathcal{E}_{per}^{x}$ and $\mathcal{E}_{per}^{y}$ are called
$x-direction$ and $y-direction$ $edge-centered$ functions, respectively. We
also define the mean zero space

$\mathring{\mathcal{C}}_{per}:=\left\\{\nu\in\mathcal{C}_{per}~{}|~{}0=\bar{\nu}:=\frac{h^{2}}{|\Omega|}\sum_{i,j=1}^{m}\nu_{i,j}\right\\}.$

Additionally, we denote
$\vec{\mathcal{E}}_{per}:=\mathcal{E}_{per}^{x}\times\mathcal{E}_{per}^{y}$.

Next, the important difference and average operators are introduced on the
function spaces:

$\displaystyle
A_{x}\nu_{i+\frac{1}{2},j}:=\frac{1}{2}(\nu_{i+1,j}+\nu_{i,j}),\quad
D_{x}\nu_{i+\frac{1}{2},j}:=\frac{1}{h}(\nu_{i+1,j}-\nu_{i,j}),$
$\displaystyle
A_{y}\nu_{i,j+\frac{1}{2}}:=\frac{1}{2}(\nu_{i,j+1}+\nu_{i,j}),\quad
D_{y}\nu_{i,j+\frac{1}{2}}:=\frac{1}{h}(\nu_{i,j+1}-\nu_{i,j}),$

with $A_{x},D_{x}:\mathcal{C}_{per}\to\mathcal{E}_{per}^{x}$ and
$A_{y},D_{y}:\mathcal{C}_{per}\to\mathcal{E}_{per}^{y}$. Likewise,

$\displaystyle
a_{x}\nu_{i,j}:=\frac{1}{2}(\nu_{i+\frac{1}{2},j}+\nu_{i-\frac{1}{2},j}),\quad
d_{x}\nu_{i,j}:=\frac{1}{h}(\nu_{i+\frac{1}{2},j}-\nu_{i-\frac{1}{2},j}),$
$\displaystyle
a_{y}\nu_{i,j}:=\frac{1}{2}(\nu_{i,j+\frac{1}{2}}+\nu_{i,j-\frac{1}{2}}),\quad
d_{y}\nu_{i,j}:=\frac{1}{h}(\nu_{i,j+\frac{1}{2}}-\nu_{i,j-\frac{1}{2}}),$

with $a_{x},d_{x}:\mathcal{C}_{per}\to\mathcal{C}_{per}^{x}$ and
$a_{y},d_{y}:\mathcal{C}_{per}\to\mathcal{C}_{per}^{y}$. The discrete gradient
$\nabla_{h}:\mathcal{C}_{per}\to\vec{\mathcal{E}}_{per}$ is defined via

$\nabla_{h}\nu_{i,j}=(D_{x}\nu_{i+\frac{1}{2},j},D_{y}\nu_{i,j+\frac{1}{2}}),$

and the discrete divergence
$\nabla_{h}\cdot:\mathcal{E}_{per}\to\vec{\mathcal{C}}_{per}$ becomes

$\nabla_{h}\cdot\vec{f}_{i,j}=d_{x}f_{i,j}^{x}+d_{y}f_{i,j}^{y},\quad\mbox{for}\,\,\,\vec{f}=(f^{x},f^{y})\in\vec{\mathcal{E}}_{per}.$

The standard 2-D discrete Laplacian,
$\Delta_{h}:\mathcal{C}_{per}\to\mathcal{C}_{per}$, is given by

$\displaystyle\Delta_{h}\nu_{i,j}:$
$\displaystyle=\nabla_{h}\cdot(\nabla_{h}\nu)_{i,j}=d_{x}(D_{x}\nu)_{i,j}+d_{y}(D_{y}\nu)_{i,j}$
$\displaystyle=\frac{1}{h^{2}}(\nu_{i+1,j}+\nu_{i-1,j}+\nu_{i,j+1}+\nu_{i,j-1}-4\nu_{i,j}).$

More generally, if $\mathcal{D}$ is a periodic $scalar$ function that is
defined at all of the face center points and
$\vec{f}\in\vec{\mathcal{E}}_{per}$, assuming point-wise multiplication, we
may define

$\nabla_{h}\cdot(\mathcal{D}\vec{f})_{i,j}=d_{x}(\mathcal{D}f^{x})_{i,j}+d_{y}(\mathcal{D}f^{y})_{i,j}.$

Specifically, if $\nu\in\mathcal{C}_{per}$, then
$\nabla_{h}\cdot(\mathcal{D}\nabla_{h}):\mathcal{C}_{per}\to\mathcal{C}_{per}$
is defined point-wise via

$\nabla_{h}\cdot(\mathcal{D}\nabla_{h}\nu)_{i,j}=d_{x}(\mathcal{D}D_{x}\nu)_{i,j}+d_{y}(\mathcal{D}D_{y})_{i,j}.$

Now we are ready to introduce the following grid inner products:

$\displaystyle\langle\nu,\xi\rangle_{\Omega}:=h^{2}\sum_{i,j=1}^{N}\nu_{i,j}\xi_{i,j},\quad\nu,\xi\in\mathcal{C}_{per},$
$\displaystyle[\nu,\xi]_{x}:=\langle
a_{x}(\nu\xi),1\rangle_{\Omega},\quad\nu,\xi\in\mathcal{E}_{per}^{x},$
$\displaystyle[\nu,\xi]_{y}:=\langle
a_{y}(\nu\xi),1\rangle_{\Omega},\quad\nu,\xi\in\mathcal{E}_{per}^{y}.$
$\displaystyle[\vec{f}_{1},\vec{f}_{2}]_{\Omega}:=[f_{1}^{x},f_{2}^{x}]_{x}+[f_{1}^{y},f_{2}^{y}]_{y},\quad\vec{f}_{i}=(f_{i}^{x},f_{i}^{y})\in\vec{\mathcal{E}}_{per},\quad
i=1,2.$

We define the following norms for cell-centered functions. If
$\nu\in\mathcal{C}_{per}$, then
$\|\nu\|_{2}^{2}:=\langle\nu,\nu\rangle_{\Omega}$;
$\|\nu\|_{p}^{p}:=\langle|\nu|^{p},1\rangle_{\Omega}$, for $1\leq
p\leq\infty$, and $\|\nu\|_{\infty}:=\max\limits_{1\leq i,j\leq
N}|\nu_{i,j}|$. The norms of the gradient are defined as follows: for
$\nu\in\mathcal{C}_{per}$,

$\|\nabla_{h}\nu\|_{2}^{2}:=[\nabla_{h}\nu,\nabla_{h}\nu]_{\Omega}=[D_{x}\nu,D_{x}\nu]_{x}+[D_{y}\nu,D_{y}\nu]_{y},$

and, more generally, for $1\leq p\leq\infty$,

$\|\nabla_{h}\nu\|_{p}:=\left([|D_{x}\nu|^{p},1]_{x}+[|D_{y}\nu|^{p},1]_{y}\right)^{\frac{1}{p}},$

Higher order norms can be similarly introduced; for example,

$\|\nu\|_{H_{h}^{1}}^{2}:=\|\nu\|_{2}^{2}+\|\nabla_{h}\nu\|_{2}^{2},\quad\|\nu\|_{H_{h}^{2}}^{2}:=\|\nu\|_{H_{h}^{1}}^{2}+\|\Delta_{h}\nu\|_{2}^{2}.$

To facilitate the convergence analysis, we need to introduce a discrete
analogue of the space $H_{per}^{-1}(\Omega)$, as outlined in [56]. Suppose
that $\mathcal{D}$ is a positive, periodic scalar function defined at all of
the face center points. For any $\phi\in\mathcal{C}_{per}$, there exists a
unique $\psi\in\mathring{\mathcal{C}}_{per}$, that solves

$\mathcal{L}_{\mathcal{D}}(\psi):=-\nabla_{h}\cdot\left(\mathcal{D}\nabla_{h}\psi\right)=\phi-\bar{\phi},$
(3.8)

where, recall, $\bar{\phi}:=|\Omega|^{-1}\langle\phi,1\rangle_{\Omega}$. We
equip this space with a bilinear form: for any
$\phi_{1},\phi_{2}\in\mathring{\mathcal{C}}_{per}$, define

$\langle\phi_{1},\phi_{2}\rangle_{\mathcal{L}_{\mathcal{D}}^{-1}}:=[\mathcal{D}\nabla_{h}\psi_{1},\nabla_{h}\psi_{2}]_{\Omega},$
(3.9)

where $\psi_{i}\in\mathring{\mathcal{C}}_{per}$ is the unique solution to

$\mathcal{L}_{\mathcal{D}}(\psi_{i}):=-\nabla_{h}\cdot(\mathcal{D}\nabla_{h}\psi_{i})=\phi_{i},\quad
i=1,2.$ (3.10)

The following identity is easy to prove via summation-by-parts:

$\langle\phi_{1},\phi_{2}\rangle_{\mathcal{L}_{\mathcal{D}}^{-1}}=\langle\phi_{1},\mathcal{L}_{\mathcal{D}}^{-1}(\phi_{2})\rangle_{\Omega}=\langle\mathcal{L}_{\mathcal{D}}^{-1}(\phi_{1}),\phi_{2}\rangle_{\Omega},$
(3.11)

and since $\mathcal{L}_{\mathcal{D}}$ is symmetric positive definite,
$\langle\cdot,\cdot\rangle_{\mathcal{L}_{\mathcal{D}}^{-1}}$ is an inner
product on $\mathring{\mathcal{C}}_{per}$. When $\mathcal{D}\equiv 1$, we drop
the subscript and write $\mathcal{L}_{1}=\mathcal{L}$, and in this case we
usually write
$\langle\cdot,\cdot\rangle_{\mathcal{L}_{\mathcal{D}}^{-1}}=:\langle\cdot,\cdot\rangle_{-1,h}$.
In the general setting, the norm associated to this inner product is denoted
as
$\|\phi\|_{\mathcal{L}_{\mathcal{D}}^{-1}}:=\sqrt{\langle\phi,\phi\rangle_{\mathcal{L}_{\mathcal{D}}^{-1}}}$,
for all $\phi\in\mathring{\mathcal{C}}_{per}$, but, if $\mathcal{D}\equiv 1$,
we write $\|\cdot\|_{\mathcal{L}_{\mathcal{D}}^{-1}}=:\|\cdot\|_{-1,h}$.

With the preparations above, we turn to discuss the discrete energy and the
fully discrete scheme. Define the discrete energies
$E,E_{c},E_{e}:\mathcal{C}_{per}\times\mathcal{C}_{per}\to\mathbb{R}$ as

$\displaystyle E\left(\phi,\rho\right)=$ $\displaystyle
h^{2}\sum_{i,j=1}^{N}\Big{(}\frac{f\left(\phi_{i,j}\right)}{\varepsilon}+\frac{\varepsilon}{2}\left|\nabla_{h}\phi_{i,j}\right|^{2}+\frac{\eta^{2}}{2}\left|\Delta_{h}\phi_{i,j}\right|^{2}+\frac{\xi}{2}\left|\nabla_{h}\rho_{i,j}\right|^{2}$
$\displaystyle\qquad\quad+\frac{\varepsilon}{2}\left(\rho_{i,j}-{\cal
A}|\nabla_{h}^{\delta}\phi|_{i,j}\right)^{2}+\beta H(\rho_{i,j})\Big{)},$

in which ${\cal A}|\nabla_{h}^{\delta}\phi|$ is defined as

$({\cal
A}|\nabla_{h}^{\delta}\phi|_{i,j})^{2}=\frac{1}{2}((D_{x}\phi)_{i+\frac{1}{2},j}^{2}+(D_{x}\phi)_{i-\frac{1}{2},j}^{2}+(D_{y}\phi)_{i,j+\frac{1}{2}}^{2}+(D_{y}\phi)_{i,j-\frac{1}{2}}^{2})+\delta^{2}.$

The relevant discrete convex and concave energy functionals
$E_{c},E_{e}:\mathcal{C}_{per}\times\mathcal{C}_{per}\to\mathbb{R}$ are given
by

$\displaystyle E_{c}\left(\phi,\rho\right)=$ $\displaystyle
h^{2}\sum_{i,j=1}^{N}\left(\frac{1}{4\varepsilon}(\phi_{i,j}-\frac{1}{2})^{4}+\frac{1}{64\varepsilon}+\frac{\varepsilon}{2}|\nabla_{h}\phi_{i,j}|^{2}+\frac{\eta^{2}}{2}|\Delta_{h}\phi_{i,j}|^{2}+\beta
H\left(\rho_{i,j}\right)+\frac{\xi}{2}|\nabla_{h}\rho_{i,j}|^{2}\right.$
$\displaystyle\left.+\frac{\alpha}{2}\left(\left(\rho_{i,j}-{\cal
A}|\nabla_{h}^{\delta}\phi|_{i,j}\right)^{2}+(\sqrt{2}-1)\rho_{i,j}^{2}+\frac{1}{\delta}|\nabla_{h}\phi_{i,j}|^{2}\right)\right),$
(3.12) $\displaystyle E_{e}\left(\phi,\rho\right)=$ $\displaystyle
h^{2}\sum_{i,j=1}^{N}\left(\frac{1}{8\varepsilon}(\phi_{i,j}-\frac{1}{2})^{2}+\frac{\alpha}{2}\left((\sqrt{2}-1)\rho_{i,j}^{2}+\frac{1}{\delta}|\nabla_{h}\phi_{i,j}|^{2}\right)\right).$
(3.13)

In particular, a Hessian matrix could be similar formulated as the one given
by (3.3), for the following discrete function

$e_{c_{2},h}(\rho,v_{1},v_{2},v_{3},v_{4})=\frac{\alpha}{2}\Big{(}(\rho-\sqrt{\frac{1}{2}(v_{1}^{2}+v_{2}^{2}+v_{3}^{2}+v_{4}^{2})+\delta^{2}})^{2}+(\sqrt{2}-1)\rho^{2}+\frac{1}{2\delta}(v_{1}^{2}+v_{2}^{2}+v_{3}^{2}+v_{4}^{2})\Big{)}.$
(3.14)

A careful calculation reveals that, the corresponding $5\times 5$ Hessian
matrix is diagonally dominated, therefore non-negative definite. This in turn
leads to the convexity of the following discrete functional:

$E_{c_{2},h}(\rho,v_{1},v_{2},v_{3},v_{4})=\frac{\alpha}{2}h^{2}\sum_{i,j=1}^{N}\left((\rho_{i,j}-{\cal
A}|\nabla_{h}^{\delta}\phi|_{i,j})^{2}+(\sqrt{2}-1)\rho_{i,j}^{2}+\frac{1}{\delta}|\nabla_{h}\phi_{i,j}|^{2}\right).$
(3.15)

The convexity analysiss for the other parts of $E_{c}$ and $E_{e}$ is more
straightforward.

We follow the idea of convexity splitting and consider the following semi-
implicit, fully discrete scheme: given
$\phi^{n},\rho^{n}\in\mathcal{C}_{per}$, find
$\phi^{n+1},\rho^{n+1},\mu_{\phi}^{n+1},\mu_{\rho}^{n+1}\in\mathcal{C}_{per}$,
such that

$\displaystyle\frac{\phi^{n+1}-\phi^{n}}{\Delta t}$
$\displaystyle=\mathcal{M}\Delta_{h}\mu_{\phi}^{n+1},$ (3.16a)
$\displaystyle\frac{\rho^{n+1}-\rho^{n}}{\Delta t}$
$\displaystyle=\mathcal{M}\Delta_{h}\mu_{\rho}^{n+1},$ (3.16b)
$\displaystyle\mu_{\phi}^{n+1}$
$\displaystyle=\delta_{\phi}E_{c}\left(\phi^{n+1},\rho^{n+1}\right)-\delta_{\phi}E_{e}\left(\phi^{n},\rho^{n}\right),$
(3.16c) $\displaystyle\mu_{\rho}^{n+1}$
$\displaystyle=\delta_{\rho}E_{c}\left(\phi^{n+1},\rho^{n+1}\right)-\delta_{\rho}E_{e}\left(\phi^{n},\rho^{n}\right),$
(3.16d)

where

$\displaystyle\mu_{\phi}^{n+1}=$
$\displaystyle\frac{1}{\varepsilon}(\phi^{n+1}-\frac{1}{2})^{3}-(\varepsilon+\alpha+\frac{\alpha}{\delta})\Delta_{h}\phi^{n+1}+\eta^{2}\Delta_{h}^{2}\phi^{n+1}+\alpha\nabla_{h}\cdot\left({\cal
A}(\frac{\rho^{n+1}}{{\cal
A}|\nabla_{h}^{\delta}\phi^{n+1}|})\nabla_{h}\phi^{n+1}\right)$
$\displaystyle-\frac{1}{4\varepsilon}(\phi^{n}-\frac{1}{2})+\frac{\alpha}{\delta}\Delta_{h}\phi^{n},$
(3.17) $\displaystyle\mu_{\rho}^{n+1}=$
$\displaystyle-\xi\Delta_{h}\rho^{n+1}+\beta\left(\ln\rho^{n+1}-\ln\left(1-\rho^{n+1}\right)\right)+\sqrt{2}\alpha\rho^{n+1}-\alpha{\cal
A}\left|\nabla_{h}\phi^{n+1}\right|-\alpha(\sqrt{2}-1)\rho^{n}.$ (3.18)

Notice that $\nabla_{h}\cdot\left({\cal A}(\frac{\rho}{{\cal
A}|\nabla_{h}\phi|})\nabla_{h}\phi\right)$ is evaluated as follows

$\displaystyle\nabla_{h}\cdot\left({\cal A}(\frac{\rho}{{\cal
A}|\nabla_{h}^{\delta}\phi|})\nabla_{h}\phi\right)$ $\displaystyle=$
$\displaystyle\frac{1}{h}\Big{(}\frac{1}{2}(\frac{\rho_{i,j}}{{\cal
A}|\nabla_{h}^{\delta}\phi_{i,j}|}+\frac{\rho_{i+1,j}}{{\cal
A}|\nabla_{h}^{\delta}\phi_{i+1,j}|})(D_{x}\phi)_{i+\frac{1}{2},j}-\frac{1}{2}(\frac{\rho_{i,j}}{{\cal
A}|\nabla_{h}^{\delta}\phi_{i,j}|}+\frac{\rho_{i-1,j}}{{\cal
A}|\nabla_{h}^{\delta}\phi_{i-1,j}|})(D_{x}\phi)_{i-\frac{1}{2},j}\Big{)}$
$\displaystyle+\frac{1}{h}\Big{(}\frac{1}{2}(\frac{\rho_{i,j}}{{\cal
A}|\nabla_{h}^{\delta}\phi_{i,j}|}+\frac{\rho_{i,j+1}}{{\cal
A}|\nabla_{h}^{\delta}\phi_{i,j+1}|})(D_{y}\phi)_{i,j+\frac{1}{2}}-\frac{1}{2}(\frac{\rho_{i,j}}{{\cal
A}|\nabla_{h}^{\delta}\phi_{i,j}|}+\frac{\rho_{i,j-1}}{{\cal
A}|\nabla_{h}^{\delta}\phi_{i,j-1}|})(D_{y}\phi)_{i,j-\frac{1}{2}}\Big{)}.$
(3.19)

It is observed that the finite difference scheme is a system of nonlinear
equations with respect to $\phi^{n+1}$ and $\rho^{n+1}$, so we that have to
solve it iteratively. The theoretical properties of this scheme are analyzed
in the next few sections.

###  3.3 The positivity-preserving property

Of course, a point-wise bound for the grid function $\rho^{n+1}$, namely,
$0<\rho_{i,j}^{n+1}<1$, is needed to make sure the numerical scheme is well-
defined. The main theoretical result is stated below, which assures that there
exists a unique numerical solution for $\eqref{cs_1st_phi}$ and
$\eqref{cs_1st_rho}$, so that the given bound is satisfied.

######  Theorem 3.1.

Given $\phi^{n},\rho^{n}\in\mathcal{C}_{per}$, with $\|\rho^{n}\|_{\infty}\leq
M$, for some $M\geq 0$, and $0<\overline{\rho^{n}}<1$, there exists a unique
solution $\phi^{n+1},\rho^{n+1}\in\mathcal{C}_{per}$ to $\eqref{cs_1st}$, with
$0<\rho^{n+1}<1$ at a point-wise level.

Before the proof of the positivity-preserving property, we recall the
following lemma, cited from [5].

######  Lemma 3.3.

[5] Suppose that $\phi\in\mathring{\mathcal{C}}_{\mathrm{per}}$ and
$\|\phi\|_{\infty}\leq C_{1}$, then we have the following estimate:

$\displaystyle\|(-\Delta_{h})^{-1}\phi\|_{\infty}\leq C_{2}C_{1},$ (3.20)

where $C_{2}$ depends only on $\Omega$.

In addition, a few more preliminary estimates are needed for the positivity-
preserving analysis. The following discrete energy functional is introduced

$\displaystyle\mathcal{J}^{n}\left(\phi,\rho\right):=$
$\displaystyle\frac{1}{2\mathcal{M}\Delta
t}\|\phi-\phi^{n}\|_{-1,h}^{2}+\frac{1}{2\mathcal{M}\Delta
t}\|\rho-\rho^{n}\|_{-1,h}^{2}$ (3.21)
$\displaystyle+\frac{1}{4\varepsilon}\|\phi-\frac{1}{2}\|_{4}^{4}+\frac{\varepsilon}{2}\left\|\nabla_{h}\phi\right\|_{2}^{2}+\frac{\eta^{2}}{2}\left\|\Delta_{h}\phi\right\|_{2}^{2}+\frac{\xi}{2}\left\|\nabla_{h}\rho\right\|_{2}^{2}$
$\displaystyle+\frac{\alpha}{2}\left\\{\left\|\rho-{\cal
A}\left|\nabla_{h}^{\delta}\phi\right|\right\|_{2}^{2}+(\sqrt{2}-1)\|\rho\|_{2}^{2}+\frac{1}{\delta}\left\|\nabla_{h}\phi\right\|_{2}^{2}\right\\}$
$\displaystyle+\beta\left\langle\rho,\ln\rho\right\rangle_{\Omega}+\beta\left\langle
1-\rho,\ln\left(1-\rho\right)\right\rangle_{\Omega}+\langle\phi,f_{\phi}^{n}\rangle_{\Omega}+\langle\rho,f_{\rho}^{n}\rangle_{\Omega},$
$\displaystyle f_{\phi}^{n}=$
$\displaystyle-\frac{1}{4\varepsilon}(\phi^{n}-\frac{1}{2})-\frac{\alpha}{\delta}\Delta_{h}\phi^{n},\quad
f_{\rho}^{n}=-(\sqrt{2}-1)\alpha\rho^{n}.$

######  Lemma 3.4.

Set $M_{1}^{*}=\max|f_{\phi}^{n}|$,
$M_{2}^{*}=\max|f_{\rho}^{n}|\leq(\sqrt{2}-1)\alpha M$. We notice that
$M_{1}^{*}$ may be $h$, $\varepsilon$, $\delta$ and $\phi^{n}$ dependent. The
following inequality is available for the lower bound of
$\mathcal{J}^{n}(\phi,\rho)$, for $0<\rho<1$:

$\displaystyle\mathcal{J}^{n}(\phi,\rho)\geq\frac{1}{8\varepsilon}\|\phi-\frac{1}{2}\|_{4}^{4}-M_{3}^{n},\quad
M_{3}^{n}=(\frac{1}{2}(M_{1}^{*})^{2}+\frac{1}{2}M_{1}^{*}+\frac{1}{2}\varepsilon+\beta\ln
2+M_{2}^{*})|\Omega|,$ (3.22)

###### Proof.

In the expansion of $\mathcal{J}^{n}(\phi,\rho)$, it is observed that all the
terms are non-negative, except for the last four term. Then we get

$\mathcal{J}^{n}(\phi,\rho)\geq\frac{1}{4\varepsilon}\|\phi-\frac{1}{2}\|_{4}^{4}+\langle\phi,f_{\phi}^{n}\rangle_{\Omega}+\beta(\left\langle\rho,\ln\rho\right\rangle_{\Omega}+\left\langle
1-\rho,\ln\left(1-\rho\right)\right\rangle_{\Omega})+\langle\rho,f_{\rho}^{n}\rangle_{\Omega}.$
(3.23)

For the third and fourth terms on the right hand side of (3.23), the following
point-wise lower bound is available:

$\displaystyle\rho\ln\rho+(1-\rho)\ln(1-\rho)\geq
2\cdot\frac{1}{2}\ln\frac{1}{2}=-\ln 2,\quad\mbox{for $0<\rho<1$},$ (3.24)
$\displaystyle\mbox{which in turn leads
to}\quad\left\langle\rho,\ln\rho\right\rangle_{\Omega}+\left\langle
1-\rho,\ln\left(1-\rho\right)\right\rangle_{\Omega}\geq-|\Omega|\ln 2.$

For the last term on the right hand side of (3.23), the following lower bound
is valid, based on the fact that $0<\rho<1$:

$\langle\rho,f_{\rho}^{n}\rangle_{\Omega}\geq-\|f_{\rho}^{n}\|_{\infty}\cdot\|{\bf
1}\|_{1}\geq-M_{2}^{*}|\Omega|.$ (3.25)

For the first two terms on the right hand side of (3.23), we begin with the
following estimate:

$\displaystyle\begin{aligned}
\langle\phi,f_{\phi}^{n}\rangle_{\Omega}=&\langle\phi-\frac{1}{2},f_{\phi}^{n}\rangle_{\Omega}+\frac{1}{2}\langle{\bf
1},f_{\phi}^{n}\rangle_{\Omega}\geq-\frac{1}{2}(\|\phi-\frac{1}{2}\|_{2}^{2}+\|f_{\phi}^{n}\|_{2}^{2})-\frac{1}{2}M_{1}^{*}|\Omega|\\\
\geq&-\frac{1}{2}\|\phi-\frac{1}{2}\|_{2}^{2}-\frac{1}{2}(M_{1}^{*})^{2}|\Omega|-\frac{1}{2}M_{1}^{*}|\Omega|.\end{aligned}$
(3.26)

Meanwhile, the following lower bound is a direct application of Cauchy
inequality:

$\frac{1}{8\varepsilon}\|\phi-\frac{1}{2}\|_{4}^{4}-\frac{1}{2}\|\phi-\frac{1}{2}\|_{2}^{2}\geq-\frac{1}{2}\varepsilon|\Omega|,$
(3.27)

and its combination with (3.26) yields

$\displaystyle\frac{1}{8\varepsilon}\|\phi-\frac{1}{2}\|_{4}^{4}+\langle\phi,f_{\phi}^{n}\rangle_{\Omega}\geq-(\frac{1}{2}(M_{1}^{*})^{2}+\frac{1}{2}M_{1}^{*}+\frac{1}{2}\varepsilon)|\Omega|.$
(3.28)

Finally, a substitution of (3.24), (3.25) and (3.28) into (3.23) results in
(3.22). This completes the proof of Lemma 3.4. ∎

Now we proceed into the proof of Theorem 3.1.

###### Proof.

We denote $M_{0}^{n}=\mathcal{J}^{n}(\phi^{n},\rho^{n})$, a fixed constant
with given $(\phi^{n},\rho^{n})$. The numerical solution of $\eqref{cs_1st}$
is a minimizer of the discrete energy functional
$\mathcal{J}^{n}\left(\phi,\rho\right)$ (defined in (3.21)), over the
admissible set

$\displaystyle
A_{h}:=\left\\{\phi,\rho\in\mathcal{C}_{per}|-A^{*}<\phi<A^{*},\,0<\rho<1,\left\langle\phi-\bar{\phi}_{0}\right\rangle_{\Omega}=0,\left\langle\rho-\bar{\rho}_{0}\right\rangle_{\Omega}=0\right\\}\subset\mathbb{R}^{N^{2}}\times\mathbb{R}^{N^{2}},$
(3.29) $\displaystyle
A^{*}:=\max\Big{(}(8\varepsilon(M_{0}^{n}+M_{3}^{n}))^{\frac{1}{4}}h^{-\frac{1}{2}},\|\phi^{n}\|_{\infty}\Big{)}+1.$

We can observe that $\mathcal{J}^{n}$ is a strictly convex function over this
set.

To facilitate the analysis below, we transform the minimization problem into
an equivalent one. Consider the functional

$\displaystyle\mathcal{F}^{n}\left(\phi,\rho\right):=$
$\displaystyle\mathcal{J}^{n}\left(\phi+\bar{\phi}_{0},\rho+\bar{\rho}_{0}\right)$
$\displaystyle=$ $\displaystyle\frac{1}{2\mathcal{M}\Delta
t}\|\phi+\bar{\phi}_{0}-\phi^{n}\|_{-1,h}^{2}+\frac{1}{2\mathcal{M}\Delta
t}\|\rho+\bar{\rho}_{0}-\rho^{n}\|_{-1,h}^{2}$
$\displaystyle+\frac{1}{4\varepsilon}\|(\phi+\bar{\phi}_{0}-\frac{1}{2})^{2}\|_{2}^{2}+\frac{\varepsilon}{2}\left\|\nabla_{h}\phi\right\|_{2}^{2}+\frac{\eta^{2}}{2}\left\|\Delta_{h}\phi\right\|_{2}^{2}+\frac{\xi}{2}\left\|\nabla_{h}\rho\right\|_{2}^{2}$
$\displaystyle+\frac{\alpha}{2}\left\\{\left\|\rho+\bar{\rho}_{0}-{\cal
A}|\nabla_{h}^{\delta}\phi|\right\|_{2}^{2}+(\sqrt{2}-1)\left\|\rho+\bar{\rho}_{0}\right\|_{2}^{2}+\frac{1}{\delta}\left\|\nabla_{h}\phi\right\|_{2}^{2}\right\\}$
$\displaystyle+\beta\left\langle\rho+\bar{\rho}_{0},\ln\left(\rho+\bar{\rho}_{0}\right)\right\rangle_{\Omega}+\beta\left\langle
1-\rho-\bar{\rho}_{0},\ln\left(1-\rho-\bar{\rho}_{0}\right)\right\rangle_{\Omega}$
$\displaystyle+\langle\phi+\bar{\phi}_{0},f_{\phi}^{n}\rangle_{\Omega}+\langle\rho+\bar{\rho}_{0},f_{\rho}^{n}\rangle_{\Omega},$
(3.30)

defined on the set

$\mathring{A}_{h}:=\left\\{\left(\phi,\rho\right)\in\mathring{\mathcal{C}}_{per}\times\mathring{\mathcal{C}}_{per}|-A^{*}-\bar{\phi}_{0}<\phi<A^{*}-\bar{\phi}_{0},-\bar{\rho}_{0}<\rho<1-\bar{\rho}_{0}\right\\}\subset\mathbb{R}^{N^{2}}\times\mathbb{R}^{N^{2}}.$

If $(\phi,\rho)\in\mathring{A}_{h}$ minimizes $\mathcal{F}^{n}$, then
$(\varphi,\varrho):=(\phi+\bar{\phi}_{0},\rho+\bar{\rho}_{0})\in A_{h}$
minimizes $\mathcal{J}^{n}$ and $vice~{}versa$. Next, we prove that there
exists a minimizer of $\mathcal{F}^{n}$ over the domain $\mathring{A}_{h}$.
Consider the following closed domain: for
$\delta_{0}\in\left(0,\frac{1}{2}\right)$

$\displaystyle\mathring{A}_{h,\delta_{0}}:=$
$\displaystyle\left\\{\left(\phi,\rho\right)\in\mathring{\mathcal{C}}_{per}\times\mathring{\mathcal{C}}_{per}|-A^{*}-\bar{\phi}_{0}+\delta_{0}\leq\phi\leq
A^{*}-\bar{\phi}_{0}-\delta_{0},\,\delta_{0}-\bar{\rho}_{0}\leq\rho\leq
1-\delta_{0}-\bar{\rho}_{0}\right\\}$
$\displaystyle\subset\mathbb{R}^{N^{2}}\times\mathbb{R}^{N^{2}}.$

Since $\mathring{A}_{h,\delta_{0}}$ is a bounded, compact, and convex set in
the subspace $\mathring{\mathcal{C}}_{per}$, there exists a (not necessarily
unique) minimizer of $\mathcal{F}^{n}$ over $\mathring{A}_{h,\delta_{0}}$. The
key point of the positivity analysis is that, such a minimizer could not occur
on the boundary of $\mathring{A}_{h,\delta_{0}}$, if $\delta_{0}$ is
sufficiently small. To be more explicit, by the boundary of
$\mathring{A}_{h,\delta_{0}}$, we mean the locus of points
$\left(\phi,\psi\right)\in\mathring{A}_{h,\delta_{0}}$ such that
$\psi+\bar{\rho}_{0}=\delta_{0}$ or $1-\delta_{0}$, or
$\phi+\bar{\phi}_{0}=-A^{*}+\delta_{0}$ or $A^{*}-\delta_{0}$, precisely.

To get a contradiction, suppose that the minimizer of $\mathcal{F}^{n}$, call
it $\left(\phi^{\star},\rho^{\star}\right)$, occurs at a boundary point of
$\mathring{A}_{h,\delta_{0}}$. There is at least one grid point
$\vec{\alpha}_{0}=\left(i_{0},j_{0}\right)$ such that
$\rho_{\vec{\alpha}_{0}}^{\star}+\bar{\rho}_{0}=\delta_{0}$ or $1-\delta_{0}$,
or $\phi_{\vec{\alpha}_{0}}^{\star}+\bar{\phi}_{0}=-A^{*}+\delta_{0}$ or
$A^{*}-\delta_{0}$. Let us assume that
$\rho_{\vec{\alpha}_{0}}^{\star}+\bar{\rho}_{0}=\delta_{0}$, and denote by
$\vec{\alpha}_{1}=\left(i_{1},j_{1}\right)$ as the grid point at which
$\rho^{\star}$ achieves its maximum. By the fact that $\bar{\rho^{\star}}=0$,
it is obvious that $\rho_{\vec{\alpha}_{1}}^{\star}\geq 0$ and

$1-\delta_{0}\geq\rho_{\vec{\alpha}_{1}}^{\star}+\bar{\rho}_{0}\geq\bar{\rho}_{0}.$

Since $\mathcal{F}^{n}$ is smooth over $\mathring{A}_{h,\delta}$, for all
$\left(\varphi,\psi\right)\in\mathring{\mathcal{C}}_{per}$, the directional
derivative is

$\displaystyle\mathrm{d}_{s}\mathcal{F}^{n}\left(\phi^{\star}+s\varphi,\rho^{\star}+s\psi\right)|_{s=0}$
$\displaystyle=$ $\displaystyle\frac{1}{\mathcal{M}\Delta
t}\left\langle\left(-\Delta_{h}\right)^{-1}\left(\phi^{\star}+\bar{\phi}_{0}-\phi^{n}\right),\varphi\right\rangle_{\Omega}+\frac{1}{\mathcal{M}\Delta
t}\left\langle\left(-\Delta_{h}\right)^{-1}\left(\rho^{\star}+\bar{\rho}_{0}-\rho^{n}\right),\psi\right\rangle_{\Omega}$
$\displaystyle+\frac{1}{2\varepsilon}\langle\phi^{\star}+\bar{\phi}_{0}-\frac{1}{2},\varphi\rangle_{\Omega}-\varepsilon\left\langle\Delta_{h}\phi^{\star},\varphi\right\rangle_{\Omega}+\eta^{2}\left\langle\Delta_{h}^{2}\phi^{\star},\varphi\right\rangle_{\Omega}-\xi\left\langle\Delta_{h}\rho^{\star},\psi\right\rangle_{\Omega}$
$\displaystyle+\beta\left\langle\ln\left(\rho^{\star}+\bar{\rho}_{0}\right)-\ln\left(1-\rho^{\star}-\bar{\rho}_{0}\right),\psi\right\rangle_{\Omega}-\alpha\left\langle\nabla_{h}\cdot\Big{(}{\cal
A}(\frac{\rho^{\star}+\bar{\rho}_{0}}{{\cal
A}|\nabla_{h}^{\delta}\phi^{\star}|})\nabla_{h}\phi^{\star}\Big{)},\varphi\right\rangle_{\Omega}$
$\displaystyle-\alpha\left\langle{\cal
A}|\nabla_{h}^{\delta}\phi^{\star}|,\psi\right\rangle_{\Omega}+\alpha(\sqrt{2}-1)\left\langle\rho^{\star}+\bar{\rho}_{0},\psi\right\rangle_{\Omega}-\frac{\alpha}{\delta}\left\langle\Delta_{h}\phi^{\star},\varphi\right\rangle_{\Omega}+\langle
f_{\phi}^{n},\varphi\rangle_{\Omega}+\langle
f_{\rho}^{n},\psi\rangle_{\Omega}.$ (3.31)

Here, we take the direction $\varphi,\psi\in\mathring{C}_{per}$, such that

$\varphi=0,\quad\psi=\delta_{i,i_{0}}\delta_{j,j_{0}}-\delta_{i,i_{1}}\delta_{j,j_{1}}.$

Then the derivative may be expressed as

$\displaystyle\frac{1}{h^{2}}\mathrm{d}_{s}\mathcal{F}^{n}\left(\phi^{\star},\rho^{\star}+s\psi\right)|_{s=0}$
$\displaystyle=$ $\displaystyle\frac{1}{\mathcal{M}\Delta
t}\left(-\Delta_{h}\right)^{-1}\left(\rho^{\star}+\bar{\rho}_{0}-\rho^{n}\right)_{\vec{\alpha}_{0}}-\frac{1}{\mathcal{M}\Delta
t}\left(-\Delta_{h}\right)^{-1}\left(\rho^{\star}+\bar{\rho}_{0}-\rho^{n}\right)_{\vec{\alpha}_{1}}-\xi\left(\Delta_{h}\rho_{\vec{\alpha}_{0}}^{\star}-\Delta_{h}\rho_{\vec{\alpha}_{1}}^{\star}\right)$
$\displaystyle+\beta\ln\left(\rho_{\vec{\alpha}_{0}}^{\star}+\bar{\rho}_{0}\right)-\beta\ln\left(1-\rho_{\vec{\alpha}_{0}}^{\star}-\bar{\rho}_{0}\right)-\beta\ln\left(\rho_{\vec{\alpha}_{1}}^{\star}+\bar{\rho}_{0}\right)+\beta\ln\left(1-\rho_{\vec{\alpha}_{1}}^{\star}-\bar{\rho}_{0}\right)$
$\displaystyle-\alpha\left({\cal
A}|\nabla_{h}^{\delta}\phi_{\vec{\alpha}_{0}}^{\star}|-{\cal
A}|\nabla_{h}^{\delta}\phi_{\vec{\alpha}_{1}}^{\star}|\right)+\alpha(\sqrt{2}-1)\left(\rho_{\vec{\alpha}_{0}}^{\star}-\rho_{\vec{\alpha}_{1}}^{\star}\right)+(f_{\rho}^{n})_{\vec{\alpha}_{0}}-(f_{\rho}^{n})_{\vec{\alpha}_{1}}.$
(3.32)

For simplicity, now let us write
$\varrho^{\star}:=\rho^{\star}+\bar{\rho}_{0}$. Since
$\varrho_{\vec{\alpha}_{0}}^{\star}=\delta_{0}$ and
$\varrho_{\vec{\alpha}_{1}}^{\star}\geq\bar{\rho}_{0}$, we have

$\ln\left(\varrho_{\vec{\alpha}_{0}}^{\star}\right)-\ln\left(1-\varrho_{\vec{\alpha}_{0}}^{\star}\right)-\ln\left(\varrho_{\vec{\alpha}_{1}}^{\star}\right)+\ln\left(1-\varrho_{\vec{\alpha}_{1}}^{\star}\right)\leq\ln\frac{\delta_{0}}{1-\delta_{0}}-\ln\frac{\bar{\rho}_{0}}{1-\bar{\rho}_{0}}.$
(3.33)

Since $\varrho^{\star}$ takes a minimum at the grid point $\vec{\alpha}_{0}$,
with $\varrho_{\vec{\alpha}_{0}}^{\star}=\delta\leq\varrho_{i,j}^{\star}$, for
any $\left(i,j\right)$, and a maximum at the grid point $\vec{\alpha}_{1}$,
with $\varrho_{\vec{\alpha}_{1}}^{\star}\geq\varrho_{i,j}^{\star}$, for any
$\left(i,j\right)$,

$\Delta_{h}\rho_{\vec{\alpha}_{0}}^{\star}\geq
0,\quad\Delta_{h}\rho_{\vec{\alpha}_{1}}^{\star}\leq
0,\quad\rho_{\vec{\alpha}_{0}}^{\star}-\rho_{\vec{\alpha}_{1}}^{\star}\leq 0.$
(3.34)

For the numerical solution $\rho^{n}$ at the previous time step, the a priori
assumption $\left\|\rho^{n}\right\|_{\infty}\leq M$ indicates that

$-2M\leq\rho_{\vec{\alpha}_{0}}^{n}-\rho_{\vec{\alpha}_{1}}^{n}\leq 2M.$
(3.35)

According to Lemma 3.3, we obtain

$-4MC_{2}\leq\left(-\Delta_{h}\right)^{-1}\left(\varrho^{\star}-\rho^{n}\right)_{\vec{\alpha}_{0}}-\left(-\Delta_{h}\right)^{-1}\left(\varrho^{\star}-\rho^{n}\right)_{\vec{\alpha}_{1}}\leq
4MC_{2}.$ (3.36)

Denote $C_{3}=\max\left\\{{\cal
A}|\nabla_{h}^{\delta}\phi_{\vec{\alpha}_{0}}^{\star}|,{\cal
A}|\nabla_{h}^{\delta}\phi_{\vec{\alpha}_{1}}^{\star}|\right\\}$. Based on the
fact that $-A^{*}-\bar{\phi}_{0}<\phi^{*}<A^{*}-\bar{\phi}_{0}$ at a point-
wise level, we conclude that

${\cal A}|\nabla_{h}^{\delta}\phi_{\vec{\alpha}_{0}}^{\star}|,\,\,{\cal
A}|\nabla_{h}^{\delta}\phi_{\vec{\alpha}_{1}}^{\star}|\leq\frac{2A^{*}}{h}+1,\quad\mbox{so
that}\,\,\,C_{3}\leq\frac{2A^{*}}{h}+1.$ (3.37)

Then we have

$-\alpha C_{3}\leq-\alpha{\cal
A}\left(|\nabla_{h}^{\delta}\phi_{\vec{\alpha}_{0}}^{\star}|-|\nabla_{h}^{\delta}\phi_{\vec{\alpha}_{1}}^{\star}|\right)\leq\alpha
C_{3}.$ (3.38)

Consequently, a substitution of
$\eqref{estimate_log},\eqref{estimate_star},\eqref{estimate_n},\eqref{estimate_lap},\eqref{estimate_phi-1}$
into $\eqref{estimate_direction}$ yields the following bound on the
directional derivative:

$\frac{1}{h^{2}}\mathrm{d}_{s}\mathcal{F}^{n}\left(\phi^{\star},\rho^{\star}+s\psi\right)|_{s=0}\leq\beta\ln\frac{\delta_{0}}{1-\delta_{0}}-\beta\ln\frac{\bar{\rho}_{0}}{1-\bar{\rho}_{0}}+4MC_{2}(\mathcal{M}\Delta
t)^{-1}+2\alpha C_{3}+2M_{2}.$

We denote $D_{0}=4MC_{2}(\mathcal{M}\Delta t)^{-1}+2\alpha C_{3}+2M_{2}$.
Notice that $D_{0}$ is a constant for fixed $\Delta t$ and $h$, though it
becomes singular as $\Delta t\to 0$ and $h\to 0$. On the other hand, for any
fixed $\Delta t$ and $h$, we may choose $\delta_{0}\in\left(0,1/2\right)$
sufficiently small so that

$\beta\ln\frac{\delta_{0}}{1-\delta_{0}}-\beta\ln\frac{\bar{\rho}_{0}}{1-\bar{\rho}_{0}}+D_{0}<0.$
(3.39)

This in turn leads to the following inequality, provided $\delta_{0}$
satisfies $\eqref{estimate_C}$,

$\frac{1}{h^{2}}\mathrm{d}_{s}\mathcal{F}^{n}\left(\phi^{\star},\rho^{\star}+s\psi\right)|_{s=0}<0.$
(3.40)

As before, this contradicts the assumption that $\mathcal{F}^{n}$ has a
minimum at $\left(\phi^{\star},\rho^{\star}\right)$, since the directional
derivative is negative in a direction pointing into the interior of
$\mathring{A}_{h,\delta_{0}}$.

Using very similar arguments, we can also prove that the global minimum of
$\mathcal{F}^{n}$ over $\mathring{A}_{h,\delta_{0}}$ could not occur at a
boundary point $\left(\phi^{\star},\rho^{\star}\right)$ such that
$\rho_{\vec{\alpha}_{0}}^{\star}+\bar{\rho}_{0}=1-\delta_{0}$, for some
$\vec{\alpha}_{0}$, so that the grid function $\rho^{\star}$ has a global
maximum at $\vec{\alpha}_{0}$. The details are left to interested readers.

Moreover, if the global minimum of $\mathcal{F}^{n}$ over
$\mathring{A}_{h,\delta_{0}}$ could occurs at a boundary point
$\left(\phi^{\star},\rho^{\star}\right)$ such that
$\varphi_{\vec{\alpha}_{0}}^{*}=\phi_{\vec{\alpha}_{0}}^{\star}+\bar{\phi}_{0}=A^{*}-\delta_{0}$.
In turn, we apply Lemma 3.4 and obtain

$\displaystyle\begin{aligned}
\mathcal{F}^{n}(\phi^{*},\rho^{*})=&\mathcal{J}^{n}(\varphi^{*},\varrho^{*})\geq\frac{1}{8\varepsilon}\|\varphi^{*}-\frac{1}{2}\|_{4}^{4}-M_{3}^{n}\geq\frac{1}{8\varepsilon}h^{2}\cdot(\varphi_{\vec{\alpha}_{0}}^{*}-\frac{1}{2})^{4}-M_{3}^{n}\\\
\geq&\frac{1}{8\varepsilon}h^{2}(A^{*}-\delta_{0}-\frac{1}{2})^{4}-M_{3}^{n}>\frac{1}{8\varepsilon}h^{2}(A^{*}-1)^{4}-M_{3}^{n}\\\
\geq&\frac{1}{8\varepsilon}h^{2}\cdot(8\varepsilon(M_{0}^{n}+M_{3}^{n}))h^{-2}-M_{3}^{n}=M_{0}^{n}=\mathcal{J}^{n}(\phi^{n},\rho^{n}),\end{aligned}$
(3.41)

in which the definition of $A^{*}$ (in (3.29)) has been recalled. This
contradicts the assumption that $\mathcal{F}^{n}$ has a minimum at
$\left(\phi^{\star},\rho^{\star}\right)$.

Using similar arguments, a minimization point cannot occur at a boundary point
$\left(\phi^{\star},\rho^{\star}\right)$ such that
$\varphi_{\vec{\alpha}_{0}}^{*}=\phi_{\vec{\alpha}_{0}}^{\star}+\bar{\phi}_{0}=-A^{*}+\delta_{0}$.
The details are left to interested readers.

A combination of above four facts have indicated that, the global minimum of
$\mathcal{F}^{n}$ over $\mathring{A}_{h,\delta_{0}}$ could only possibly occur
at interior point
$\left(\phi,\rho\right)\in\left(\mathring{A}_{h,\delta_{0}}\right)^{\circ}\subset\left(\mathring{A}_{h}\right)^{\circ}$.
We conclude that there must be a solution $\left(\phi,\rho\right)\in A_{h}$
that minimizes $\mathcal{J}^{n}$ over $A_{h}$, which is equivalent to the
numerical solution of $\eqref{cs_1st}$. The existence of the numerical
solution is established.

In addition, since $\mathcal{J}^{n}$ is strictly convex function over $A_{h}$,
the uniqueness analysis for this numerical solution is straightforward. The
proof of Theorem 3.1 is completed. ∎

###  3.4 Unconditional energy stability

######  Theorem 3.2 (Energy stability).

For $n\geq 1$, the numerical scheme $\eqref{cs_1st}$ is unconditionally energy
stable, i.e.

$E\left(\phi^{n+1},\rho^{n+1}\right)\leq E\left(\phi^{n},\rho^{n}\right).$

###### Proof.

Let $\mathcal{L}=-\Delta_{h}$. Due to the mass conservation,
$\mathcal{L}^{-1}\left(\phi^{n+1}-\phi^{n}\right)$ and
$\mathcal{L}^{-1}\left(\rho^{n+1}-\rho^{n}\right)$ are well-defined. Taking a
discrete inner product with $\eqref{cs_1st_phi}$, $\eqref{cs_1st_rho}$,
$\eqref{cs_1st_mu1}$, $\eqref{cs_1st_mu2}$ by
$\mathcal{L}^{-1}\left(\phi^{n+1}-\phi^{n}\right)$,
$\mathcal{L}^{-1}\left(\rho^{n+1}-\rho^{n}\right)$, $\phi^{n+1}-\phi^{n}$ and
$\rho^{n+1}-\rho^{n}$, respectively, yields the following estimate

$\displaystyle 0=$ $\displaystyle\frac{1}{\mathcal{M}\Delta
t}\left\langle\phi^{n+1}-\phi^{n},\mathcal{L}^{-1}\left(\phi^{n+1}-\phi^{n}\right)\right\rangle_{\Omega}+\frac{1}{\mathcal{M}\Delta
t}\left\langle\rho^{n+1}-\rho^{n},\mathcal{L}^{-1}\left(\rho^{n+1}-\rho^{n}\right)\right\rangle_{\Omega}$
$\displaystyle+\left\langle\delta_{\phi}E_{c}\left(\phi^{n+1},\rho^{n+1}\right)-\delta_{\phi}E_{e}\left(\phi^{n},\rho^{n}\right),\phi^{n+1}-\phi^{n}\right\rangle_{\Omega}$
$\displaystyle+\left\langle\delta_{\rho}E_{c}\left(\phi^{n+1},\rho^{n+1}\right)-\delta_{\rho}E_{e}\left(\phi^{n},\rho^{n}\right),\rho^{n+1}-\rho^{n}\right\rangle_{\Omega}$
$\displaystyle\geq$ $\displaystyle
E\left(\phi^{n+1},\rho^{n+1}\right)-E\left(\phi^{n},\rho^{n}\right).$ (3.42)

Hence that

$E\left(\phi^{n+1},\rho^{n+1}\right)\leq E\left(\phi^{n},\rho^{n}\right).$

This completes the proof. ∎

###  3.5 Optimal rate convergence analysis

Let $\Phi$ and $\Psi$ be the exact solution for the binary fluid-surfactant
system $\eqref{bfs}$. With initial data with sufficient regularity, we could
assume that the exact solution has regularity of class $\mathcal{R}_{1}$ and
$\mathcal{R}_{2}$:

$\displaystyle\Phi\in\mathcal{R}_{1}:=H^{2}\left(0,T;C_{per}\left(\Omega\right)\right)\cap
L^{\infty}\left(0,T;C_{per}^{8}\left(\Omega\right)\right),$ (3.43a)
$\displaystyle\Psi\in\mathcal{R}_{2}:=H^{2}\left(0,T;C_{per}\left(\Omega\right)\right)\cap
L^{\infty}\left(0,T;C_{per}^{6}\left(\Omega\right)\right).$ (3.43b)

Define
$\Phi_{N}\left(\cdot,t\right):=\mathcal{P}_{N}\Phi\left(\cdot,t\right)$, and
$\Psi_{N}\left(\cdot,t\right):=\mathcal{P}_{N}\Psi\left(\cdot,t\right)$, the
spatial Fourier projection of the exact solutions into $\mathcal{B}^{K}$, the
space of trigonometric polynomials of degree to and including $K$ (with
$N=2K+1$). The following projection approximation is standard: if $\Phi\in
L^{\infty}\left(0,T;H_{per}^{l}\left(\Omega\right)\right)$ for some
$l\in\mathbb{N}$,

$\left\|\Phi_{N}-\Phi\right\|_{L^{\infty}\left(0,T;H^{k}\left(\Omega\right)\right)}\leq
Ch^{l-k}\left\|\Phi\right\|_{L^{\infty}\left(0,T;H^{l}\left(\Omega\right)\right)},\quad
0\leq k\leq l.$ (3.44)

By $\Phi_{N}^{m},\Phi^{m}$, we denote $\Phi_{N}\left(\cdot,t_{m}\right)$ and
$\Phi\left(\cdot,t_{m}\right)$, respectively, with $t_{m}=m\cdot\Delta t$.
Since $\Phi_{N}\in\mathcal{B}^{m}$, the mass conservative property is
available at the discrete level:

$\overline{\Phi_{N}^{m}}=\frac{1}{|\Omega|}\int_{\Omega}\Phi_{N}\left(\cdot,t_{m}\right)\mathrm{d}\bm{x}=\frac{1}{|\Omega|}\int_{\Omega}\Phi_{N}\left(\cdot,t_{m+1}\right)\mathrm{d}\bm{x}=\overline{\Phi_{N}^{m+1}},\quad
m\in\mathbb{N}.$

We have a similar result about $\Psi$. On the other hand, the solution of
$\eqref{bfs}$ is also mass conservative at the discrete level:

$\overline{\phi^{m}}=\overline{\phi^{m+1}},\quad\overline{\rho^{m}}=\overline{\rho^{m+1}},\quad
m\in\mathbb{N}.$ (3.45)

As indicated before, we use the mass conservative interpolation for the
initial data: $\phi^{0}=\mathcal{P}_{h}\Phi_{N}\left(\cdot,t=0\right)$ and
$\rho^{0}=\mathcal{P}_{h}\Psi_{N}\left(\cdot,t=0\right)$, that is

$\phi^{0}_{i,j}=\mathcal{P}_{h}(\Phi_{N}(\cdot,t=0))_{i,j}:=\Phi_{N}\left(x_{i},y_{j},t=0\right),\quad\rho^{0}_{i,j}=\mathcal{P}_{h}(\Psi_{N}(\cdot,t=0))_{i,j}:=\Psi_{N}\left(x_{i},y_{j},t=0\right).$
(3.46)

The error grid function is defined as

$\tilde{\phi}^{m}:=\mathcal{P}_{h}\Phi_{N}^{m}-\phi^{m},\quad\tilde{\rho}^{m}:=\mathcal{P}_{h}\Psi_{N}^{m}-\rho^{m},\quad
m\in\mathbb{N},$ (3.47)

in which a similar interpolation formula could be applied to ${\cal P}_{h}$ as
in (3.46). Therefore, it follows that $\overline{\tilde{\phi}^{m}}=0$ and
$\overline{\tilde{\rho}^{m}}=0$, for any $m\in\mathbb{N}$, so that the
discrete norm $\left\|\cdot\right\|_{-1,h}$ is well defined for the numerical
error grid function.

######  Theorem 3.3.

Given initial data $\Phi(\cdot,t=0)\in C_{per}^{8}(\Omega)$ and
$\Psi(\cdot,t=0)\in C_{per}^{6}(\Omega)$, suppose the exact solutin for binary
fluid-surfactant system $\eqref{bfs}$ is of regularity class
$\mathcal{R}=\mathcal{R}_{1}\times\mathcal{R}_{2}$. Then, provided that
$\Delta t$ is sufficiently small, for all positive integers n, such that
$t_{n}\leq T$, we have

$\|\tilde{\phi}^{n+1}\|_{-1,h}+\|\tilde{\rho}^{n+1}\|_{-1,h}+\left(\mathcal{M}\Delta
t\sum_{k=0}^{n+1}\left(\eta^{2}\|\Delta_{h}\tilde{\phi}^{k}\|_{2}^{2}+\varepsilon\|\nabla_{h}\tilde{\phi}^{k}\|_{2}^{2}+\xi\|\nabla_{h}\tilde{\rho}^{k}\|_{2}^{2}\right)\right)^{\frac{1}{2}}\leq
C\left(\Delta t+h^{2}\right),$ (3.48)

where $C>0$ is independent of n, $\Delta t$, and h.

###### Proof.

A carefully consistency analysis indicates the following truncation error
estimate:

$\displaystyle\frac{\Phi^{n+1}_{N}-\Phi^{n}_{N}}{\Delta t}=$
$\displaystyle\mathcal{M}\Delta_{h}\left(\eta^{2}\Delta_{h}^{2}\Phi^{n+1}_{N}-\varepsilon\Delta_{h}\Phi^{n+1}_{N}+\frac{1}{\varepsilon}(\Phi^{n+1}_{N}-\frac{1}{2})^{3}+\alpha\nabla_{h}\cdot\left({\cal
A}(\frac{\Psi^{n+1}_{N}}{{\cal
A}|\nabla_{h}^{\delta}\Phi^{n+1}_{N}|})\nabla_{h}\Phi^{n+1}_{N}\right)\right.$
$\displaystyle\left.-\alpha\left(1+\frac{1}{\delta}\right)\Delta_{h}\Phi^{n+1}_{N}-\frac{1}{4\varepsilon}(\Phi^{n}_{N}-\frac{1}{2})-\frac{\alpha}{\delta}\Delta_{h}\Phi^{n}_{N}\right)+\tau_{\phi},$
(3.49a) $\displaystyle\frac{\Psi^{n+1}_{N}-\Psi^{n}_{N}}{\Delta t}=$
$\displaystyle\mathcal{M}\Delta_{h}\left(-\xi\Delta_{h}\Psi_{N}^{n+1}+\beta
H^{\prime}\left(\Psi^{n+1}_{N}\right)+\sqrt{2}\alpha\Psi^{n+1}_{N}-\alpha{\cal
A}|\nabla_{h}^{\delta}\Phi^{n+1}_{N}|-\alpha\left(\sqrt{2}-1\right)\Psi^{n}_{N}\right)+\tau_{\rho},$
(3.49b)

with $\left\|\tau^{n}\right\|_{-1,h}\leq C\left(\Delta t+h^{2}\right)$.
Observe that we have dropped the operator $\mathcal{P}_{h}$, which should
appear in front of $\Phi_{N}$, for simplicity of presentation.

Subtracting the numerical scheme $\eqref{cs_1st}$ from $\eqref{estimate_1st}$
gives

$\displaystyle\frac{\tilde{\phi}^{n+1}_{N}-\tilde{\phi}^{n}_{N}}{\Delta t}=$
$\displaystyle\mathcal{M}\Delta_{h}\left(\eta^{2}\Delta_{h}^{2}\tilde{\phi}^{n+1}-\varepsilon\Delta_{h}\tilde{\phi}^{n+1}+\frac{1}{\varepsilon}(\Phi^{n+1}_{N}-\frac{1}{2})^{3}-\frac{1}{\varepsilon}(\phi^{n+1}-\frac{1}{2})^{3}\right.$
$\displaystyle\left.+\alpha\nabla_{h}\cdot\left({\cal
A}(\frac{\Psi^{n+1}_{N}}{{\cal
A}|\nabla_{h}^{\delta}\Phi^{n+1}_{N}|})\nabla_{h}\Phi^{n+1}_{N}\right)-\alpha\nabla_{h}\cdot\left({\cal
A}(\frac{\rho^{n+1}}{{\cal
A}|\nabla_{h}^{\delta}\phi^{n+1}|})\nabla_{h}\phi^{n+1}\right)\right.$
$\displaystyle\left.-\alpha(1+\frac{1}{\delta})\Delta_{h}\tilde{\phi}^{n+1}-\frac{1}{4\varepsilon}\tilde{\phi}^{n}-\frac{\alpha}{\delta}\Delta_{h}\tilde{\phi}^{n}\right)+\tau_{\phi},$
(3.50a) $\displaystyle\frac{\tilde{\rho}^{n+1}-\tilde{\rho}^{n}}{\Delta t}=$
$\displaystyle\mathcal{M}\Delta_{h}\left(-\xi\Delta_{h}\tilde{\rho}^{n+1}+\beta
H^{\prime}\left(\Psi^{n+1}_{N}\right)-\beta
H^{\prime}\left(\rho^{n+1}\right)+\sqrt{2}\alpha\tilde{\rho}^{n+1}\right.$
$\displaystyle\left.-\alpha{\cal
A}|\nabla_{h}^{\delta}\Phi^{n+1}_{N}|+\alpha{\cal
A}|\nabla_{h}^{\delta}\phi^{n+1}|-\alpha(\sqrt{2}-1)\tilde{\rho}^{n}\right)+\tau_{\rho}.$
(3.50b)

Since the numerical error function has zero-mean, we see that both
$\left(-\Delta_{h}\right)^{-1}\tilde{\phi}^{m}$ and
$\left(-\Delta_{h}\right)^{-1}\tilde{\rho}^{m}$ are well-defined, for any
$k\geq 0$. Taking a discrete inner product with $\eqref{estimate_phi}$ and
$\eqref{estimate_rho}$ by
$2\left(-\mathcal{M}\Delta_{h}\right)^{-1}\tilde{\phi}^{n+1}$ and
$2\left(-\mathcal{M}\Delta_{h}\right)^{-1}\tilde{\rho}^{n+1}$, respectively,
yields

$\displaystyle\frac{1}{\mathcal{M}\Delta
t}\left(\|\tilde{\phi}^{n+1}\|_{-1,h}^{2}-\|\tilde{\phi}^{n}\|_{-1,h}^{2}+\|\tilde{\phi}^{n+1}-\tilde{\phi}^{n}\|_{-1,h}^{2}+\|\tilde{\rho}^{n+1}\|_{-1,h}^{2}-\|\tilde{\rho}^{n}\|_{-1,h}^{2}+\|\tilde{\rho}^{n+1}-\tilde{\rho}^{n}\|_{-1,h}^{2}\right)$
$\displaystyle+2\eta^{2}\langle\tilde{\phi}^{n+1},\Delta_{h}^{2}\tilde{\phi}^{n+1}\rangle_{\Omega}-2\varepsilon\langle\tilde{\phi}^{n+1},\Delta_{h}\tilde{\phi}^{n+1}\rangle_{\Omega}-2\xi\langle\tilde{\rho}^{n+1},\Delta_{h}\tilde{\rho}^{n+1}\rangle_{\Omega}$
$\displaystyle+\frac{2}{\varepsilon}\left\langle\tilde{\phi}^{n+1},(\Phi^{n+1}_{N}-\frac{1}{2})^{3}-(\phi^{n+1}-\frac{1}{2})^{3}\right\rangle_{\Omega}-2\alpha(1+\frac{1}{\delta})\langle\tilde{\phi}^{n+1},\Delta_{h}\tilde{\phi}^{n+1}\rangle_{\Omega}$
$\displaystyle+2\alpha\left\langle\tilde{\phi}^{n+1},\nabla_{h}\cdot\left({\cal
A}(\frac{\Psi^{n+1}_{N}}{{\cal
A}|\nabla_{h}^{\delta}\Phi^{n+1}_{N}|})\nabla_{h}\Phi^{n+1}_{N}-{\cal
A}(\frac{\rho^{n+1}}{{\cal
A}|\nabla_{h}^{\delta}\phi^{n+1}|})\nabla_{h}\phi^{n+1}\right)\right\rangle_{\Omega}$
$\displaystyle+2\beta\left\langle\tilde{\rho}^{n+1},H^{\prime}\left(\Psi^{n+1}_{N}\right)-H^{\prime}\left(\rho^{n+1}_{N}\right)\right\rangle_{\Omega}+2\sqrt{2}\alpha\left\langle\tilde{\rho}^{n+1},\tilde{\rho}^{n+1}\right\rangle_{\Omega}-2\alpha\left\langle\tilde{\rho}^{n+1},{\cal
A}|\nabla_{h}^{\delta}\Phi^{n+1}_{N}|-{\cal
A}|\nabla_{h}^{\delta}\phi^{n+1}|\right\rangle_{\Omega}$ $\displaystyle=$
$\displaystyle\frac{1}{2\varepsilon}\langle\tilde{\phi}^{n+1},\tilde{\phi}^{n}\rangle_{\Omega}+\frac{2\alpha}{\delta}\langle\tilde{\phi}^{n+1},\Delta_{h}\tilde{\phi}^{n}\rangle_{\Omega}+2\alpha(\sqrt{2}-1)\left\langle\tilde{\rho}^{n+1},\tilde{\rho}^{n}\right\rangle_{\Omega}$
$\displaystyle+\frac{2}{\mathcal{M}}(\langle\tilde{\phi}^{n+1},\tau_{\phi}^{n}\rangle_{-1,h}+\langle\tilde{\rho}^{n+1},\tau_{\rho}^{n}\rangle_{-1,h}).$

The estimate for the terms associated with the surface diffusion is
straightforward:

$\displaystyle\langle\tilde{\phi}^{n+1},\Delta_{h}^{2}\tilde{\phi}^{n+1}\rangle_{\Omega}$
$\displaystyle=\|\Delta_{h}\tilde{\phi}^{n+1}\|_{2}^{2},\quad-\langle\tilde{\phi}^{n+1},\Delta_{h}\tilde{\phi}^{n+1}\rangle_{\Omega}=\|\nabla_{h}\tilde{\phi}^{n+1}\|_{2}^{2},$
$\displaystyle-\langle\tilde{\rho}^{n+1},\Delta_{h}\tilde{\rho}^{n+1}\rangle_{\Omega}$
$\displaystyle=\|\nabla_{h}\tilde{\rho}^{n+1}\|_{2}^{2}.$

For the nonlinear inner product, we have the following result

$2\beta\left\langle\tilde{\rho}^{n+1},H^{\prime}\left(\Psi^{n+1}_{N}\right)-H^{\prime}\left(\rho^{n+1}_{N}\right)\right\rangle_{\Omega}\geq
0,$ (3.51)

due to the fact that the logarithmic function is an increasing function.
Similarly, the convexity of the nonlinear functional
$g_{1}(\phi)=h^{2}\sum_{i,j=1}^{N}(\phi_{i,j}-\frac{1}{2})^{4}$ and
$g_{2}(\phi,\rho)=h^{2}\sum_{i,j=1}^{N}(\rho_{i,j}-{\cal
A}|\nabla_{h}^{\delta}\phi_{i,j}|)^{2}+(\sqrt{2}-1)\|\rho\|_{2}^{2}+\frac{1}{\delta}\|\nabla_{h}\phi\|_{2}^{2}$
leads to the following inequalities:

$\displaystyle\frac{2}{\varepsilon}\left\langle\tilde{\phi}^{n+1},(\Phi^{n+1}_{N}-\frac{1}{2})^{3}-(\phi^{n+1}-\frac{1}{2})^{3}\right\rangle_{\Omega}\geq
0,$ $\displaystyle
2\alpha\left\langle\tilde{\phi}^{n+1},\nabla_{h}\cdot\left({\cal
A}(\frac{\Psi^{n+1}_{N}}{{\cal
A}|\nabla_{h}^{\delta}\Phi^{n+1}_{N}|})\nabla_{h}\Phi^{n+1}_{N}-{\cal
A}(\frac{\rho^{n+1}}{{\cal
A}|\nabla_{h}^{\delta}\phi^{n+1}|})\nabla_{h}\phi^{n+1}\right)\right\rangle_{\Omega}-2\alpha(1+\frac{1}{\delta})\langle\tilde{\phi}^{n+1},\Delta_{h}\tilde{\phi}^{n+1}\rangle_{\Omega}$
$\displaystyle+2\sqrt{2}\alpha\left\langle\tilde{\rho}^{n+1},\tilde{\rho}^{n+1}\right\rangle_{\Omega}-2\alpha\left\langle\tilde{\rho}^{n+1},{\cal
A}|\nabla_{h}^{\delta}\Phi^{n+1}_{N}|-{\cal
A}|\nabla_{h}^{\delta}\phi^{n+1}|\right\rangle_{\Omega}\geq 0.$ (3.52)

Then we arrive at the following estimate:

$\displaystyle\frac{1}{\mathcal{M}\Delta
t}\left(\|\tilde{\phi}^{n+1}\|_{-1,h}^{2}-\|\tilde{\phi}^{n}\|_{-1,h}^{2}+\|\tilde{\phi}^{n+1}-\tilde{\phi}^{n}\|_{-1,h}^{2}+\|\tilde{\rho}^{n+1}\|_{-1,h}^{2}-\|\tilde{\rho}^{n}\|_{-1,h}^{2}+\|\tilde{\rho}^{n+1}-\tilde{\rho}^{n}\|_{-1,h}^{2}\right)$
$\displaystyle+2\eta^{2}\|\Delta_{h}\tilde{\phi}^{n+1}\|_{2}^{2}+2\varepsilon\|\nabla_{h}\tilde{\phi}^{n+1}\|_{2}^{2}+2\xi\|\nabla_{h}\tilde{\rho}^{n+1}\|_{2}^{2}$
$\displaystyle\leq$
$\displaystyle\frac{1}{2\varepsilon}\left\langle\tilde{\phi}^{n+1},\tilde{\phi}^{n}\right\rangle_{\Omega}+\frac{2\alpha}{\delta}\left\langle\tilde{\phi}^{n+1},\Delta_{h}\tilde{\phi}^{n}\right\rangle_{\Omega}+2\alpha\left(\sqrt{2}-1\right)\left\langle\tilde{\rho}^{n+1},\tilde{\rho}^{n}\right\rangle_{\Omega}$
$\displaystyle+\frac{2}{\mathcal{M}}\left\langle\tilde{\phi}^{n+1},\tau_{\phi}^{n}\right\rangle_{-1,h}+\frac{2}{\mathcal{M}}\left\langle\tilde{\rho}^{n+1},\tau_{\rho}^{n}\right\rangle_{-1,h}.$

Meanwhile, for the inner product associated with the concave part, the
following inequalities could be derived:

$\displaystyle\frac{1}{2\varepsilon}\left\langle\tilde{\phi}^{n+1},\tilde{\phi}^{n}\right\rangle_{\Omega}$
$\displaystyle\leq\frac{\varepsilon}{2}\|\nabla_{h}\tilde{\phi}^{n+1}\|_{2}^{2}+\frac{1}{8\varepsilon^{3}}\|\tilde{\phi}^{n}\|_{-1,h}^{2},$
(3.53a) $\displaystyle
2\left(\sqrt{2}-1\right)\alpha\left\langle\tilde{\rho}^{n+1},\tilde{\rho}^{n}\right\rangle_{\Omega}$
$\displaystyle\leq\xi\|\nabla_{h}\tilde{\rho}^{n+1}\|_{2}^{2}+\frac{\alpha^{2}\left(\sqrt{2}-1\right)^{2}}{\xi}\|\tilde{\rho}^{n}\|_{-1,h}^{2},$
(3.53b)
$\displaystyle\frac{2}{\mathcal{M}}\left\langle\tilde{\phi}^{n+1},\tau_{\phi}^{n}\right\rangle_{-1,h}$
$\displaystyle\leq\frac{1}{\mathcal{M}}\|\tilde{\phi}^{n+1}\|_{-1,h}^{2}+\frac{1}{\mathcal{M}}\|\tau_{\phi}^{n}\|_{-1,h}^{2},$
(3.53c)
$\displaystyle\frac{2}{\mathcal{M}}\left\langle\tilde{\rho}^{n+1},\tau_{\rho}^{n}\right\rangle_{-1,h}$
$\displaystyle\leq\frac{1}{\mathcal{M}}\|\tilde{\rho}^{n+1}\|_{-1,h}^{2}+\frac{1}{\mathcal{M}}\|\tau_{\rho}^{n}\|_{-1,h}^{2},$
(3.53d)
$\displaystyle-\frac{2\alpha}{\delta}\left\langle\tilde{\phi}^{n+1},\Delta_{h}\tilde{\phi}^{n}\right\rangle_{\Omega}$
$\displaystyle=\frac{2\alpha}{\delta}\left\langle\nabla_{h}\tilde{\phi}^{n+1},\nabla_{h}\tilde{\phi}^{n}\right\rangle_{\Omega}$
$\displaystyle\leq\frac{\varepsilon}{2}\|\nabla_{h}\tilde{\phi}^{n+1}\|_{2}^{2}+\frac{2\alpha^{2}}{\varepsilon\delta^{2}}\|\nabla_{h}\tilde{\phi}^{n}\|_{2}^{2}$
$\displaystyle\leq\frac{\varepsilon}{2}\|\nabla_{h}\tilde{\phi}^{n+1}\|_{2}^{2}+\frac{\alpha^{2}}{\varepsilon\delta^{2}}\|\tilde{\phi}^{n}\|_{-1}^{\frac{2}{3}}\|\Delta_{h}\tilde{\phi}^{n}\|_{2}^{\frac{4}{3}}$
$\displaystyle\leq\frac{\varepsilon}{2}\|\nabla_{h}\tilde{\phi}^{n+1}\|_{2}^{2}+\frac{8\alpha^{6}}{\eta^{4}\varepsilon^{3}\delta^{6}}\|\tilde{\phi}^{n}\|_{-1,h}^{2}+\eta^{2}\|\Delta_{h}\tilde{\phi}^{n}\|_{2}^{2}.$
(3.53e)

Therefore, we obtain

$\displaystyle\frac{1}{\mathcal{M}\Delta
t}\left(\|\tilde{\phi}^{n+1}\|_{-1,h}^{2}-\|\tilde{\phi}^{n}\|_{-1,h}^{2}+\|\tilde{\phi}^{n+1}-\tilde{\phi}^{n}\|_{-1,h}^{2}+\|\tilde{\rho}^{n+1}\|_{-1,h}^{2}-\|\tilde{\rho}^{n}\|_{-1,h}^{2}+\|\tilde{\rho}^{n+1}-\tilde{\rho}^{n}\|_{-1,h}^{2}\right)$
$\displaystyle+2\eta^{2}\|\Delta_{h}\tilde{\phi}^{n+1}\|_{2}^{2}+2\varepsilon\|\nabla_{h}\tilde{\phi}^{n+1}\|_{2}^{2}+2\xi\|\nabla_{h}\tilde{\rho}^{n+1}\|_{2}^{2}$
$\displaystyle\leq$
$\displaystyle\frac{\varepsilon}{2}\|\nabla_{h}\tilde{\phi}^{n+1}\|_{2}^{2}+\frac{1}{8\varepsilon^{2}}\|\tilde{\phi}^{n}\|_{-1,h}^{2}+\xi\|\nabla_{h}\tilde{\rho}^{n+1}\|_{2}^{2}+\frac{\alpha^{2}\left(\sqrt{2}-1\right)^{2}}{\xi}\|\tilde{\rho}^{n}\|_{-1,h}^{2}+\frac{1}{\mathcal{M}}\|\tilde{\phi}^{n+1}\|_{-1,h}^{2}+\frac{1}{\mathcal{M}}\|\tau_{\phi}^{n}\|_{-1,h}^{2}$
$\displaystyle+\frac{1}{\mathcal{M}}\|\tilde{\rho}^{n+1}\|_{-1,h}^{2}+\frac{1}{\mathcal{M}}\|\tau_{\rho}^{n}\|_{-1,h}^{2}+\frac{\varepsilon}{2}\|\nabla_{h}\tilde{\phi}^{n+1}\|_{2}^{2}+\frac{8\alpha^{6}}{\eta^{4}\varepsilon^{3}\delta^{6}}\|\tilde{\phi}^{n}\|_{-1,h}^{2}+\eta^{2}\|\Delta_{h}\tilde{\phi}^{n}\|_{2}^{2},$
(3.54)

which in turn gives

$\displaystyle\frac{1}{\mathcal{M}\Delta
t}\left(\|\tilde{\phi}^{n+1}\|_{-1,h}^{2}-\|\tilde{\phi}^{n}\|_{-1,h}^{2}+\|\tilde{\rho}^{n+1}\|_{-1,h}^{2}-\|\tilde{\rho}^{n}\|_{-1,h}^{2}\right)$
$\displaystyle+2\eta^{2}\|\Delta_{h}\tilde{\phi}^{n+1}\|_{2}^{2}+\varepsilon\|\nabla_{h}\tilde{\phi}^{n+1}\|_{2}^{2}+\xi\|\nabla_{h}\tilde{\rho}^{n+1}\|_{2}^{2}$
$\displaystyle\leq$
$\displaystyle\left(\frac{1}{8\varepsilon^{2}}+\frac{8\alpha^{6}}{\eta^{4}\varepsilon^{3}\delta^{6}}\right)\|\tilde{\phi}^{n}\|_{-1,h}^{2}+\frac{\alpha^{2}\left(\sqrt{2}-1\right)^{2}}{\xi}\|\tilde{\rho}^{n}\|_{-1,h}^{2}+\eta^{2}\|\Delta_{h}\tilde{\phi}^{n}\|_{2}^{2}$
$\displaystyle+\frac{1}{\mathcal{M}}\|\tilde{\phi}^{n+1}\|_{-1,h}^{2}+\frac{1}{\mathcal{M}}\|\tau_{\phi}^{n}\|_{-1,h}^{2}+\frac{1}{\mathcal{M}}\|\tilde{\rho}^{n+1}\|_{-1,h}^{2}+\frac{1}{\mathcal{M}}\|\tau_{\rho}^{n}\|_{-1,h}^{2}.$
(3.55)

Finally, an application of a discrete Gronwall inequality results in the
desired convergence estimate:

$\|\tilde{\phi}^{n+1}\|_{-1,h}+\|\tilde{\rho}^{n+1}\|_{-1,h}+\left(\mathcal{M}\Delta
t\sum_{k=0}^{n+1}\left(\eta^{2}\|\Delta_{h}\tilde{\phi}^{k}\|_{2}^{2}+\varepsilon\|\nabla_{h}\tilde{\phi}^{k}\|_{2}^{2}+\xi\|\nabla_{h}\tilde{\rho}^{k}\|_{2}^{2}\right)\right)^{\frac{1}{2}}\leq
C\left(\Delta t+h^{2}\right),$ (3.56)

where $C>0$ is independent of $\Delta t,h$ and $n$. This completes the proof
of Theorem 3.3. ∎

######  Remark 3.1.

In the application of the discrete Gronwall inequality, we see that the growth
constants for $\|\tilde{\phi}^{n}\|_{-1,h}^{2}$ and
$\|\tilde{\rho}^{n}\|_{-1,h}^{2}$ terms, given by
$\frac{1}{8\varepsilon^{2}}+\frac{8\alpha^{6}}{\eta^{4}\varepsilon^{3}\delta^{6}}$
and $\frac{\alpha^{2}\left(\sqrt{2}-1\right)^{2}}{\xi}$, respectively, depend
singularly on $\varepsilon$, $\eta$, $\delta$ and $\xi$. In turn, it would be
reasonable to require that

$\Big{(}\frac{1}{8\varepsilon^{2}}+\frac{8\alpha^{6}}{\eta^{4}\varepsilon^{3}\delta^{6}}\Big{)}{\cal
M}\Delta t\leq 1,\quad\frac{\alpha^{2}\left(\sqrt{2}-1\right)^{2}}{\xi}{\cal
M}\Delta t\leq 1,$

so that a singular convergence constant is avoided. In other words, the time
step size $\Delta t$ should be bounded by a given constant, dependent on
$\varepsilon$, $\eta$, $\delta$ and $\xi$, to present a singular convergence
constant at a theoretical level; this requirement refers to the condition that
“provided that $\Delta t$ is sufficiently small” in the statement of Theorem
3.3. Meanwhile, such a requirement is only associated with a theoretical
analysis, and this requirement may not be necessary in the practical
computations to preserve a numerical convergence.

######  Remark 3.2.

As the regularization parameter $\delta\to 0$, the positivity-preserving
property and the energy stability estimates, as established in Theorems 3.1
and 3.2, are still valid. In fact, these two theoretical properties are
available even with $\delta=0$. On the other hand, the optimal rate
convergence estimate, as established in Theorem 3.3, is only available for a
fixed $\delta>0$, due to the singularly-$\delta$-dependent convergence
constant appearing in (3.55). In other words, the convergence constant in
Theorem 3.3 depends singularly on $\delta$, and such a convergence estimate
would not be theoretically justified as $\delta\to 0$, although the numerical
convergence has also been verified in various numerical experiments.

######  Remark 3.3.

In the proposed numerical scheme (3.16), we take $M_{1}=M_{2}={\cal M}$, and
$M(\rho)\equiv 1$, for simplicity of presentation. In case of a
$\rho$-dependent mobility function $M(\rho)$, the positivity-preserving
property and energy stability are still valid, as long as $M(\rho)>0$ is
available at a point-wise level. Meanwhile, the corresponding convergence
analysis and error estimate are expected to face certain theoretical
difficulties in the case of a non-constant mobility function, due to the
highly nonlinear and singular nature of the chemical potential. The
theoretical justification of this convergence analysis will be left to the
future works, and some techniques of rough error estimate and refined error
estimate, as reported in a recent work [34] to analyze the non-constant-
mobility Poisson-Nernst-Planck system, may have to be applied in this future
work.

##  4 Numerical experiments

In this section, we preform a few two-dimensional numerical simulations using
the proposed scheme $\eqref{cs_1st}$. The mass conservation, energy decay,
positivity of the numerical solution, as well as the numerical accuracy, will
be demonstrated in these computations. To achieve this goal, we will present
two numerical examples with different initial conditions.

###  4.1 Accuracy test

Here, we take the domain $\Omega=(0,8)^{2}$, and choose the parameters as
follows

$\displaystyle\varepsilon=0.05,\quad\alpha=0.001,\quad\beta=0.02,\quad\delta=0.001,$
$\displaystyle\eta=0.05,\quad\xi=0.05,\quad M_{1}=0.01,\quad M_{2}=0.01.$

The initial data are set as

$\left\\{\begin{aligned} &\phi_{0}(x,y)=0.5+0.2\cos\frac{4\pi
x}{8}\cos\frac{4\pi y}{8},\\\ &\rho_{0}(x,y)=0.5+0.2\sin\frac{4\pi
x}{8}\sin\frac{4\pi y}{8}.\end{aligned}\right.$ (4.1)

It is obvious that the initial data are subject to periodic boundary
condition. This example is designed to study the numerical accuracy in time
and space. In order to test the first order convergence rate in time and
second order convergence rate in space, we use a linear refinement pate, i.e.
$\Delta t=Ch^{2}$, $C=0.01$. The global error is expected to be $O(\Delta
t)+O(h^{2})=O(h^{2})$ under the discrete $L^{2}$ norm. Since an exact solution
is not available, we compute the Cauchy difference instead of directly
calculating the numerical error, which is defined as
$\delta_{u}=u_{h_{f}}-\mathcal{I}_{c}^{f}(u_{h_{c}})$, where
$\mathcal{I}_{c}^{f}$ is a bilinear interpolation operator. This requires a
relatively coarse solution, parametrized by $h_{c}$, and a relatively fine
solution, parametrized by $h_{f}$, where $h_{c}=2h_{f}$, at the same final
time. The discrete $L^{2}$ norms of Cauchy difference and the convergence
rates are displayed in Table 1. These results confirm the expected convergence
rate.

Table 1: The discrete $L^{2}$ error and convergence rate at $t=0.1$ with initial data $\eqref{CS:test}$ and the given parameters. Grid sizes | Error$(\phi)$ | Rate | Error$(\rho)$ | Rate
---|---|---|---|---
$16\times 16$ | 1.93E–01 | – | 1.88E–01 | –
$32\times 32$ | 5.07E–02 | 1.93 | 4.86E–02 | 1.95
$64\times 64$ | 1.28E–02 | 1.98 | 1.23E–02 | 1.99
$128\times 128$ | 3.21E–03 | 2.00 | 3.07E–03 | 2.00
$256\times 256$ | 8.04E–04 | 2.00 | 7.68E–04 | 2.00

###  4.2 Spinodal decomposition

In this example, we study the phase separation phenomenon, so called spinodal
decomposition. Usually we describe this process as a thermal quench, which is
considered that an initially homogeneous mixture is thrust into a two-phase
region. In this case, the spinodal decomposition occurs and leads the system
from the homogeneous to two-phase state. We take the domain as
$\Omega=(0,2\pi)^{2}$. The initial data are given by

$\left\\{\begin{aligned} &\phi_{0}(x,y)=0.4+0.1\mbox{rand}(x,y),\\\
&\rho_{0}(x,y)=0.4+0.1\mbox{rand}(x,y),\end{aligned}\right.$ (4.2)

where rand$(x,y)$ is a random number in $[-1,1]$ and has zero mean. The
parameters are chosen as follows

$\varepsilon=0.02,~{}\alpha=0.02,~{}\beta=0.02,~{}\eta=0.02,~{}\delta=0.01,~{}\xi=0.02,~{}\mathcal{M}=0.01.$
(4.3)

From Figures 1 to 6, we display the snapshots of coarsening dynamics.
Initially, the two fluids are well mixed, and they sooner start to decompose
and accumulate. We observe that a relatively high value of the concentration
variable $\rho$ gathers at the interface between the two different fluids. A
monotone decay evolution of the physical energy is illustrated in Figure 7.

(a) $\phi(t=0)$

(b) $\rho(t=0)$

Figure 1: Snapshots of the phase variables $\phi$ and $\rho$, taken at $t=0$
for Example 4.2.

(a) $\phi(t=0.5)$

(b) $\rho(t=0.5)$

Figure 2: Snapshots of the phase variables $\phi$ and $\rho$, taken at $t=0.5$
for Example 4.2.

(a) $\phi(t=2)$

(b) $\rho(t=2)$

Figure 3: Snapshots of the phase variables $\phi$ and $\rho$, taken at $t=2$
for Example 4.2.

(a) $\phi(t=10)$

(b) $\rho(t=10)$

Figure 4: Snapshots of the phase variables $\phi$ and $\rho$, taken at $t=10$
for Example 4.2.

(a) $\phi(t=40)$

(b) $\rho(t=40)$

Figure 5: Snapshots of the phase variables $\phi$ and $\rho$, taken at $t=40$
for Example 4.2.

(a) $\phi(t=700)$

(b) $\rho(t=700)$

Figure 6: Snapshots of the phase variables $\phi$ and $\rho$, taken at $t=700$
for Example 4.2. Figure 7: Time evolution of the free energy functional for
Example 4.2, which shows a monotone decay.

##  5 Conclusions

In this paper, we propose and analyze a first order (in time) accurate, convex
splitting scheme for the binary fluid-surfactant phase field model. The multi-
phase structure and the singularity associated with the 1-Laplacian part makes
the whole system very challenging, at both the theoretical and numerical
levels. To overcome this subtle difficulty, we make an observation for a non-
standard convex-concave decomposition of the free energy. In addition, the
singular nature of the logarithmic function around the limit values prevents
the numerical solution approaching these limit values, so that the positivity
property is preserved for the numerical scheme. As a result, the convex
structure of the implicit part guarantees the unique solvability and energy
stability of the proposed numerical scheme. Furthermore, an optimal rate
convergence analysis is carefully derived, which is the first such result in
this area. Some numerical experiments are performed to validate the accuracy
and energy stability of the proposed scheme.

##  6 Acknowledgements

The authors greatly appreciate many helpful discussions with Professor Hui
Zhang, in particular for his insightful suggestions and comments. This work is
supported in part by the grants NSF DMS-2012669 (C. Wang), NSFC-11871105, the
Science Challenge Project TZ2018002 and the Fundamental Research Funds for the
Central Universities (Z. Zhang).

## References

  * [1] F. Boyer and S. Minjeaud.  Numerical schemes for a three component Cahn-Hilliard model.  ESAIM Mathematical Modeling and Numerical Analysis, 45(4):697–738, 2011.
  * [2] A. B. Branger and D. M. Eckmann.  Accelerated arteriolar gas embolism reabsorption by an exogenous surfactant.  Anesthesiology, 96(4):971–979, 2002.
  * [3] W. Chen, Y. Liu, C. Wang, and S. M. Wise.  Convergence analysis of a fully discrete finite difference scheme for the Cahn-Hilliard-Hele-Shaw equation.  Mathematics of Computation, 85(301):2231–2257, 2015.
  * [4] W. Chen, C. Wang, X. Wang, and S. M. Wise.  A linear iteration algorithm for a second-order energy stable scheme for a thin film model without slope selection.  Journal of Scientific Computing, 59:574–601, 2014.
  * [5] W. Chen, C. Wang, X. Wang, and S. M. Wise.  Positivity-preserving, energy stable numerical schemes for the Cahn-Hilliard equation with logarithmic potential.  Journal of Computational Physics: X, 3:1–29, 2019.
  * [6] Q. Cheng, X. Yang, and J. Shen.  Efficient and accurate numerical schemes for a hydro-dynamically coupled phase field diblock copolymer model.  Journal of Computational Physics, 341(15):44–60, 2017.
  * [7] A. E. Diegel and S. W. Walker.  A finite element method for a phase field model of nematic liquid crystal droplets.  Communications in Computational Physics, 25(1):155–188, 2017.
  * [8] A. E. Diegel, C. Wang, X. Wang, and S. M. Wise.  Convergence analysis and error estimates for a second order accurate finite element method for the Cahn-Hilliard-Navier-Stokes system.  Numerische Mathematik, 137:495–534, 2017.
  * [9] L. Dong, C. Wang, S. Wise, and Z. Zhang.  A positivity-preserving, energy stable scheme for a ternary Cahn-Hilliard system with the singular interfacial parameters.  Journal of Computational Physics, 2021.  Submitted and in review.
  * [10] L. Dong, C. Wang, H. Zhang, and Z. Zhang.  A positivity-preserving, energy stable and convergent numerical scheme for the Cahn-Hilliard equation with a Flory-Huggins-Degennes energy.  Communications in Computational Physics, 17(4):921–939, 2019.
  * [11] L. Dong, C. Wang, H. Zhang, and Z. Zhang.  A positivity-preserving second-order BDF scheme for the Cahn-Hilliard equation with variable interfacial parameters.  Communications in Computational Physics, 28(3):967–998, 2020.
  * [12] Q. Du, L. Ju, X. Li, and Z. Qiao.  Maximum bound principles for a class of semi-linear parabolic equations and exponential time differencing schemes.  SIAM Review, 2021.  Accepted and in press: https://arxiv.org/abs/2015.114.65.
  * [13] D. Eyre.  Unconditionally gradient stable time marching the Cahn-Hilliard equation.  Materials Research Society, Warrendale, PA, USA, 53:1686–1712, 1998.
  * [14] W. Feng, C. Wang, S. Wise, and Z. Zhang.  A second-order energy stable backward differentiation formula method for the epitaxial thin film equation with slope selection.  Numerical Methods for Partial Differential Equations, 34(6):1975–2007, 2018.
  * [15] X. Feng and A. Prohl.  Numerical analysis of the Allen-Cahn equation and approximation for mean curvature flows.  Numerische Mathematik, 94(1):33–65, 2003.
  * [16] G. J. Fix.  Phase field models for free boundary problems.  In A. Fasano and M. Primicerio, editors, Free boundary problems: theory and applications, volume 2, pages 580–589. Pitman, 1983.
  * [17] I. Fonseca, M. Morini, and V. Slastikov.  Surfactants in foam stability: a phase-field model.  Archive for Rational Mechanics and Analysis, 183(3):411–456, 2007\.
  * [18] Lord Rayleigh F.R.S.  On the theory of surface forces.II.Compressible fluids.  Philosophical Magazine Series 5, 33(201):209–220, 1892.
  * [19] J. W. Gibbs.  On the equilibrium of heterogeneous substances.  Transactions of the Connecticut Academy of Arts and Sciences, 3:108–248 and 343–524, 1874-1878.
  * [20] S. Gu, H. Zhang, and Z. Zhang.  An energy-stable finite-difference scheme for the binary fluid-surfactant system.  Journal of Computational Physics, 270(3):416–431, 2014.
  * [21] F. Guillén-González and G. Tierra.  On linear schemes for a Cahn-Hilliard diffuse interface model.  Journal of Computational Physics, 234(2):140–171, 2013.
  * [22] J. Guo, C. Wang, and S. Wise.  An $H^{2}$ convergence of a second-order convex-splitting, finite difference scheme for the three-dimensional Cahn-Hilliard equation.  Communications in Mathematical Sciences, 14:489–515, 2016.
  * [23] D. Han, A. Brylev, X. Yang, and Z. Tan.  Numerical analysis of second order, fully discrete energy stable schemes for phase field models of two-phase incompressible flows.  Journal of Scientific Computing, 70(3):965–989, 2017.
  * [24] A. J. James and J. Lowengrub.  A surfactant-conserving volume-of-fluid method for interfacial flows with insoluble surfactant.  Journal of Computational Physics, 201(2):685–722, 2004.
  * [25] G. Ji.  A BDF2 energy-stable scheme for a general tensor based model of liquid crystals.  East Asian Journal of Applied Mathematics, 10(1):57–71, 2020.
  * [26] S. Khatri and A. K. Tornberg.  An embedded boundary method for soluble surfactants with interface tracking for two-phase flows.  Journal of Computational Physics, 256(1):768–790, 2014.
  * [27] S. Komura and H. Kodama.  Two-order-parameter model for an oil-water-surfactant system.  Physical Review E, 55(2):1722–1727, 1997.
  * [28] J. S. Langer.  Models of pattern formation in first-order phase transitions.  In G. Grinstein and G. Mazenko, editors, Directions in Condensed Matter Physics. Directions in Condensed Matter Physics. Series: Series on Directions in Condensed Matter Physics. Series on Directions in Condensed Matter Physics, volume 1, pages 165–186. World Scientific, 1986.
  * [29] M. Laradji, H. Guo, M. Grant, and M. J. Zuckermann.  The effect of surfactants on the dynamics of phase separation.  Journal of Physics Condensed Matter, 4(4):6715–6728, 1992.
  * [30] D. Li and Z. Qiao.  On second order semi-implicit Fourier spectral methods for 2D Cahn-Hilliard equations.  Journal of Scientific Computing, 70:301–341, 2017.
  * [31] D. Li and Z. Qiao.  On the stabilization size of semi-implicit Fourier-spectral methods for 3D Cahn-Hilliard equations.  Communications in Mathematical Sciences, 15:1489–1506, 2017.
  * [32] D. Li, Z. Qiao, and T. Tang.  Characterizing the stabilization size for semi-implicit Fourier-spectral method to phase field equations.  SIAM Journal on Numerical Analysis, 54:1653–1681, 2016.
  * [33] X. Li, J. Shen, and H. Rui.  Energy stability and convergence of SAV Block-Centered finite difference method for gradient flows.  Mathematics of Computation, 88(319):2047–2068, 2019.
  * [34] C. Liu, C. Wang, S. Wise, X. Yue, and S. Zhou.  A positivity-preserving, energy stable and convergent numerical scheme for the Poisson-Nernst-Planck system.  Mathematics of Computation, 2021.  Accepted and in press.
  * [35] H. Liu and Y. Zhang.  Phase-field modeling droplet dynamics with soluble surfactants.  Journal of Computational Physics, 229(24):9166–9187, 2010.
  * [36] J.-G. Liu, R. Samelson, and C. Wang.  Global weak solution of the planetary geostrophic equations with inviscid geostrophic balance.  Applicable Analysis, 85:593–606, 2007.
  * [37] D. Myers.  Fluid surfaces an interfaces.  In Surfactant science and technology, pages 80–106. Wiley-Interscience, 2005.
  * [38] Y. Qian, C. Wang, and S. Zhou.  A positive and energy stable numerical scheme for the Poisson-Nernst-Planck-Cahn-Hilliard equations with steric interactions.  Journal of Computational Physics, 426:109908, 2021.
  * [39] Y. Qin, Z. Xu, H. Zhang, and Z. Zhang.  Fully decoupled, linear and unconditionally energy stable schemes for the binary fluid-surfactant model.  Communications in Computational Physics, 28(4):1389–1414, 2020\.
  * [40] R. Samelson, R. Temam, C. Wang, and S. Wang.  A fourth order numerical method for the planetary geostrophic equations with inviscid geostrophic balance.  Numerische Mathematik, 107:669–705, 2007.
  * [41] R. Samelson, R. Temam, and S. Wang.  Some mathematical properties of the planetary geostrophic equations for large scale ocean circulation.  Applicable Analysis, 70:147–173, 1998.
  * [42] R. Samelson, R. Temam, and S. Wang.  Remarks on the planetary geostrophic model of grye scale ocean circulation.  Differential and Integral Equations, 13:1–14, 2000.
  * [43] J. Shen, C. Wang, X. Wang, and S. M. Wise.  Second-order convex splitting schemes for gradient flows with Ehrlich-Schwoebel type energy: application to thin film epitaxy.  SIAM Journal on Numerical Analysis, 50(1):105–125, 2012.
  * [44] J. Shen and J. Xu.  Convergence and error analysis for the scalar auxiliary variable (SAV) schemes to gradient flows.  SIAM Journal on Numerical Analysis, 56(5):2895–2912, 2018.
  * [45] J. Shen, J. Xu, and J. Yang.  The scalar auxiliary variable (SAV) approach for gradient flows.  Journal of Computational Physics, 353(15):407–416, 2018.
  * [46] J. Shen, J. Xu, and J. Yang.  A new class of efficient and robust energy stable schemes for gradient flows.  SIAM Review, 61(3):474–506, 2019.
  * [47] J. Shen and X. Yang.  Numerical approximations of Allen-Cahn and Cahn-Hilliard equations.  Discrete and Continuous Dynamical Systems, 28(4):1669–1691, 2010\.
  * [48] J. Stefan.  Über die theorie der eisbildung.  Monatshefte Für Mathematik Und Physik, 1(1):1–6, 1890.
  * [49] K. E. Teigen, P. Song, J. Lowengrub, and A. Voigt.  A diffuse-interface method for two-phase flows with soluble surfactants.  Journal of Computational Physics, 230(2):375–393, 2011.
  * [50] C. Teng, I. Chern, and M. Lai.  Simulating binary fluid-surfactant dynamics by a phase field model.  Discrete and Continuous Dynamical Systems - Series B, 17(4):1289–1307, 2012.
  * [51] T. Teramoto and F. Yonezawa.  Droplet growth dynamics in a water/oil/surfactant system.  Journal of Colloid and Interface Science, 235(2):329–333, 2001\.
  * [52] O. Theissen and G. Gompper.  Lattice-Boltzmann study of spontaneous emulsification.  The European Physical Journal B - Condensed Matter and Complex Systems, 11(1):91–100, 1999.
  * [53] R. G. M. van der Sman and S. van der Graaf.  Diffuse interface model of surfactant adsorption onto flat and droplet interfaces.  Rheologica Acta, 46(1):3–11, 2006.
  * [54] J. D. van der Waals.  The thermodynamic theory of capillarity under the hypothesis of a continuous variation of density.  Journal of Statistical Physics, 20(2):200–244, 1979.
  * [55] C. Wang and S. M. Wise.  Global smooth solutions of the three-dimensional modified phase field crystal equation.  Methods and Applications of Analysis, 17(2):191–212, 2010.
  * [56] C. Wang and S. M. Wise.  An energy stable and convergent finite-difference scheme for the modified phase field crystal equation.  SIAM Journal of Numerical Analysis, 49(3):945–969, 2011.
  * [57] S.M. Wise, C. Wang, and J. Lowengrub.  An energy-stable and convergent finite-stable scheme for the phase field crystal equation.  SIAM Journal of Numerical Analysis, 47(3):2269–2288, 2009.
  * [58] X. Wu, G. J. van Zwieten, and K. G. van der Zee.  Stabilized second-order convex splitting schemes for Cahn-Hilliard models with application to diffuse-interface tumor-growth models.  International Journal for Numerical Methods in Biomedical Engineering, 30(2):180–203, 2017.
  * [59] Y. Yan, W. Chen, C. Wang, and S. M. Wise.  A second-order energy stable BDF numerical scheme for the Cahn-Hilliard equation.  Communications in Computational Physics, 23(2):572–602, 2018.
  * [60] X. Yang.  Linear, first and second-order, unconditionally energy stable numerical schemes for the phase field model of homopolymer blends.  Journal of Computational Physics, 327(15):294–316, 2016.
  * [61] X. Yang.  Numerical approximations for the Cahn-Hilliard phase field model of the binary fluid-surfactant system.  Journal of Scientific Computing, 74(3):1533–1553, 2018.
  * [62] X. Yang and L. Ju.  Efficient linear schemes with unconditional energy stability for the phase field elastic bending energy model.  Computer Methods in Applied Mechanics and Engineering, 315(1):691–712, 2017.
  * [63] X. Yang and L. Ju.  Linear and unconditionally energy stable schemes for the binary fluid-surfactant phase field model.  Computer Methods in Applied Mechanics and Engineering, 318(1):1005–1029, 2017.
  * [64] X. Yang, J. Zhao, and Q. Wang.  Numerical approximations for the molecular beam epitaxial growth model based on the invariant energy quadratization method.  Journal of Computational Physics, 333(15):104–127, 2017.
  * [65] X. Yang, J. Zhao, Q. Wang, and J. Shen.  Numerical approximations for a three-component Cahn-Hilliard phase-field model based on the invariant energy quadratization method.  Mathematical Models and Methods in Applied Sciences, 27(11):1993–2030, 2017.
  * [66] C. Zhang, H. Li, X. Zhang, and L. Ju.  Linear and unconditionally energy stable schemes for the multi-component two-phase diffuse interface model with Peng-Robinson equation of state.  Communications in Computational Physics, 26(4):1071–1097, 2019\.
  * [67] J. Zhang, C. Wang, S.M. Wise, and Z. Zhang.  Structure-preserving, energy stable numerical schemes for a droplet liquid film coarsening model.  SIAM Journal on Scientific Computing, 2021.  Accepted and in press.
  * [68] J. Zhao, Q. Wang, and X. Yang.  Numerical approximations for a phase field dendritic crystal growth model based on the invariant energy quadratization approach.  International Journal for Numerical Methods in Engineering, 110(3):279–300, 2017.
  * [69] G. Zhu, J. Kou, S. Sun, J. Yao, and A. Li.  Decoupled, energy stable schemes for a phase-field surfactant model.  Computer Physics Communications, 233:67–77, 2018.