# Structure-preserving, energy stable numerical schemes for a liquid thin film
coarsening model

Juan Zhang School of Mathematical Sciences, Beijing Normal University, Beijing
100875, P.R. China<EMAIL_ADDRESS>Cheng Wang Department of
Mathematics, The University of Massachusetts, North Dartmouth, MA 02747, USA
(Corresponding Author<EMAIL_ADDRESS>Steven M. Wise Department of
Mathematics, The University of Tennessee, Knoxville, TN 37996, USA
<EMAIL_ADDRESS>Zhengru Zhang Laboratory of Mathematics and Complex Systems,
Beijing Normal University, Beijing 100875, P.R. China<EMAIL_ADDRESS>
###### Abstract

In this paper, two finite difference numerical schemes are proposed and
analyzed for the droplet liquid film model, with a singular Leonard-Jones
energy potential involved. Both first and second order accurate temporal
algorithms are considered. In the first order scheme, the convex potential and
the surface diffusion terms are implicitly, while the concave potential term
is updated explicitly. Furthermore, we provide a theoretical justification
that this numerical algorithm has a unique solution, such that the positivity
is always preserved for the phase variable at a point-wise level, so that a
singularity is avoided in the scheme. In fact, the singular nature of the
Leonard-Jones potential term around the value of 0 prevents the numerical
solution reaching such singular value, so that the positivity structure is
always preserved. Moreover, an unconditional energy stability of the numerical
scheme is derived, without any restriction for the time step size. In the
second order numerical scheme, the BDF temporal stencil is applied, and an
alternate convex-concave decomposition is derived, so that the concave part
corresponds to a quadratic energy. In turn, the combined Leonard-Jones
potential term is treated implicitly, and the concave part the is approximated
by a second order Adams-Bashforth explicit extrapolation, and an artificial
Douglas-Dupont regularization term is added to ensure the energy stability.
The unique solvability and the positivity-preserving property for the second
order scheme could be similarly established. In addition, optimal rate
convergence analysis is provided for both the first and second order accurate
schemes. A few numerical simulation results are also presented, which
demonstrate the robustness of the numerical schemes.

Key words and phrases: droplet liquid film model; unique solvability;
positivity-preserving; energy stability; optimal rate convergence analysis

AMS subject classification: 35K35, 35K55, 49J40, 65M06, 65M12

##  1 Introduction

Certain liquids on a solid, chemo-attractive substrate spontaneously form a
droplet structure connected by a very thin precursor (or wetting) layer. After
the droplets appear, coarsening will occur, whereby smaller droplets will
shrink and larger droplets will grow. The coarsening behavior, especially the
rate of coarsening, of droplets has been of great scientific interest. The
average droplet size increases with the decrease of the number of droplets and
the increase of the characteristic distance. The coarsening process can be
mediated by two mechanisms: collapse or collision. Collapse occurs through
mass exchange between droplets in the precursor/wetting layer, and the
droplets can also drift and collide with each other like particles in the
traditional Ostwald ripening process. Thin layers of viscous liquid are well
described by the lubrication approximation, which capitalizes on the
separation of horizontal and vertical length scales. Under the assumption that
the liquid film does not evaporate, lubrication theory leads to a single
equation for the height function, $\phi=\phi(x,t)>0$, of a time-dependent film
[36], in the form of an $H^{-1}$ gradient flow:

$\partial_{t}\phi=\nabla\cdot\left(\mathcal{M}(\phi)\nabla\delta_{\phi}F\right).$

Here $F$ is the free energy of the film/substrate system and is given by

$F(\phi)=\int_{\Omega}\left(\frac{1}{2}|\nabla\phi|^{2}+\mathcal{U}(\phi)\right)d{\bf
x},$

and the term $\delta_{\phi}F$ is the variational derivative of $F$. In more
detail, the gradient term describes the linearized contribution of the liquid-
air surface tension, while $\mathcal{U}(\phi)$ models the intermolecular
forces between the substrate and the film. See the related research on the
coarsening dynamics [1, 2, 3, 4, 5, 21, 22, 26, 27, 28, 32, 33, 37, 38, 39,
41, 42].

In this paper, we consider a well-known Lennard-Jones-type potential [25],
which is expressed as
$\mathcal{U}(\phi)=\frac{1}{3}\phi^{-8}-\frac{4}{3}\phi^{-2}$. Other
potentials are also physically relevant [39]. The free energy is

$F(\phi)=\int_{\Omega}\left(\frac{1}{3}\phi^{-8}-\frac{4}{3}\phi^{-2}+\frac{\varepsilon^{2}}{2}|\nabla\phi|^{2}\right)d{\bf
x},$ (1.1)

where $\Omega\subset\mathbb{R}^{3}$ is a bounded domain. The $H^{-1}$ gradient
flow associated with the given free energy functional (1.1) with constant
mobility $\mathcal{M}\equiv 1$ (non-constant mobility case could be handled in
a similar way) is

$\partial_{t}\phi=\Delta\mu,\quad\mu:=\delta_{\phi}F=-\frac{8}{3}(\phi^{-9}-\phi^{-3})-\varepsilon^{2}\Delta\phi.$
(1.2)

Mass is conserved, and, due to the gradient structure, the following energy
dissipation law is formally available:

$\frac{d}{dt}F(\phi(t))=-\int_{\Omega}\,|\nabla\mu|^{2}\,d{\bf x}\leq 0.$

For simplicity of presentation, we assume periodic boundary conditions hold
over the rectangular domain $\Omega$. Other types of boundary conditions, such
as homogeneous Neumann, can also be handled.

Some numerical simulation results have been reported for the droplet liquid
film model [21, 36], while the numerical analysis for this model has been very
limited. The highly nonlinear and highly singular nature of the Leonard-Jones
potential makes the equation very challenging, at both the analytic and
numerical levels. In particular, the structure of the potential function
requires that the phase variable has to keep a fixed sign, i.e., being, either
positive or negative, to avoid a singularity. A positivity-preserving
structure is required for the numerical scheme and for physical reality. The
theoretical justification of such a property, in addition to the energy
dissipation property, turns out to be a major difficulty for the numerical
analysis. Moreover, a convergence analysis for any numerical scheme for the
liquid film equation (1.2) has also been a long-standing open problem.

In this paper we construct and analyze finite difference numerical schemes
that preserve all these three important theoretical features, in both the
first and second order temporal accuracy orders. In the first order scheme, we
make use of the convex-concave decomposition of the physical energy, and come
up with a semi-implicit numerical method. The positivity-preserving property
will be theoretically established, which is based on the subtle fact that the
numerical solution is equivalent to a minimization of a convex energy
functional, and the singular nature of the $\phi^{-8}$ term around the value
of $0$ prevents the numerical solution reaching the singular value. As a
result, the positivity structure is always preserved, and the unique
solvability analysis comes from the convex nature of the implicit parts in the
scheme. The energy stability becomes a direct consequence of this convexity
analysis. In the second order numerical scheme, we have to make use of an
alternate convex-concave decomposition of the physical energy, by observing
that $\frac{1}{3}\phi^{-8}-\frac{4}{3}\phi^{-2}+\frac{4}{3}A_{0}\phi^{2}$ is a
convex function, provided $A_{0}$ is sufficiently large. In turn, the concave
part becomes a quadratic energy, which is approximated by a second order
Adams-Bashforth explicit extrapolation. The second order BDF temporal stencil
is applied, and an artificial Douglas-Dupont regularization term is added to
ensure energy stability.

In addition, a detailed convergence analysis of the proposed numerical
schemes, including the first and second order accurate ones, will be derived
in this paper. This analysis gives an optimal rate error estimate in the
$\ell^{\infty}(0,T;H_{h}^{-1})\cap\ell^{2}(0,T;H_{h}^{1})$. A key point in the
analysis lies in the following fact: since the nonlinear potential term
corresponds to a convex energy, the corresponding nonlinear error inner
product is always non-negative. Furthermore, the error estimate associated
with the surface diffusion term indicates an $\ell^{2}(0,T;H_{h}^{1})$
convergence. It will be the first such work for this model.

The rest of the paper is organized as follows. In Section 2, we describe the
finite difference discretization of space and recall some basic facts. Section
3, we propose the first order numerical scheme and state the main theoretical
results, including positivity-preserving, unique solvability, and
unconditional energy stability of the scheme. Further, we provide a detailed
convergence analysis. The second order accurate scheme is proposed in Section
4; the unique solvability, the positivity-preserving property and an optimal
rate convergence analysis are established as well, using similar ideas. Some
numerical simulation results are presented in Section 5. Finally, the
concluding remarks are given in Section 6.

##  2 Finite difference spatial discretization

The model described above is inherently two-dimensional in nature. However,
other models of this type may be three-dimensional. Our analyses to follow,
therefore, will generically cover one through three space dimensions. For
spatial approximation we use a staggered-grid finite difference approach. We
use the notation and results for some discrete functions and operators from
[23, 44, 45], _et cetera_. Let
$\Omega=(0,L_{x})\times(0,L_{y})\times(0,L_{z}),$ where for simplicity, we
assume $L_{x}=L_{y}=L_{z}=:L>0$. Let $N\in\mathbb{N}$ be given, and define the
grid spacing $h:=\frac{L}{N}$, i.e., a uniform mesh assumption is made for
simplicity of notation. We define the following two uniform, infinite grids
with spacing $h>0$: $C:=\\{p_{i}\,|\,i\in\mathbb{Z}\\}$ and
$E:=\\{p_{i+1/2}\,|\,i\in\mathbb{Z}\\}$, where $p_{i}=p(i):=(i-1/2)h.$
Consider the following 3D discrete periodic function spaces:

$\begin{split}\mathcal{C}_{\rm per}:=\left\\{\upsilon:C\times C\times
C\rightarrow\mathbb{R}\,\middle|\,\upsilon_{i,j,k}=\upsilon_{i+\alpha
N,j+\beta N,k+\gamma N},\
\forall\,i,j,k,\alpha,\beta,\gamma\in\mathbb{Z},\right\\},\\\ \mathcal{E}_{\rm
per}^{x}:=\left\\{\upsilon:E\times C\times
C\rightarrow\mathbb{R}\,\middle|\,\upsilon_{i+\frac{1}{2},j,k}=\upsilon_{i+\frac{1}{2}+\alpha
N,j+\beta N,k+\gamma N},\
\forall\,i,j,k,\alpha,\beta,\gamma\in\mathbb{Z}\right\\},\end{split}$

with the identification $\upsilon_{i,j,k}=\upsilon(p_{i},p_{j},p_{k}),$ _et
cetera_. The spaces $\mathcal{E}_{\rm per}^{y}$ and $\mathcal{E}_{\rm
per}^{z}$ are analogously defined. The functions of $\mathcal{C}_{\rm per}$
are called cell-centered functions, and the functions of $\mathcal{E}_{\rm
per}^{x},$ $\mathcal{E}_{\rm per}^{y}$ and $\mathcal{E}_{\rm per}^{z}$ called
face-centered functions. We also introduce the following mean zero space to
facilitate the later analysis: $\mathring{\mathcal{C}}_{\rm
per}:=\\{\upsilon\in\mathcal{C}_{\rm
per}\,|\,0=\overline{\upsilon}:=\frac{h^{3}}{|\Omega|}\sum_{i,j,k=1}^{N}\upsilon_{i,j,k}\\}$.
We define $\vec{\mathcal{E}}_{\rm per}:=\mathcal{E}_{\rm
per}^{x}\times\mathcal{E}_{\rm per}^{y}\times\mathcal{E}_{\rm per}^{z}.$

The center-to-face difference and averaging operators,
$A_{x},\;D_{x}:\mathcal{C}_{\rm per}\rightarrow\mathcal{E}_{\rm per}^{x}$,
$A_{y},\;D_{y}:\mathcal{C}_{\rm per}\rightarrow\mathcal{E}_{\rm per}^{y}$, and
$A_{z},\;D_{z}:\mathcal{C}_{\rm per}\rightarrow\mathcal{E}_{\rm per}^{z}$, are
defined as follows:

$\begin{split}A_{x}\upsilon_{i+\frac{1}{2},j,k}:=\frac{1}{2}(\upsilon_{i+1,j,k}+\upsilon_{i,j,k}),\;\;\;\;D_{x}\upsilon_{i+\frac{1}{2},j,k}:=\frac{1}{h}(\upsilon_{i+1,j,k}-\upsilon_{i,j,k}),\\\
A_{y}\upsilon_{i,j+\frac{1}{2},k}:=\frac{1}{2}(\upsilon_{i,j+1,k}+\upsilon_{i,j,k}),\;\;\;\;D_{y}\upsilon_{i,j+\frac{1}{2},k}:=\frac{1}{h}(\upsilon_{i,j+1,k}-\upsilon_{i,j,k}),\\\
A_{z}\upsilon_{i,j,k+\frac{1}{2}}:=\frac{1}{2}(\upsilon_{i,j,k+1}+\upsilon_{i,j,k}),\;\;\;\;D_{z}\upsilon_{i,j,k+\frac{1}{2}}:=\frac{1}{h}(\upsilon_{i,j,k+1}-\upsilon_{i,j,k}).\end{split}$

Likewise, the face-to-center difference and averaging operators,
$a_{x},\;d_{x}:\mathcal{E}_{\rm per}^{x}\rightarrow\mathcal{C}_{\rm per}$,
$a_{y},\;d_{y}:\mathcal{E}_{\rm per}^{y}\rightarrow\mathcal{C}_{\rm per}$, and
$a_{z},\;d_{z}:\mathcal{E}_{\rm per}^{z}\rightarrow\mathcal{C}_{\rm per}$, are
defined as follows:

$\begin{split}a_{x}\upsilon_{i,j,k}:=\frac{1}{2}(\upsilon_{i+\frac{1}{2},j,k}+\upsilon_{i-\frac{1}{2},j,k}),\quad
d_{x}\upsilon_{i,j,k}:=\frac{1}{h}(\upsilon_{i+\frac{1}{2},j,k}-\upsilon_{i-\frac{1}{2},j,k}),\\\
a_{y}\upsilon_{i,j,k}:=\frac{1}{2}(\upsilon_{i,j+\frac{1}{2},k}+\upsilon_{i,j-\frac{1}{2},k}),\quad
d_{y}\upsilon_{i,j,k}:=\frac{1}{h}(\upsilon_{i,j+\frac{1}{2},k}-\upsilon_{i,j-\frac{1}{2},k}),\\\
a_{z}\upsilon_{i,j,k}:=\frac{1}{2}(\upsilon_{i,j,k+\frac{1}{2}}+\upsilon_{i,j,k-\frac{1}{2}}),\quad
d_{z}\upsilon_{i,j,k}:=\frac{1}{h}(\upsilon_{i,j,k+\frac{1}{2}}-\upsilon_{i,j,k-\frac{1}{2}}).\end{split}$

The discrete gradient $\nabla_{h}:\mathcal{C}_{\rm
per}\rightarrow\vec{\mathcal{E}}_{\rm per}$ and the discrete divergence
$\nabla_{h}\cdot:\vec{\mathcal{E}}_{\rm per}\rightarrow\mathcal{C}_{\rm per}$
are defined via

$\displaystyle\nabla_{h}\upsilon_{i,j,k}$
$\displaystyle:=\left(D_{x}\upsilon_{i+\frac{1}{2},j,k},D_{y}\upsilon_{i,j+\frac{1}{2},k},D_{z}\upsilon_{i,j,k+\frac{1}{2}}\right),$
$\displaystyle\nabla_{h}\cdot\vec{f}_{i,j,k}$
$\displaystyle:=d_{x}f_{i,j,k}^{x}+d_{y}f_{i,j,k}^{y}+d_{z}f_{i,j,k}^{z},$

where $\upsilon\in\mathcal{C}_{\rm per}$ and
$\vec{f}=(f^{x},f^{y},f^{z})\in\vec{\mathcal{E}}_{\rm per}$. The standard 3D
discrete Laplacian, $\Delta_{h}:\mathcal{C}_{\rm
per}\rightarrow\mathcal{C}_{\rm per}$ is given by

$\displaystyle\Delta_{h}\upsilon_{i,j,k}:$
$\displaystyle=d_{x}(D_{x}v)_{i,j,k}+d_{y}(D_{y}v)_{i,j,k}+d_{z}(D_{z}v)_{i,j,k})$
$\displaystyle=\frac{1}{h^{2}}(\upsilon_{i+1,j,k}+\upsilon_{i-1,j,k}+\upsilon_{i,j+1,k}+\upsilon_{i,j-1,k}+\upsilon_{i,j,k+1}+\upsilon_{i,j,k-1}-6\upsilon_{i,j,k}).$

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

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

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

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

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

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

We define the following norms for grid functions: for
$\upsilon\in\mathcal{C}_{\rm per},$ then
$\|\upsilon\|_{2}^{2}:=\langle\upsilon,\upsilon\rangle_{\Omega};$
$\|\upsilon\|_{p}^{p}:=\langle|\upsilon|^{p},1\rangle_{\Omega},$ for $1\leq
p<\infty,$ and $\|\upsilon\|_{\infty}:=\max_{1\leq i,j,k\leq
N}|\upsilon_{i,j,k}|$. We define norms of the gradient as follows: for
$\upsilon\in\mathcal{C}_{\rm per},$

$\begin{split}\|\nabla_{h}\upsilon\|_{2}^{2}:=[\nabla_{h}\upsilon,\nabla_{h}\upsilon]_{\Omega}=[D_{x}\upsilon,D_{x}\upsilon]_{x}+[D_{y}\upsilon,D_{y}\upsilon]_{y}+[D_{z}\upsilon,D_{z}\upsilon]_{z},\end{split}$

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

$\begin{split}\|\nabla_{h}\upsilon\|_{p}:=([|D_{x}\upsilon|^{p},1]_{x}+[|D_{y}\upsilon|^{p},1]_{y}+[|D_{z}\upsilon|^{p},1]_{z})^{\frac{1}{p}}.\end{split}$

Higher order norms can be defined. For example,

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

######  Lemma 2.1 ([45]).

Let $\mathcal{D}$ be an arbitrary periodic, scalar function defined on all of
the face center points. For any $\psi,\nu\in\mathcal{C}_{\rm per}$ and any
$\vec{f}\in\vec{\mathcal{E}}_{\rm per},$ the following summation by parts
formulas are valid:

$\langle\psi,\nabla_{h}\cdot\vec{f}\rangle_{\Omega}=-[\nabla_{h}\psi,\vec{f}]_{\Omega},\quad\langle\psi,\nabla_{h}\cdot(\mathcal{D}\nabla_{h}\nu)\rangle_{\Omega}=-[\nabla_{h}\psi,\mathcal{D}\nabla_{h}\nu]_{\Omega}.$

To facilitate the convergence analysis, we need to introduce a discrete
analogue of the space $H_{\rm per}^{-1},$ as outlined in [43]. For any
$f\in\mathcal{C}_{\rm per}$, with $\overline{f}:=|\Omega|^{-1}\langle
f,1\rangle_{\Omega}=0$, we define
$\psi:=(-\Delta_{h})^{-1}f\in\mathcal{C}_{\rm per}$ as the periodic solution
of

$-\Delta_{h}\psi=f,\quad\overline{\psi}=0.$

In turn, the inner product $\langle\,\cdot,\,\cdot\,\rangle_{-1,h}$ and the
$\|\cdot\|_{-1,h}$ norm are introduced as

$\langle f,g\rangle_{-1,h}:=\langle
f,(-\Delta)^{-1}g\rangle_{\Omega}=\langle(-\Delta)^{-1}f,g\rangle_{\Omega},\quad\|f\|_{-1,h}:=\sqrt{\langle
f,f\rangle_{-1,h}}.$

######  Lemma 2.2 ([11]).

For any $f\in\mathring{\mathcal{C}}_{\rm per}$, we have

$\|f\|_{4}\leq
C\|f\|_{-1,h}^{\frac{1}{8}}\cdot\|\nabla_{h}f\|_{2}^{\frac{7}{8}},$

for some constant $C>0$ that is independent of $h$ and $f$.

To prove the positivity of our methods, we will need the following result.

######  Lemma 2.3 ([7]).

Suppose that $\phi_{1},\phi_{2}\in\mathcal{C}_{\rm per},$ with
$\phi_{1}-\phi_{2}\in\mathring{\mathcal{C}}_{\rm per}$. Assume that
$0<\phi_{1},\phi_{2}<M_{h}$, for any $M_{h}>0$ that may depend on $h$. Then we
have the following estimate:

$\|(-\Delta_{h})^{-1}(\phi_{1}-\phi_{2})\|_{\infty}\leq CM_{h},$ (2.1)

for some constant $C>0$ that only depends on $\Omega$.

##  3 A first order numerical scheme

We follow a convex-concave decomposition methodology in deriving our scheme,
and consider the following semi-implicit, fully discrete schemes: given
$\phi^{n}\in\mathcal{C}_{\rm per},$ find
$\phi^{n+1},\mu^{n+1}\in\mathcal{C}_{\rm per},$ such that

$\frac{\phi^{n+1}-\phi^{n}}{\Delta
t}=\Delta_{h}\mu^{n+1},\quad\mu^{n+1}=-\frac{8}{3}(\phi^{n+1})^{-9}+\frac{8}{3}(\phi^{n})^{-3}-\varepsilon^{2}\Delta_{h}\phi^{n+1}.$
(3.1)

If solutions exist, it is obvious that the numerical scheme is mass
conservative, i.e.,
$\overline{\phi^{n+1}}=\overline{\phi^{n}}=\cdots=\overline{\phi^{0}}:=\beta_{0}$,
since

$\left\langle\phi^{n+1}-\phi^{n},1\right\rangle_{\Omega}=\Delta
t\left\langle\Delta_{h}\mu^{n+1},1\right\rangle_{\Omega}=-\Delta
t\left[\nabla_{h}\mu^{n+1},\nabla_{h}1\right]_{\Omega}=0.$

We will describe in a later section how to construct the initial data for the
discrete scheme, $\phi^{0}$, from the initial data for the PDE system (1.2).

We define the discrete energy $F_{h}(\phi):\mathcal{C}_{\rm per}\to\mathbb{R}$
as

$F_{h}(\phi)=F_{h,c}(\phi)-F_{h,e}(\phi),\quad
F_{h,c}(\phi)=\frac{1}{3}\langle\phi^{-8},1\rangle_{\Omega}+\frac{\varepsilon^{2}}{2}\|\nabla_{h}\phi\|_{2}^{2},\quad
F_{h,e}(\phi)=\frac{4}{3}\langle\phi^{-2},1\rangle_{\Omega}.$ (3.2)

Observe that both $F_{h,c}(\phi)$ and $F_{h,e}(\phi)$ are convex, and both
require $\phi$ to be positive to be well defined. Next, we prove that positive
numerical solutions for (3.1) exist and are unique.

###  3.1 Existence and uniqueness of positive solutions

######  Theorem 3.1.

Let $\phi^{n}\in\mathcal{C}_{\rm per}$, with $\delta_{0}\leq\phi^{n}$, point-
wise, for some $\delta_{0}>0$. Then, there exists a unique solution
$\phi^{n+1}\in\mathcal{C}_{\rm per}$ to the scheme (3.1), with
$\overline{\phi^{n+1}}=\overline{{\phi}^{n}},$ and $\phi^{n+1}>0$, point-wise.

###### Proof.

If a positive numerical solution exists, then it must satisfy the mass
conservation property. Thus, we will only look for such solutions. Set
$\beta_{0}:=\overline{\phi^{n}}$ and $\varphi^{n}:=\phi^{n}-\beta_{0}$. The
numerical solution of (3.1) must be a minimizer of the following discrete
energy functional:

$\mathcal{J}(\varphi):=\frac{1}{2\Delta
t}\|\varphi-\varphi^{n}\|_{-1,h}^{2}+\frac{1}{3}\langle(\varphi+\beta_{0})^{-8},1\rangle_{\Omega}+\frac{\varepsilon^{2}}{2}\|\nabla_{h}\varphi\|_{2}^{2}+\frac{8}{3}\langle(\varphi+\beta_{0}),(\phi^{n})^{-3}\rangle_{\Omega},$
(3.3)

over the admissible set

$\mathring{V}_{h}:=\\{\varphi\in\mathring{\mathcal{C}}_{\rm
per}\,|\,-\beta_{0}<\varphi\leq M_{h}-\beta_{0}\\}\subset\mathbb{R}^{N^{3}},$
(3.4)

where $M_{h}=\frac{\beta_{0}|\Omega|}{h^{3}}$. In fact, such a bound comes
from the fact that, the numerical value of $\varphi$ at any numerical grid
point $(i,j,k)$ has to be bounded by $M_{h}-\beta_{0}$, due to the fact that
$\sum_{i,j,k}^{N}\phi_{i,j,k}=\frac{\beta_{0}|\Omega|}{h^{3}}$, as well as the
positivity assumption for $\phi$: $\phi_{i,j,k}>0$, for any $(i,j,k)$. The
function $\mathcal{J}$ is a strictly convex function over this domain.
Consider the following closed domain: for any $\delta\in(0,1)$,

$\mathring{V}_{\delta,h}:=\\{\varphi\in\mathring{\mathcal{C}}_{\rm
per}\,|\,\delta-\beta_{0}\leq\varphi\leq
M_{h}-\beta_{0}\\}\subset\mathbb{R}^{N^{3}}.$ (3.5)

Since $\mathring{V}_{\delta,h}$ is a bounded, compact, and convex set in the
subspace $\mathring{\mathcal{C}}_{\rm per},$ there exists a (not necessarily
unique) minimizer of $\mathcal{J}$ over $\mathring{V}_{\delta,h}.$ The key
point of the positivity analysis is that such a minimizer can not occur on the
left boundary, when $\delta$ is sufficiently small.

To get a contradiction, assume that the minimizer of $\mathcal{J}$, denoted
$\varphi^{\star}$, occurs at the left boundary of $\mathring{V}_{\delta,h}$.
In particular, assume that

$\varphi_{\vec{\alpha_{0}}}^{\star}=\delta-\beta_{0},$

where $\vec{\alpha_{0}}:=(i_{0},j_{0},k_{0})$. Thus, the grid function
$\varphi^{\star}$ has a global minimum at the grid point $\vec{\alpha_{0}}$.
Suppose that $\varphi^{\star}$ reaches its global maximum at the grid point
$\vec{\alpha_{1}}:=(i_{1},j_{1},k_{1})$. By the the fact that
$\overline{\varphi^{\star}}=0$, it is obvious that
$M_{h}\geq\varphi_{\vec{\alpha_{1}}}^{\star}+\beta_{0}\geq\beta_{0}$.

The function $\mathcal{J}$ is smooth over $\mathring{V}_{\delta,h}$, and, for
any $\psi\in\mathring{\mathcal{C}}_{\rm per}$, the directional derivative is
precisely

$\begin{split}d_{s}\mathcal{J}(\varphi^{\star}+s\psi)|_{s=0}&:=\frac{1}{\Delta
t}\langle(-\Delta_{h})^{-1}(\varphi^{\star}-\varphi^{n}),\psi\rangle_{\Omega}-\frac{8}{3}\langle(\varphi^{\star}+\beta_{0})^{-9},\psi\rangle_{\Omega}\\\
&\quad-\varepsilon^{2}\langle\Delta_{h}\varphi^{\star},\psi\rangle_{\Omega}+\frac{8}{3}\langle(\phi^{n})^{-3},\psi\rangle_{\Omega}.\end{split}$
(3.6)

Let us pick the particular direction $\psi\in\mathring{\mathcal{C}}_{\rm
per}$, satisfying

$\psi_{i,j,k}=\delta_{i,i_{0}}\delta_{j,j_{0}}\delta_{k,k_{0}}-\delta_{i,i_{1}}\delta_{j,j_{1}}\delta_{k,k_{1}},$

where $\delta_{i,j}$ is the Dirac delta function. Then the derivative may be
expressed as

$\begin{split}\frac{1}{h^{3}}d_{s}\mathcal{J}(\varphi^{\star}+s\psi)|_{s=0}:=&\frac{1}{\Delta
t}(-\Delta_{h})^{-1}(\varphi^{*}_{\vec{\alpha_{0}}}+\beta_{0}-\phi^{n}_{\vec{\alpha_{0}}})-\frac{1}{\Delta
t}(-\Delta_{h})^{-1}(\varphi^{\star}_{\vec{\alpha_{1}}}+\beta_{0}-\phi^{n}_{\vec{\alpha_{1}}})\\\
&\quad-\frac{8}{3}(\varphi^{\star}_{\vec{\alpha_{0}}}+\beta_{0})^{-9}+\frac{8}{3}(\varphi^{\star}_{\vec{\alpha_{1}}}+\beta_{0})^{-9}\\\
&\quad+\frac{8}{3}(\phi^{n}_{\vec{\alpha_{0}}})^{-3}-\frac{8}{3}(\phi^{n}_{\vec{\alpha_{1}}})^{-3}-\varepsilon^{2}(\Delta_{h}\varphi^{\star}_{\vec{\alpha_{0}}}-\Delta_{h}\varphi^{\star}_{\vec{\alpha_{1}}}).\end{split}$
(3.7)

Let us define $\phi^{\star}:=\varphi^{\star}+\beta_{0}$. Because
$\phi^{\star}_{\vec{\alpha_{0}}}=\delta$ and
$\phi^{\star}_{\vec{\alpha_{1}}}\geq\beta_{0}$, we have

$-\frac{8}{3}(\phi^{\star}_{\vec{\alpha_{0}}})^{-9}+\frac{8}{3}(\phi^{\star}_{\vec{\alpha_{1}}})^{-9}\leq-\frac{8}{3}(\delta^{-9}-\beta_{0}^{-9}).$
(3.8)

Since $\phi^{\star}$ has a global minimum at the grid point
$\vec{\alpha_{0}}$, with,
$\phi^{\star}_{\vec{\alpha_{0}}}=\delta\leq\phi^{\star}_{i,j,k}$, and a
maximum at the grid point $\vec{\alpha_{1}}$, with
$\phi^{\star}_{\vec{\alpha_{1}}}\geq\phi^{\star}_{i,j,k}$, for any $(i,j,k)$,
the following results are valid:

$\Delta_{h}\varphi^{\star}_{\vec{\alpha_{0}}}\geq
0,\quad\Delta_{h}\varphi^{\star}_{\vec{\alpha_{1}}}\leq 0.$ (3.9)

For the numerical solution $\phi^{n}$ at the previous time step, the a-prior
assumption $\delta_{0}\leq\phi^{n}\leq M_{h}$ indicates that

$-\delta_{0}^{-3}<(\phi^{n}_{\vec{\alpha_{0}}})^{-3}-(\phi^{n}_{\vec{\alpha_{1}}})^{-3}<\delta_{0}^{-3}.$
(3.10)

For the first two terms appearing in (3.7), we apply Lemma 2.3 and obtain

$-2CM_{h}\leq(-\Delta_{h})^{-1}(\phi^{\star}-\phi^{n})_{\vec{\alpha_{0}}}-(-\Delta_{h})^{-1}(\phi^{\star}-\phi^{n})_{\vec{\alpha_{1}}}\leq
2CM_{h}.$ (3.11)

Consequently, a substitution of (3.8) – (3.11) into (3.7) yields

$\frac{1}{h^{3}}d_{s}\mathcal{J}(\varphi^{\star}+s\psi)|_{s=0}\leq-\frac{8}{3}(\delta^{-9}-\beta_{0}^{-9})+\delta_{0}^{-3}+2C\Delta
t^{-1}M_{h}.$ (3.12)

Define $D_{0}:=\delta_{0}^{-3}+2C\Delta t^{-1}M_{h}$, and note that $D_{0}$ is
a constant for fixed $\Delta t$ and $h$, even though it becomes singular as
$\Delta t,h\to 0$. In particular, $D_{0}=O(\Delta
t^{-1}h^{-2}+\delta_{0}^{-3})$. In any case, for any fixed $\Delta t$ and $h$,
we can choose $\delta$ small enough so that

$-\frac{8}{3}(\delta^{-9}-\beta_{0}^{-9})+D_{0}<0.$ (3.13)

This shows that

$\frac{1}{h^{3}}d_{s}\mathcal{J}(\varphi^{\star}+s\psi)|_{s=0}<0,$ (3.14)

for $\delta$ satisfying (3.13). This contradicts the assumption that
$\mathcal{J}$ has a minimum at $\varphi^{\star}$. Therefore, the global
minimum of $\mathcal{J}(\varphi)$ over $\mathring{V}_{\delta,h}$ could only
possibly occur at interior point, if $\delta>0$ is small. We conclude that
there must be a solution $\varphi\in\mathring{\mathcal{C}}_{\rm per}$, so that

$\frac{1}{h^{3}}d_{s}\mathcal{J}(\varphi+s\psi)|_{s=0}=0,$ (3.15)

such that $\phi=\varphi+\beta_{0}$ is positive. Such a minimizer is equivalent
to the numerical solution of (3.1). The existence of a positive numerical
solution is established.

In addition, since $\mathcal{J}(\varphi)$ is a strictly convex function over
$V_{h}$, the uniqueness analysis for this numerical solution is
straightforward. The proof of Theorem 3.1 is complete. ∎

######  Remark 3.1.

This positivity-preserving analysis has been very successfully applied to the
Cahn-Hilliard equation with Flory-Huggins-type logarithmic potentials [7, 13,
14]. It is the first time, to our knowledge, that it has been applied to a
gradient flow with singular Leonard-Jones style potential. Also see related
works [15, 16] for the related maximum bound principles for a class of semi-
linear parabolic equations.

######  Remark 3.2.

In the positivity analysis, it is observed that $\phi^{n+1}\geq\delta$, while
$\delta$ depends on $\Delta t$ and $h$ in a form of $\delta=O(\Delta
t^{\frac{1}{9}}+h^{\frac{2}{9}})$, so that $\delta\to 0$ as $\Delta t,h\to 0$.
in other words, a directly numerical analysis would not be able to ensure a
uniform distance between the numerical solution and the singular limit value
of $0$. On the other hand, in the 1-D and 2-D cases, the phase separation is
expected to be available at the PDE analysis level, i.e., a uniform distance
between the phase variable and the singular limit value $0$, with such a
uniform distance dependent on $\varepsilon$, $\theta_{0}$ and the initial
data. Such a theoretical property has been established for the Cahn-Hilliard
model with logarithmic Flory-Huggins energy potential; see the related
existing works [12, 20, 35], etc. With such a theoretical bound available for
the PDE solution, we are also able to derive a similar bound for the numerical
solution, in combination with the convergence analysis and error estimate
presented in the later section.

Because of the convex-concave decomposition structure of the numerical scheme,
an unconditional energy stability is available.

######  Theorem 3.2.

The scheme (3.1) is unconditional energy stable: for any time step $\Delta
t>0$,

$\begin{split}F_{h}(\phi^{k+1})+\Delta t\|\nabla_{h}\mu^{k+1}\|_{2}^{2}\;\leq
F_{h}(\phi^{k}).\end{split}$

###### Proof.

The following estimate is valid:

$\begin{split}F_{h}(\phi^{k+1})-F_{h}(\phi^{k})&=F_{h,c}(\phi^{k+1})-F_{h,e}(\phi^{k+1})-F_{h,c}(\phi^{k})+F_{h,e}(\phi^{k})\\\
&=F_{h,c}(\phi^{k+1})-F_{h,c}(\phi^{k})-(F_{h,e}(\phi^{k+1})-F_{h,e}(\phi^{k}))\\\
&\leq\langle\delta_{\phi}F_{h,c}(\phi^{k+1}),\phi^{k+1}-\phi^{k}\rangle_{\Omega}-\langle\delta_{\phi}F_{h,e}(\phi^{k}),\phi^{k+1}-\phi^{k}\rangle_{\Omega}\\\
&=\langle\delta_{\phi}F_{h,c}(\phi^{k+1})-\delta_{\phi}F_{h,e}(\phi^{k}),\phi^{k+1}-\phi^{k}\rangle_{\Omega}\\\
&=\langle\mu^{k+1},\phi^{k+1}-\phi^{k}\rangle_{\Omega}\\\ &=\Delta
t\langle\mu^{k+1},\Delta_{h}\mu^{k+1}\rangle_{\Omega}=-\Delta
t\|\nabla_{h}\mu^{k+1}\|_{2}^{2}\leq 0.\end{split}$

∎

###  3.2 The $\ell^{\infty}(0,T;H_{h}^{-1})\cap\ell^{2}(0,T;H_{h}^{1})$
convergence analysis

By $\Phi$ we denote the exact solution for the PDE (1.2). Define
$\Phi_{N}(\,\cdot\,,t):=\mathcal{P}_{N}\Phi(\,\cdot\,,t)$, where
$\mathcal{P}_{N}:L^{2}(\Omega)\to\mathbb{T}_{K}^{d}$ is the (spatial) Fourier
projection into the space of trigonometric polynomials of degree up to and
including $K$, and $N=2K+1$. In particular, in 3D

$\mathcal{P}_{N}f({\bf x})=\sum_{\left\|{\bf k}\right\|_{1}\leq
K}\hat{f}_{{\bf k}}\exp\left(\frac{2\pi\mathfrak{i}}{L}\,{\bf k}\cdot{\bf
x}\right),\quad\hat{f}_{{\bf
k}}=\frac{1}{L^{3}}\int_{\Omega}f(x)\exp\left(-\frac{2\pi\mathfrak{i}}{L}\,{\bf
k}\cdot{\bf x}\right)d{\bf x}.$

The following projection approximation is standard: if $\Phi\in
L^{\infty}(0,T;H_{\rm per}^{\ell}(\Omega))$, for some $\ell\in\mathbb{N}$,

$\|\Phi_{N}-\Phi\|_{L^{\infty}(0,T;H^{k})}\leq
Ch^{\ell-k}\|\Phi\|_{L^{\infty}(0,T;H^{\ell})},\quad\forall\,k:\,0\leq
k\leq\ell.$ (3.16)

By $\Phi_{N}^{m},\Phi^{m}$ we denote $\Phi_{N}(\,\cdot\,,t_{m}),$ and
$\Phi(\,\cdot\,,t_{m})$, respectively, with $t_{m}=m\cdot\Delta t$. Since
$\Phi_{N}\in\mathbb{T}_{K}^{d}$, the mass conservative property is available
at the time-discrete level:

$\overline{\Phi_{N}^{m}}=\frac{1}{L^{3}}\int_{\Omega}\Phi_{N}(\,\cdot\,,t_{m})\,d{\bf
x}=\frac{1}{L^{3}}\int_{\Omega}\Phi_{N}(\,\cdot\,,t_{m-1})\,d{\bf
x}=\overline{\Phi_{N}^{m-1}},\quad\forall\,m\in\mathbb{N}.$ (3.17)

As shown before, we know that the numerical solution (3.1) is mass
conservative at the fully discrete level:

$\overline{\phi^{m}}=\overline{\phi^{m-1}},\quad\forall\,m\in\mathbb{N}.$
(3.18)

Define the grid projection operator, $\mathcal{P}_{h}:C^{0}_{\rm
per}(\Omega)\to\mathcal{C}_{\rm per}$, by
$\mathcal{P}_{h}f_{i,j,k}=f(p_{i},p_{j},p_{k})$, for all $f\in C^{0}_{\rm
per}(\Omega)$. For the initial data, we take

$\phi^{0}=\mathcal{P}_{h}\Phi_{N}(\,\cdot\,,t=0).$ (3.19)

The error grid function is defined as

$\tilde{\phi}_{m}:=\mathcal{P}_{h}\Phi^{m}_{N}-\phi^{m},\quad\forall\,m\in\\{0,1,2,3,\cdots\\}.$
(3.20)

Because the initial data are chosen using composition of operators
$\mathcal{P}_{h}\mathcal{P}_{N}$, it is not hard to show that

$\overline{\tilde{\phi}_{m}}:=\frac{1}{L^{3}}\left\langle\tilde{\phi}_{m},1\right\rangle_{\Omega}=0,\quad\forall\,m\in\\{0,1,2,3,\cdots\\},$

so that the discrete norm $\|\,\cdot\,\|_{-1,h}$ is well defined for the error
grid function.

In the following theorem, the convergence result for the numerical scheme
(3.1) is stated.

######  Theorem 3.3.

Given initial data $\Phi(\,\cdot\,,t=0)\in C^{6}_{\rm per}(\Omega),$ suppose
the solution of the PDE (1.2), $\Phi$, is in the regularity class

$\mathcal{R}:=H^{2}(0,T;C_{\rm per}(\Omega))\cap H^{1}(0,T;C^{2}_{\rm
per}(\Omega))\cap L^{\infty}(0,T;C^{6}_{\rm per}(\Omega)).$

Then, provided $\Delta t$ and $h$ are sufficiently small, for all positive
integers n, such that $t_{n}\leq T,$ we have the following convergence
estimate for the numerical solution

$\|\tilde{\phi}^{n}\|_{-1,h}+\left(\varepsilon^{2}\Delta
t\sum_{m=1}^{n}\|\nabla_{h}\tilde{\phi}^{m}\|_{2}^{2}\right)^{\frac{1}{2}}\leq
C(\Delta t+h^{2}),$ (3.21)

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

###### Proof.

A careful consistency analysis indicates the following truncation error
estimate:

$\frac{\Phi^{n+1}_{N}-\Phi^{n}_{N}}{\Delta
t}=\Delta_{h}\Big{(}-\frac{8}{3}(\Phi^{n+1}_{N})^{-9}+\frac{8}{3}(\Phi^{n}_{N})^{-3}-\varepsilon^{2}\Delta_{h}\Phi^{n+1}_{N}\Big{)}+\tau^{n},$
(3.22)

with $\|\tau^{n}\|_{-1,h}\leq C(\Delta t+h^{2}).$ Observe that in equation
(3.22), and from this point forward, we drop the operator $\mathcal{P}_{h}$,
which should appear in front of $\Phi_{N}$, for simplicity. In turn,
subtracting the numerical scheme (3.1) from (3.22) gives

$\frac{\tilde{\phi}^{n+1}-\tilde{\phi}^{n}}{\Delta
t}=\Delta_{h}\left(-\frac{8}{3}(\Phi^{n+1}_{N})^{-9}+\frac{8}{3}(\phi^{n+1})^{-9}+\frac{8}{3}(\Phi^{n}_{N})^{-3}-\frac{8}{3}(\phi^{n})^{-3}-\varepsilon^{2}\Delta_{h}\tilde{\phi}^{n+1}\right)+\tau^{n}.$
(3.23)

Taking the discrete inner product of (3.23) with
$2(-\Delta_{h})^{-1}\tilde{\phi}^{n+1}$ yields

$\begin{split}&\langle\tilde{\phi}^{n+1}-\tilde{\phi}^{n},2(-\Delta_{h})^{-1}\tilde{\phi}^{n+1}\rangle_{\Omega}+2\varepsilon^{2}\Delta
t\langle\Delta_{h}^{2}\tilde{\phi}^{n+1},(-\Delta_{h})^{-1}\tilde{\phi}^{n+1}\rangle_{\Omega}\\\
&\quad=\frac{16}{3}\Delta
t\langle(\Phi^{n+1}_{N})^{-9}-(\Phi^{n}_{N})^{-3}-(\phi^{n+1})^{-9}+(\phi^{n})^{-3},\tilde{\phi}^{n+1}\rangle_{\Omega}+2\Delta
t\langle\tau^{n},\tilde{\phi}^{n+1}\rangle_{-1,h}.\end{split}$ (3.24)

For the time difference, the following identity is valid:

$2\langle\tilde{\phi}^{n+1}-\tilde{\phi}^{n},(-\Delta_{h})^{-1}\tilde{\phi}^{n+1}\rangle_{\Omega}=\|\tilde{\phi}^{n+1}\|_{-1,h}^{2}-\|\tilde{\phi}^{n}\|_{-1,h}^{2}+\|\tilde{\phi}^{n+1}-\tilde{\phi}^{n}\|_{-1,h}^{2}.$
(3.25)

The estimate for the surface diffusion term is straightforward:

$\langle\Delta_{h}^{2}\tilde{\phi}^{n+1},(-\Delta_{h})^{-1}\tilde{\phi}^{n+1}\rangle_{\Omega}=\|\nabla_{h}\tilde{\phi}^{n+1}\|_{2}^{2}.$
(3.26)

For the nonlinear inner product, the convex nature of $\frac{1}{3}\phi^{-8}$
implies the following result:

$\langle(\Phi^{n+1}_{N})^{-9}-(\phi^{n+1})^{-9},\tilde{\phi}^{n+1}\rangle_{\Omega}\leq
0.$ (3.27)

For the concave nonlinear inner product, we begin with an application of mean
value theorem:

$(\Phi^{n}_{N})^{-3}-(\phi^{n})^{-3}=-3(\xi^{n})^{-4}\tilde{\phi}^{n},$ (3.28)

where $\xi^{n}\in\mathcal{C}_{\rm per}$ is a grid function that is between
$\Phi^{n}_{N}$ and $\phi^{n}$ in a point-wise sense. Moreover, we see that

$\displaystyle\|(\xi^{n})^{-4}\|_{2}\leq\max\left(\|(\Phi^{n}_{N})^{-1}\|_{8}^{4},\|(\phi^{n})^{-1}\|_{8}^{4}\right)\leq
C^{*},$ (3.29)

which is based on the smoothness of the projection solution $\Phi_{N}$, as
well as the $\|\cdot\|_{8}$ bound of the numerical solution $(\phi^{n})^{-1}$,
coming from the energy stability. Then we arrive at

$\displaystyle-\frac{16}{3}\langle(\Phi^{n}_{N})^{-3}-(\phi^{n})^{-3},\tilde{\phi}^{n+1}\rangle_{\Omega}$
$\displaystyle\leq
16\|(\xi^{n})^{-4}\|_{2}\cdot\|\tilde{\phi}^{n}\|_{4}\cdot\|\tilde{\phi}^{n+1}\|_{4}$
$\displaystyle\leq
16C^{*}\|\tilde{\phi}^{n}\|_{4}\cdot\|\tilde{\phi}^{n+1}\|_{4}\leq
8C^{*}\left(\|\tilde{\phi}^{n}\|_{4}^{2}+\|\tilde{\phi}^{n+1}\|_{4}^{2}\right).$
(3.30)

Meanwhile, for the $\|\,\cdot\,\|_{4}$ norm, an application of Lemma 2.2
indicates that

$\|\tilde{\phi}^{k}\|_{4}^{2}\leq
C\|\tilde{\phi}^{k}\|_{-1,h}^{\frac{1}{4}}\cdot\|\nabla_{h}\tilde{\phi}^{k}\|_{2}^{\frac{7}{4}}\leq\tilde{C}_{2}\varepsilon^{-2}\|\tilde{\phi}^{k}\|_{-1,h}^{2}+\frac{\varepsilon^{2}}{16C^{*}}\|\nabla_{h}\tilde{\phi}^{k}\|_{2}^{2},\quad
k=n,n+1,$ (3.31)

where Young’s inequality has been applied. Then we get

$\displaystyle-\frac{16}{3}\langle(\Phi^{n}_{N})^{-3}-(\phi^{n})^{-3},\tilde{\phi}^{n+1}\rangle_{\Omega}$
$\displaystyle\leq
8C^{*}\tilde{C}_{2}\varepsilon^{-2}\left(\|\tilde{\phi}^{n+1}\|_{-1,h}^{2}+\|\tilde{\phi}^{n}\|_{-1,h}^{2}\right)$
$\displaystyle\quad+\frac{\varepsilon^{2}}{2}\left(\|\nabla_{h}\tilde{\phi}^{n+1}\|_{2}^{2}+\|\nabla_{h}\tilde{\phi}^{n}\|_{2}^{2}\right).$
(3.32)

The term associated with the truncation error can be controlled by the
standard Cauchy inequality:

$\langle\tau^{n},\tilde{\phi}^{n+1}\rangle_{-1,h}\leq\|\tau^{n}\|_{-1,h}\cdot\|\tilde{\phi}^{n+1}\|_{-1,h}\leq\frac{1}{2}\left(\|\tau^{n}\|_{-1,h}^{2}+\|\tilde{\phi}^{n+1}\|_{-1,h}^{2}\right).$
(3.33)

Subsequently, a substitution of (3.25)-(3.33) into (3.24) yields

$\begin{split}&\|\tilde{\phi}^{n+1}\|_{-1,h}^{2}-\|\tilde{\phi}^{n}\|_{-1,h}^{2}+\frac{3}{2}\Delta
t\varepsilon^{2}\|\nabla_{h}\tilde{\phi}^{n+1}\|_{2}^{2}\\\ &\quad\leq\Delta
t\|\tau^{n}\|_{-1,h}^{2}+(8C^{*}\tilde{C}_{2}\varepsilon^{-2}+1)\Delta
t\left(\|\tilde{\phi}^{n+1}\|_{-1,h}^{2}+\|\tilde{\phi}^{n}\|_{-1,h}^{2}\right)+\frac{\varepsilon^{2}}{2}\Delta
t\|\nabla_{h}\tilde{\phi}^{n}\|_{2}^{2}.\end{split}$ (3.34)

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

$\begin{split}\|\tilde{\phi}^{n+1}\|_{-1,h}+\left(\varepsilon^{2}\Delta
t\sum_{k=0}^{n+1}\|\nabla_{h}\tilde{\phi}^{m}\|_{2}^{2}\right)^{\frac{1}{2}}\leq
C(\Delta t+h^{2}),\end{split}$ (3.35)

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

######  Remark 3.3.

It is observed that, a constant in the form of
$8C^{*}\tilde{C}_{2}\varepsilon^{-2}$ has appeared in the coefficient of
$\|\tilde{\phi}^{n+1}\|_{-1,h}^{2}$ on the right hand side of (3.34). As a
result, a constraint for the time step size, namely in the form of
$C^{*}\tilde{C}_{2}\varepsilon^{-2}\Delta t\leq\frac{1}{32}$, is needed to
obtain an inequality in the following format

$\begin{split}&\|\tilde{\phi}^{n+1}\|_{-1,h}^{2}+\varepsilon^{2}\Delta
t\sum_{k=1}^{n+1}\|\nabla_{h}\tilde{\phi}^{k}\|_{2}^{2}\leq\Delta
t\sum_{k=0}^{n}\|\tau^{n}\|_{-1,h}^{2}+(CC^{*}\tilde{C}_{2}\varepsilon^{-2}+1)\Delta
t\sum_{k=0}^{n}\|\tilde{\phi}^{k}\|_{-1,h}^{2},\end{split}$

so that a discrete Gronwall inequality could be effectively applied. In fact,
such a constraint could be viewed as the condition, “provided $\Delta t$ and
$h$ are sufficiently small”, stated in Theorem 3.3.

##  4 A second order numerical scheme

In the first order scheme (3.1), the nonlinear concave term is treated
explicitly. This approach may lead to a difficulty to derive a second order
accurate, energy stable scheme. Instead, we make use of an alternate convexity
analysis. The following preliminary estimate is needed.

######  Lemma 4.1.

For $x>0$, the function
$f(x)=\frac{1}{3}x^{-8}-\frac{4}{3}x^{-2}+\frac{4}{3}A_{0}x^{2}$ is convex,
provided $A_{0}\geq A_{0}^{\star}:=\frac{9}{5}(\frac{2}{15})^{\frac{2}{3}}$.

###### Proof.

A direct calculation gives
$f^{\prime\prime}(x)=\frac{8}{3}(9x^{-10}-3x^{-4}+A_{0})$. Meanwhile, a
careful application of Young inequality reveals that

$y\leq 3y^{\frac{5}{2}}+\frac{1}{3}A_{0}^{\star},\quad\forall\,y>0.$

By setting $y=x^{-4},$ we conclude that $f^{{}^{\prime\prime}}(x)\geq 0$, for
any $x>0$. This, in turn, implies that $f(x)$ is convex for $x>0$. ∎

As a result, we obtain the following alternate decomposition of $F(\phi)$:

$F(\phi)=F_{c}(\phi)-F_{e}(\phi),$

where

$F_{c}(\phi)=\int_{\Omega}\left(\frac{1}{3}\phi^{-8}+\frac{4}{3}A_{0}\phi^{2}-\frac{4}{3}\phi^{-2}+\frac{\varepsilon^{2}}{2}|\nabla\phi|^{2}\right)d{\bf
x},\quad F_{e}(\phi)=\int_{\Omega}\,\frac{4}{3}A_{0}\phi^{2}d{\bf x}.$

We propose the following second order scheme for the wetting and droplet
coarsening equation: given $\phi^{n-1},\phi^{n}\in\mathcal{C}_{\rm per},$ find
$\phi^{n+1},\mu^{n+1}\in\mathcal{C}_{\rm per},$ such that

$\displaystyle\frac{\frac{3}{2}\phi^{n+1}-2\phi^{n}+\frac{1}{2}\phi^{n-1}}{\Delta
t}$ $\displaystyle=\Delta_{h}\mu^{n+1},$ (4.1) $\displaystyle\mu^{n+1}$
$\displaystyle=-\frac{8}{3}((\phi^{n+1})^{-9}-(\phi^{n+1})^{-3})+\frac{8}{3}A_{0}(\phi^{n+1}-\check{\phi}^{n+1})$
$\displaystyle\quad-A\Delta
t\Delta_{h}(\phi^{n+1}-\phi^{n})-\varepsilon^{2}\Delta_{h}\phi^{n+1},$ (4.2)

where $\check{\phi}^{n+1}:=2\phi^{n}-\phi^{n-1}$ and $A\geq 0$ is a
stabilization parameter to be determined. The unique solvability and
positivity-preserving properties for the second order numerical scheme are
established in the following theorem. The initial data,
$\phi^{0},\phi^{-1}\in\mathcal{C}_{\rm per}$, will be specified later.

######  Theorem 4.1.

Given $\phi^{k}\in\mathcal{C}_{\rm per}$, with $\|\phi^{k}\|_{\infty}\leq
M_{h}$, $k=n,n-1$, for some $M_{h}>0$, and
$\overline{\phi^{n}}=\overline{\phi^{n-1}}=\beta_{0}$, there exists a unique
solution $\phi^{n+1}\in\mathcal{C}_{\rm per}$ to (4.1), with
$\phi^{n+1}-\phi^{n}\in\mathring{\mathcal{C}}_{\rm per}$, and $\phi^{n+1}>0$
at a point-wise level.

###### Proof.

We follow the notation in the proof of Theorem 3.1. Solving (4.1) is
equivalent to minimizing the discrete energy functional

$\begin{split}\mathcal{J}(\varphi)&=\frac{1}{3\Delta
t}\left\|\frac{3}{2}(\varphi+\beta_{0})-2\phi^{n}+\frac{1}{2}\phi^{n-1}\right\|_{-1,h}^{2}+\frac{1}{3}\langle(\varphi+\beta_{0})^{-8}-4(\phi+\beta_{0})^{-2},1\rangle_{\Omega}\\\
&\
+4A_{0}\langle(\varphi+\beta_{0})^{2},1\rangle_{\Omega}+\frac{\varepsilon^{2}+A\Delta
t}{2}\|\nabla_{h}\varphi\|_{2}^{2}+\langle(\varphi+\beta_{0}),A\Delta
t\Delta_{h}\phi^{n}-\frac{8}{3}A_{0}\check{\phi}^{n+1}\rangle_{\Omega}.\end{split}$
(4.3)

over the admissible set $\mathring{V}_{h}$, defined in (3.4).

To obtain the existence of a minimizer for $\mathcal{J}$ over
$\mathring{V}_{h}$, we consider the closed domain $\mathring{V}_{\delta,h}$,
defined in (3.5). Since $\mathring{V}_{\delta,h}$ is a bounded, compact, and
convex set in the subspace $\mathring{\mathcal{C}}_{\rm per}$, there exists a
minimizer of $\mathcal{J}$ over $\mathring{V}_{\delta,h}$. We will show that
such a minimizer could not occur on the left boundary, if $\delta$ is small
enough.

To get a contradiction, assume that the minimizer of $\mathcal{J}$ occurs at
the left boundary of $\mathring{V}_{\delta,h}$. Denote the minimizer by
$\varphi^{\star}$, and suppose that
$\varphi_{\vec{\alpha_{0}}}^{\star}=\delta-\beta_{0}$, where
$\vec{\alpha_{0}}=(i_{0},j_{0},k_{0})$. Then the grid function
$\varphi^{\star}$ has a global minimum at $\vec{\alpha_{0}}$. Suppose that
$\varphi^{\star}$ reaches its maximum value at the grid point
$\vec{\alpha_{1}}=(i_{1},j_{1},k_{1})$. Since $\overline{\varphi^{\star}}=0$,
it is obvious that
$M_{h}\geq\varphi_{\vec{\alpha_{1}}}^{\star}+\beta_{0}\geq\beta_{0}$.

Since $\mathcal{J}$ is smooth over $\mathring{V}_{\delta,h},$ for any
$\psi\in\mathring{\mathcal{C}}_{\rm per}$, the directional derivative is

$\begin{split}d_{s}\mathcal{J}(\varphi^{\star}+s\psi)|_{s=0}&:=\frac{1}{\Delta
t}\left\langle(-\Delta_{h})^{-1}\left(\frac{3}{2}(\varphi^{\star}+\beta_{0})-2\phi^{n}+\frac{1}{2}\phi^{n-1}\right),\psi\right\rangle_{\Omega}\\\
&\quad+\frac{8}{3}\langle(-\varphi^{\star}+\beta_{0})^{-9}+(\varphi^{\star}+\beta_{0})^{-3}+A_{0}(\varphi^{\star}+\beta_{0}),\psi\rangle_{\Omega}\\\
&\quad-(\varepsilon^{2}+A\Delta
t)\langle\Delta_{h}\varphi^{\star},\psi\rangle_{\Omega}+\left\langle A\Delta
t\Delta_{h}\phi^{n}-\frac{8}{3}A_{0}\check{\phi}^{n+1},\psi\right\rangle_{\Omega}.\end{split}$

As before, let us select the direction
$\psi_{i,j,k}=\delta_{i,i_{0}}\delta_{j,j_{0}}\delta_{k,k_{0}}-\delta_{i,i_{1}}\delta_{j,j_{1}}\delta_{k,k_{1}}$.
Then, the derivative may be expressed as

$\begin{split}\frac{1}{h^{3}}d_{s}\mathcal{J}(\varphi^{\star}+s\psi)|_{s=0}&:=\frac{1}{\Delta
t}(-\Delta_{h})^{-1}\left(\frac{3}{2}(\varphi^{\star}_{\vec{\alpha_{0}}}+\beta_{0})-2\phi^{n}_{\vec{\alpha_{0}}}+\frac{1}{2}\phi^{n-1}_{\vec{\alpha_{0}}}\right)\\\
&\quad-\frac{1}{\Delta
t}(-\Delta_{h})^{-1}\left(\frac{3}{2}(\varphi^{\star}_{\vec{\alpha_{1}}}+\beta_{0})-2\phi^{n}_{\vec{\alpha_{1}}}+\frac{1}{2}\phi^{n-1}_{\vec{\alpha_{1}}}\right)\\\
&\quad+\frac{8}{3}\left(-(\varphi^{\star}_{\vec{\alpha_{0}}}+\beta_{0})^{-9}+(\varphi^{\star}_{\vec{\alpha_{0}}}+\beta_{0})^{-3}+A_{0}(\varphi^{\star}_{\vec{\alpha_{0}}}+\beta_{0})\right)\\\
&\quad-\frac{8}{3}\left(-(\varphi^{\star}_{\vec{\alpha_{1}}}+\beta_{0})^{-9}+(\varphi^{\star}_{\vec{\alpha_{1}}}-\beta_{0})^{-3}+A_{0}(\varphi^{\star}_{\vec{\alpha_{1}}}+\beta_{0})\right)\\\
&\quad-(\varepsilon^{2}+A\Delta
t)(\Delta_{h}\varphi^{\star}_{\vec{\alpha_{0}}}-\Delta_{h}\varphi^{\star}_{\vec{\alpha_{1}}})+A\Delta
t(\Delta_{h}\phi^{n}_{\vec{\alpha_{0}}}-\Delta_{h}\phi^{n}_{\vec{\alpha_{1}}})\\\
&\quad-\frac{8}{3}A_{0}(\check{\phi}^{n+1}_{\vec{\alpha_{0}}}-\check{\phi}^{n+1}_{\vec{\alpha_{1}}}).\end{split}$
(4.4)

The following estimates are derived

$\begin{split}-5\tilde{C}M_{h}\Delta t^{-1}&\leq\frac{1}{\Delta
t}(-\Delta_{h})^{-1}\left(\frac{3}{2}\phi^{\star}_{\vec{\alpha_{0}}}-2\phi^{n}_{\vec{\alpha_{0}}}+\frac{1}{2}\phi^{n-1}_{\vec{\alpha_{0}}}\right)\\\
&\quad-\frac{1}{\Delta
t}(-\Delta_{h})^{-1}\left(\frac{3}{2}\phi^{\star}_{\vec{\alpha_{1}}}-2\phi^{n}_{\vec{\alpha_{1}}}+\frac{1}{2}\phi^{n-1}_{\vec{\alpha_{1}}}\right)\leq
5\tilde{C}M_{h}\Delta t^{-1},\end{split}$ (4.5)
$\Delta_{h}\varphi^{\star}_{\vec{\alpha_{0}}}\geq
0,\quad\Delta_{h}\varphi^{\star}_{\vec{\alpha_{1}}}\leq 0,$ (4.6)
$\Delta_{h}\phi^{n}_{\vec{\alpha_{0}}}\leq\frac{6M_{h}}{h^{2}},\quad\Delta_{h}\phi^{n}_{\vec{\alpha_{1}}}<\frac{6M_{h}}{h^{2}},$
(4.7)

and

$-3M_{h}\leq\check{\phi}^{n+1}_{\vec{\alpha_{0}}}-\check{\phi}^{n+1}_{\vec{\alpha_{1}}}\leq
3M_{h},$ (4.8)

where we have made use of the fact that
$\|\phi^{k}\|_{\infty}<M_{h},\;k=n,n-1$, and where we have applied Lemma 2.3.
Similar to the argument in (3.12), we obtain the following inequalities, based
on the fact that $\varphi_{\vec{\alpha_{0}}}^{\star}+\beta_{0}=\delta$,
$\varphi_{\vec{\alpha_{1}}}^{\star}+\beta_{0}\geq\beta_{0}$:

$\displaystyle\frac{8}{3}(-(\varphi^{\star}_{\vec{\alpha_{0}}}+\beta_{0})^{-9}+(\varphi^{\star}_{\vec{\alpha_{0}}}+\beta_{0})^{-3}+A_{0}(\varphi^{\star}_{\vec{\alpha_{0}}}+\beta_{0})=\frac{8}{3}(-\delta^{-9}+\delta^{-3}+A_{0}\delta),$
(4.9)
$\displaystyle\frac{8}{3}(-(\varphi^{\star}_{\vec{\alpha_{1}}}+\beta_{0})^{-9}+(\varphi^{\star}_{\vec{\alpha_{1}}}+\beta_{0})^{-3}+A_{0}(\varphi^{\star}_{\vec{\alpha_{1}}}+\beta_{0}))\geq\frac{8}{3}(-\beta_{0}^{-9}+(\beta_{0})^{-3}+A_{0}\beta_{0}).$
(4.10)

Subsequently, a substitution of (4.5) to (4.10) into (4.4) yields the
following bound:

$\begin{split}\frac{1}{h^{3}}d_{s}\mathcal{J}(\varphi^{\star}+s\psi)|_{s=0}&\leq\frac{8}{3}(-\delta^{-9}+\delta^{-3}+\beta_{0}^{-9}+A_{0}(\delta-\beta_{0}))\\\
&\quad+8A_{0}M_{h}+12M_{h}A\Delta th^{-2}+10\tilde{C}M_{h}\Delta
t^{-1}.\end{split}$ (4.11)

We observe that $(\Delta t)^{-1}$ and $\Delta th^{-2}$ are constants, since
$\Delta t$ and $h$ are assumed fixed, though the terms become singular as
$\Delta t,h\rightarrow 0$. The rest of the analysis follows the same arguments
as in the proof of Theorem 3.1; the details are left to interested readers.
This finishes the proof of Theorem 4.1. ∎

In the following Theorem, we could prove that a modified energy stability is
available for the second order $BDF$ scheme (4.1), provided that
$A\geq\frac{4}{9}A_{0}^{2}$.

######  Theorem 4.2.

With the same assumptions as in Theorem 4.1, we have the following stability
for the proposed numerical scheme (4.1): provided $A\geq\frac{4}{9}A_{0}^{2}$,

$\tilde{F}_{h}(\phi^{n+1},\phi^{n})\leq\tilde{F}_{h}(\phi^{n},\phi^{n-1}),$
(4.12)

where

$\tilde{F}_{h}(\phi^{n+1},\phi^{n}):=F_{h}(\phi^{n+1})+\frac{1}{4\Delta
t}\|\phi^{n+1}-\phi^{n}\|_{-1,h}^{2}+\frac{4}{3}A_{0}\|\phi^{n+1}-\phi^{n}\|_{2}^{2},$
(4.13)

and

$F_{h}(\phi):=\frac{1}{3}\langle\phi^{-8},1\rangle_{\Omega}+\frac{\varepsilon^{2}}{2}\|\nabla_{h}\phi\|_{2}^{2}-\frac{4}{3}\langle\phi^{-2},1\rangle_{\Omega},$

for any $\Delta t$, $h>0$.

###### Proof.

By taking the inner products of the terms of (4.1) with
$(-\Delta_{h})^{-1}(\phi^{n+1}-\phi^{n})$, we can derive the following
inequalities: for the time difference,

$\begin{split}&\frac{1}{\Delta
t}\left\langle\frac{3}{2}\phi^{n+1}-2\phi^{n}+\frac{1}{2}\phi^{n-1},(-\Delta_{h})^{-1}\left(\phi^{n+1}-\phi^{n}\right)\right\rangle_{\Omega}\\\
&\quad=\frac{3}{2\Delta t}\|\phi^{n+1}-\phi^{n}\|_{-1,h}^{2}-\frac{1}{2\Delta
t}\langle\phi^{n+1}-\phi^{n},\phi^{n}-\phi^{n-1}\rangle_{-1,h}\\\
&\quad\geq\frac{1}{\Delta
t}\left(\frac{5}{4}\|\phi^{n+1}-\phi^{n}\|_{-1,h}^{2}-\frac{1}{4}\|\phi^{n}-\phi^{n-1}\|_{-1,h}^{2}\right);\end{split}$
(4.14)

for the singular terms,

$\begin{split}&\left\langle-\Delta_{h}\left(-\frac{8}{3}((\phi^{n+1})^{-9}-(\phi^{n+1})^{-3})+\frac{8}{3}A_{0}\phi^{n+1}\right),(-\Delta_{h})^{-1}\left(\phi^{n+1}-\phi^{n}\right)\right\rangle_{\Omega}\\\
&\quad=\left\langle-\frac{8}{3}\left((\phi^{n+1})^{-9}+(\phi^{n+1})^{-3}\right)+\frac{8}{3}A_{0}\phi^{n+1},\phi^{n+1}-\phi^{n}\right\rangle_{\Omega}\\\
&\quad\geq\frac{1}{3}\left(\left\langle(\phi^{n+1})^{-8}-4(\phi^{n+1})^{-2}+4A_{0}(\phi^{n+1})^{2},1\rangle_{\Omega}-\langle(\phi^{n})^{-8}-4(\phi^{n})^{-2}+4A_{0}(\phi^{n})^{2},1\right\rangle_{\Omega}\right);\end{split}$
(4.15)

for the highest-order diffusion term,

$\begin{split}\left\langle\Delta_{h}^{2}\phi^{n+1},(-\Delta_{h})^{-1}(\phi^{n+1}-\phi^{n})\right\rangle_{\Omega}&=\left\langle\nabla_{h}\phi^{n+1},\nabla_{h}(\phi^{n+1}-\phi^{n})\right\rangle_{\Omega}\\\
&=\frac{1}{2}\left(\|\nabla_{h}\phi^{n+1}\|_{2}^{2}-\|\nabla_{h}\phi^{n}\|_{2}^{2}+\|\nabla_{h}(\phi^{n+1}-\phi^{n})\|_{2}^{2}\right);\end{split}$
(4.16)

for stabilization terms,

$A\Delta
t\left\langle\Delta_{h}^{2}\left(\phi^{n+1}-\phi^{n}\right),(-\Delta_{h})^{-1}(\phi^{n+1}-\phi^{n})\right\rangle_{\Omega}=A\Delta
t\|\nabla_{h}(\phi^{n+1}-\phi^{n})\|_{2}^{2};$ (4.17)

and for the extrapolation terms,

$\begin{split}&\left\langle\Delta_{h}(2\phi^{n}-\phi^{n-1}),(-\Delta_{h})^{-1}(\phi^{n+1}-\phi^{n})\right\rangle_{\Omega}\\\
&\quad=-\left\langle
2\phi^{n}-\phi^{n-1},\phi^{n+1}-\phi^{n}\right\rangle_{\Omega}\geq-\frac{1}{2}(\|\phi^{n+1}\|_{2}^{2}-\|\phi^{n}\|_{2}^{2})-\frac{1}{2}\|\phi^{n}-\phi^{n-1}\|_{2}^{2}.\end{split}$
(4.18)

Estimate (4.15) is based on the convexity of

$\frac{1}{3}\phi^{-8}-\frac{4}{3}\phi^{-2}+\frac{4}{3}A_{0}\phi^{2},$

as guaranteed by Lemma 4.1. Meanwhile, an application of Cauchy’s inequality
yields

$\frac{1}{\Delta t}\|\phi^{n+1}-\phi^{n}\|_{-1,h}^{2}+A\Delta
t\|\nabla_{h}(\phi^{n+1}-\phi^{n})\|_{2}^{2}\geq
2A^{\frac{1}{2}}\|\phi^{n+1}-\phi^{n}\|_{2}^{2}.$ (4.19)

A combination of (4.14)-(4.19) yields

$\begin{split}&F_{h}(\phi^{n+1})-F_{h}(\phi^{n})+\frac{1}{4\Delta
t}\left(\|\phi^{n+1}-\phi^{n}\|_{-1,h}^{2}-\|\phi^{n}-\phi^{n-1}\|_{-1,h}^{2}\right)\\\
&+\frac{4}{3}A_{0}\left(\|\phi^{n+1}-\phi^{n}\|_{2}^{2}-\|\phi^{n}-\phi^{n-1}\|_{2}^{2}\right)\\\
&\hskip
72.26999pt\leq\left(-2A^{\frac{1}{2}}+\frac{4}{3}A_{0}\right)\|\phi^{n+1}-\phi^{n}\|_{2}^{2}\\\
&\hskip 72.26999pt\leq 0,\end{split}$ (4.20)

provided that $A\geq\frac{4}{9}A_{0}^{2}$. Therefore, we get the energy
estimate (4.12). This completes the proof of Theorem 4.2. ∎

######  Remark 4.1.

There have been a few recent works describing BDF-type schemes for certain
gradient flow models – such as the Cahn-Hilliard equation [9, 46], the
epitaxial thin film equation [19, 24, 31, 34], and the square phase field
crystal model [10] – in which an energy stability was theoretically
established. In the present case, a Douglas-Dupont-type regularization has to
be included in the BDF2 scheme, and a careful analysis reveals an energy
stability with respect to a modified energy. Meanwhile, the present analysis
for a gradient flow with a strongly singular energy potential is non-trivial;
but see the related analyses in [7, 14].

######  Remark 4.2.

In the proposed second order numerical scheme (4.1)-(4.2), an additional
quadratic term $\frac{4}{3}A_{0}\phi^{2}$ makes the nonlinear singular energy
functional convex, so that the energy stability analysis could be derived.
Meanwhile, there are some alternate ways to get a convex-concave decomposition
of the physical energy, such as taking a second order approximation to the
concave term $\phi^{-3}$ at $t^{n+1}$ using an explicit extrapolation formula:
$2(\phi^{n})^{-3}-(\phi^{n-1})^{-3}$. In fact, such an approach has been
widely applied to the standard Cahn-Hilliard equation, in which the concave
term is linear; see the related works [9, 46]. However, for the droplet liquid
film model (1.2), this approach will lead to a theoretical difficulty to
justify the energy stability of the corresponding numerical scheme, due to the
singular and nonlinear nature of the concave term $\phi^{-3}$. Instead, we
combine the singular convex term and singular concave term, as well as an
additional quadratic term, so that their combined energy becomes convex, and
this approach leads to perfect theoretical properties. Also see a few related
works [29, 30] for the analysis of the stabilization numerical methods applied
to the standard Cahn-Hilliard equation.

######  Remark 4.3.

The proposed scheme (4.1) is a two-step method. In addition to the initial
data $\phi^{0}$, as given by (3.19), we need a “ghost” point data,
$\phi^{-1}$, a numerical approximation to the phase variable at
$t^{-1}=-\Delta t$. Due to the long-stencil nature of the BDF2 temporal
discretization, we take $\phi^{-1}=\phi^{0}-\Delta t\Delta_{h}\mu_{h}^{0}$, in
which $\mu_{h}$ stands for the finite difference approximation to the chemical
potential; also see a similar treatment for the Cahn-Hilliard model in [23].
Such a “ghost” point initialization gives a local consistency estimate
$\phi^{-1}-\Phi^{-1}=O(\Delta t^{2}+h^{2})$, while its substitution into (4.1)
implies a first order truncation error, at $n=0$. On the other hand, because
of the temporal derivative involved, such a fact does not spoil the truncation
error for $n\geq 1$.

This approach does not affect the positivity-preserving property of the
numerical solution for $\phi^{n+1}$, $n\geq 0$, since $\phi^{0}$ and
$\phi^{-1}$ always stay bounded within one time step. Such an initialization
method does not affect the modified energy stability estimate (4.12), either.
Meanwhile, to derive a uniform bound for the original energy functional of the
numerical solution, we have to add two $O(\Delta t^{2})$ correction terms:

$\displaystyle F_{h}(\phi^{n+1})\leq$
$\displaystyle\tilde{F}_{h}(\phi^{n+1},\phi^{n})\leq\tilde{F}_{h}(\phi^{n},\phi^{n-1})...$
$\displaystyle\leq$
$\displaystyle\tilde{F}_{h}(\phi^{0},\phi^{-1})=F_{h}(\phi^{0})+\frac{1}{4\Delta
t}\|\phi^{0}-\phi^{-1}\|_{-1,h}^{2}+\frac{4}{3}A_{0}\|\phi^{0}-\phi^{-1}\|_{2}^{2}.$

The convergence result for the second order scheme is stated in the following
theorem.

######  Theorem 4.3.

Given initial data $\Phi(\,\cdot\,,t=0)\in C^{6}_{\rm per}(\Omega)$, suppose
the exact solution for (1.2) is of regularity class

$\mathcal{R}_{2}:=H^{3}(0,T;C_{\rm per}(\Omega))\cap H^{2}(0,T;C^{2}_{\rm
per}(\Omega))\cap L^{\infty}(0,T;C^{6}_{\rm per}(\Omega)).$

Then, provided $\Delta t$ and $h$ are sufficiently small, for all positive
integers $n$, such that $t_{n}\leq T,$ we have the following convergence
estimate for the numerical solution (4.1)

$\|\tilde{\phi}^{n}\|_{-1,h}+\Bigl{(}\varepsilon^{2}\Delta
t\sum_{m=1}^{n}\|\nabla_{h}\tilde{\phi}^{m}\|_{2}^{2}\Bigr{)}^{\frac{1}{2}}\leq
C(\Delta t^{2}+h^{2}),$ (4.21)

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

###### Proof.

A careful consistency analysis indicates the following truncation error
estimate:

$\begin{split}\frac{\frac{3}{2}\Phi_{N}^{n+1}-2\Phi_{N}^{n}+\frac{1}{2}\Phi_{N}^{n-1}}{\Delta
t}=&\Delta_{h}\Bigl{(}-\frac{8}{3}((\Phi_{N}^{n+1})^{-9}-(\Phi_{N}^{n+1})^{-3})+\frac{8}{3}A_{0}(\Phi_{N}^{n+1}-\check{\Phi}_{N}^{n+1})\\\
&-\varepsilon^{2}\Delta_{h}\Phi_{N}^{n+1}-A\Delta
t\Delta_{h}(\Phi_{N}^{n+1}-\Phi_{N}^{n})\Bigr{)}+\tau^{n},\end{split}$ (4.22)

with $\check{\Phi}_{N}^{n}=2\Phi_{N}^{n}-\Phi_{N}^{n-1}$,
$\|\tau^{n}\|_{-1,h}\leq C(\Delta t^{2}+h^{2})$. In turn, subtracting the
numerical scheme (4.1) from (4.22) gives

$\begin{split}\frac{\frac{3}{2}\tilde{\phi}^{n+1}-2\tilde{\phi}^{n}+\frac{1}{2}\tilde{\phi}^{n-1}}{\Delta
t}&=\Delta_{h}\Bigl{(}-\frac{8}{3}((\Phi_{N}^{n+1})^{-9}-(\Phi_{N}^{n+1})^{-3}-A_{0}\Phi_{N}^{n+1})\\\
&\quad+\frac{8}{3}((\phi^{n+1})^{-9}-(\phi^{n+1})^{-3}-A_{0}\phi^{n+1})-\frac{8}{3}A_{0}\tilde{\check{\phi}}^{n+1}\\\
&\quad-\varepsilon^{2}\Delta_{h}\tilde{\phi}^{n+1}-A\Delta
t\Delta_{h}(\tilde{\phi}^{n+1}-\tilde{\phi}^{n})\Bigr{)}+\tau^{n},\end{split}$
(4.23)

with $\tilde{\check{\phi}}^{n+1}=2\tilde{\phi}^{n}-\tilde{\phi}^{n-1}$. Taking
the discrete inner product of (4.23) with
$2(-\Delta_{h})^{-1}\tilde{\phi}^{n+1}$ yields

$\displaystyle\left\langle
3\tilde{\phi}^{n+1}-4\tilde{\phi}^{n}+\tilde{\phi}^{n-1},\tilde{\phi}^{n+1}\right\rangle_{-1,h}-2\varepsilon^{2}\Delta
t\left\langle\tilde{\phi}^{n+1},\Delta_{h}\tilde{\phi}^{n+1}\right\rangle_{\Omega}-2A\Delta
t\langle\Delta_{h}(\tilde{\phi}^{n+1}-\tilde{\phi}^{n}),\tilde{\phi}^{n+1}\rangle_{\Omega}$
$\displaystyle\quad+\frac{16}{3}\Delta
t\left\langle-((\Phi_{N}^{n+1})^{-9}-(\Phi_{N}^{n+1})^{-3}-A_{0}\Phi_{N}^{n+1})+((\phi^{n+1})^{-9}-(\phi^{n+1})^{-3}-A_{0}\phi^{n+1}),\tilde{\phi}^{n+1}\right\rangle_{\Omega}$
$\displaystyle\qquad=\frac{16}{3}A_{0}\Delta
t\left\langle\tilde{\check{\phi}}^{n+1},\tilde{\phi}^{n+1}\right\rangle_{\Omega}+2\Delta
t\langle\tau^{n},\tilde{\phi}^{n+1}\rangle_{-1,h}.$ (4.24)

For the temporal derivative stencil, the following identity is valid:

$\begin{split}\left\langle
3\tilde{\phi}^{n+1}-4\tilde{\phi}^{n}+\tilde{\phi}^{n-1},\tilde{\phi}^{n+1}\right\rangle_{-1,h}=&\frac{1}{2}\Bigl{(}\|\tilde{\phi}^{n+1}\|_{-1,h}^{2}-\|\tilde{\phi}^{n}\|_{-1,h}^{2}\\\
&+\|2\tilde{\phi}^{n+1}-\tilde{\phi}^{n}\|_{-1,h}^{2}-\|2\tilde{\phi}^{n}-\tilde{\phi}^{n-1}\|_{-1,h}^{2}\\\
&+\|\tilde{\phi}^{n+1}-2\tilde{\phi}^{n}+\tilde{\phi}^{n-1}\|_{-1,h}^{2}\Bigr{)}.\end{split}$
(4.25)

The estimate for the surface diffusion term is straightforward:

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

For the nonlinear inner product, the convex nature of
$\frac{1}{8}\phi^{-8}-\frac{1}{2}\phi^{-2}+\frac{1}{2}A_{0}\phi^{2}$ (at a
point-wise level) yields the following result:

$\displaystyle\left\langle-((\Phi_{N}^{n+1})^{-9}-(\Phi_{N}^{n+1})^{-3}-A_{0}\Phi_{N}^{n+1})+((\phi^{n+1})^{-9}-(\phi^{n+1})^{-3}-A_{0}\phi^{n+1}),\tilde{\phi}^{n+1}\right\rangle_{\Omega}\geq
0.$ (4.27)

The term associated with the truncation error can be controlled by the
standard Cauchy inequality:

$2\langle\tau^{n},\tilde{\phi}^{n+1}\rangle_{-1,h}\leq
2\|\tau^{n}\|_{-1,h}\cdot\|\tilde{\phi}^{n+1}\|_{-1,h}\leq\|\tau^{n}\|_{-1,h}^{2}+\|\tilde{\phi}^{n+1}\|_{-1,h}^{2}.$
(4.28)

For the concave expansive error term, a similar inequality is available:

$\displaystyle\frac{16}{3}A_{0}\left\langle\tilde{\check{\phi}}^{n+1},\tilde{\phi}^{n+1}\right\rangle_{\Omega}$
$\displaystyle\leq$
$\displaystyle\frac{16}{3}A_{0}\|\tilde{\check{\phi}}^{n+1}\|_{-1,h}\|\nabla_{h}\tilde{\phi}^{n+1}\|_{2}\leq\frac{64}{9}A_{0}^{2}\varepsilon^{-2}\|\tilde{\check{\phi}}^{n+1}\|_{-1,h}^{2}+\varepsilon^{2}\|\nabla_{h}\tilde{\phi}^{n+1}\|_{2}$
(4.29) $\displaystyle\leq$
$\displaystyle\frac{64}{9}A_{0}^{2}\varepsilon^{-2}(8\|\tilde{\phi}^{n}\|_{-1,h}^{2}+2\|\tilde{\phi}^{n-1}\|_{-1,h}^{2})+\varepsilon^{2}\|\nabla_{h}\tilde{\phi}^{n+1}\|_{2}.$

In addition, the following identity could be derived for the artificial
diffusion term:

$\begin{split}-2\langle\Delta_{h}(\tilde{\phi}^{n+1}-\tilde{\phi}^{n}),\tilde{\phi}^{n+1}\rangle_{\Omega}=&2\langle\nabla_{h}(\tilde{\phi}^{n+1}-\tilde{\phi}^{n}),\nabla_{h}\tilde{\phi}^{n+1}\rangle_{\Omega}\\\
=&\|\nabla_{h}\tilde{\phi}^{n+1}\|_{2}^{2}-\|\nabla_{h}\tilde{\phi}^{n}\|_{2}^{2}+\|\nabla_{h}(\tilde{\phi}^{n+1}-\tilde{\phi}^{n})\|_{2}^{2}.\end{split}$
(4.30)

Subsequently, a substitution of (4.25) – (4.30) into (4.24) yields

$\begin{split}&\|\tilde{\phi}^{n+1}\|_{-1,h}^{2}-\|\tilde{\phi}^{n}\|_{-1,h}^{2}+\|2\tilde{\phi}^{n+1}-\tilde{\phi}^{n}\|_{-1,h}^{2}-\|2\tilde{\phi}^{n}-\tilde{\phi}^{n-1}\|_{-1,h}^{2}\\\
&+A\Delta
t^{2}(\|\nabla_{h}\tilde{\phi}^{n+1}\|_{2}^{2}-\|\nabla_{h}\tilde{\phi}^{n}\|_{2}^{2})+\varepsilon^{2}\Delta
t\|\nabla_{h}\tilde{\phi}^{n+1}\|_{2}^{2}\\\
\leq&\frac{128}{9}A_{0}^{2}\varepsilon^{-2}\Delta
t(4\|\tilde{\phi}^{n}\|_{-1,h}^{2}+\|\tilde{\phi}^{n-1}\|_{-1,h}^{2})+2\Delta
t\|\tau^{n}\|_{-1,h}^{2}.\end{split}$ (4.31)

Similar to the arguments in Remark 3.3, under a constraint for the time step
size, namely in the form of $\frac{128}{9}A_{0}^{2}\varepsilon^{-2}\Delta
t\leq\frac{1}{32}$, we are able to apply the discrete Gronwall inequality and
obtain the desired convergence estimate:

$\|\tilde{\phi}^{n+1}\|_{-1,h}+\Bigl{(}\varepsilon^{2}\Delta
t\sum_{k=0}^{n+1}\|\nabla_{h}\tilde{\phi}^{m}\|_{2}^{2}\Bigr{)}^{1/2}\leq
C(\Delta t^{2}+h^{2}),$ (4.32)

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

##  5 The numerical solver and the numerical results

A preconditioned steepest descent (PSD) iteration algorithm is used to
implement the proposed numerical schemes (3.1) and (4.1), following the
practical and theoretical framework in [18]. We give the details for the first
order scheme (3.1); the details for the second order one (4.1) will be
similar. It is clear that (3.1) can be recast as a minimization of the
discrete convex energy functional (3.3), which becomes the discrete variation
of (3.3) set equal to zero: $\mathcal{N}_{h}[\phi]=f^{n}$, where

$\mathcal{N}_{h}[\phi]:=\frac{1}{\Delta
t}(-\Delta_{h})^{-1}(\phi-\phi^{n})-\frac{8}{3}\phi^{-9}-\varepsilon^{2}\Delta_{h}\phi\quad\mbox{and}\quad
f^{n}:=-\frac{8}{3}(\phi^{n})^{-3}.$

The essential idea of the PSD solver is to use a linearized version of the
nonlinear operator as a pre-conditioner. Specifically, the preconditioner,
$\mathcal{L}_{h}:\mathring{\mathcal{C}}_{\rm
per}\to\mathring{\mathcal{C}}_{\rm per}$, is defined as

${\mathcal{L}}_{h}[\psi]:=\frac{1}{\Delta
t}(-\Delta_{h})^{-1}\psi+\psi-\varepsilon^{2}\Delta_{h}\psi,$

and is a positive, symmetric operator. Specifically, this “metric” is used to
find an appropriate search direction for the steepest descent solver [18].
Given the current iterate $\phi_{n}\in\mathcal{C}_{\rm per}$, we define the
following _search direction_ problem: find $d_{n}\in\mathcal{C}_{\rm per}$
such that

${\mathcal{L}}_{h}[d_{n}]=r_{n}-\overline{r_{n}},\quad
r_{n}:=f^{n}-\mathcal{N}_{h}[\phi_{n}],$

where $r_{n}$ is the nonlinear residual of the $n^{\rm th}$ iterate
$\phi_{n}$. Of course, this equation can be efficiently solved using the Fast
Fourier Transform (FFT). Subsequently, the next iterate is obtained as

$\phi_{n+1}:=\phi_{n}+\alpha_{n}d_{n},$ (5.1)

where $\alpha_{n}\in\mathbb{R}$ is the unique solution to the steepest descent
line minimization problem

$\alpha_{n}:=\operatorname*{argmin}_{\alpha\in\mathbb{R}}{\cal
J}[\phi_{n}+\alpha
d_{n}]=\operatorname*{argzero}_{\alpha\in\mathbb{R}}\left\langle\mathcal{N}_{h}[\phi_{n}+\alpha
d_{n}]-f^{n},d_{n}\right\rangle_{\Omega}.$ (5.2)

Following similar techniques reported in [18], a theoretical analysis ensures
a geometric convergence of the iteration sequence. Also see [6, 10, 17, 19]
for the applications of the PSD solver to various gradient flow models.

######  Remark 5.1.

We observe that the PSD method can be viewed as a quasi-Newton method, with an
orthogonalization (line search) step. Indeed, $\mathcal{L}_{h}$ may be viewed
as an approximation of the Jacobian. To fit more neatly into the framework of
a traditional quasi-Newton method, one could just take step size equal to 1,
so that the correction is just $\phi_{n+1}:=\phi_{n}+d_{n}$. Alternatively,
one can just use quadratic line search methods to obtain an approximation of
$\alpha_{n}$, call it $\alpha_{n}^{\rm q}$, to obtain a sufficiently good
approximation that one can still prove a geometric convergence rate that is
independent of $N$.

###  5.1 Convergence test for the numerical schemes

In this subsection we perform a numerical accuracy check for the proposed
numerical schemes (3.1) and (4.1). The computational domain is chosen as
$\Omega=(0,1)^{2}$, and the exact profile for the phase variable is set to be

$\Phi(x,y,t)=1+\frac{1}{2\pi}\sin(2\pi x)\cos(2\pi y)\cos(t).$ (5.3)

In particular, we see that the exact profile is positive at a point-wise
level, so that the computation will not cause any singularity issue. To make
$\Phi$ satisfy the original PDE (1.2), we have to add an artificial, time-
dependent forcing term. Then the proposed scheme, either the first order one
(3.1), or the second order algorithm (4.1), can be implemented to solve for
(1.2).

In the accuracy check for the first order scheme (3.1), we fix the spatial
resolution as $N=256$ (with $h=\frac{1}{256}$), so that the spatial numerical
error is negligible. The final time is set as $T=1$, and the surface diffusion
parameter is taken as $\varepsilon=0.5$. Naturally, a sequence of time step
sizes are taken as $\Delta t=\frac{T}{N_{T}}$, with $N_{T}=100:100:1000$. The
expected temporal numerical accuracy assumption $e=C\Delta t$ indicates that
$\ln|e|=\ln(CT)-\ln N_{T}$, so that we plot $\ln|e|$ vs. $\ln N_{T}$ to
demonstrate the temporal convergence order. The fitted line displayed in
Figure 1 shows an approximate slope of -0.9704, which in turn verifies a nice
first order temporal convergence order, in both the discrete $\ell^{2}$ and
$\ell^{\infty}$ norms.

Figure 1: The discrete $\ell^{2}$ and $\ell^{\infty}$ numerical errors vs.
temporal resolution $N_{T}$ for $N_{T}=100:100:1000$, with a spatial
resolution $N=256$. The numerical results are obtained by the computation
using the first order scheme (3.1). The surface diffusion parameter is taken
to be $\varepsilon=0.5$. The data lie roughly on curves $CN_{T}^{-1}$, for
appropriate choices of $C$, confirming the full first-order accuracy of the
scheme.

In the accuracy test for the second order scheme (4.1), we set the time size
as $\Delta t=\frac{1}{2}h$, with $h=\frac{1}{N}$, so that the second order
accuracy in both time and space could be confirmed. Again, the final time set
as by $T=1$, and the surface diffusion parameter is given by
$\varepsilon=0.5$. The artificial regularization parameters are taken as
$A_{0}=A_{0}^{\star}=\frac{9}{5}(\frac{2}{15})^{\frac{2}{3}}$,
$A=\frac{4}{9}A_{0}^{2}$, so that the convexity condition and stability
requirement are satisfied. A sequence of spatial resolutions are taken as
$N=48:16:192$. The expected temporal numerical accuracy assumption $e=C(\Delta
t^{2}+h^{2})$ indicates that $\ln|e|=\ln C-2\ln N$, so that we plot $\ln|e|$
vs. $\ln N$ to demonstrate the temporal convergence order. The fitted line
displayed in Figure 2 shows an approximate slope of -1.9974, which in turn
verifies a perfect second convergence order, in both time and space.

Figure 2: The discrete $\ell^{2}$ and $\ell^{\infty}$ numerical errors vs.
spatial resolution $N$ for $N=48:16:192$, and the time step size is set as
$\Delta t=\frac{1}{2}h$. The numerical results are obtained by the computation
using the second order scheme (4.1). The surface diffusion parameter is taken
to be $\varepsilon=0.5$. The data lie roughly on curves $CN^{-2}$, for
appropriate choices of $C$, confirming the full second-order accuracy of the
scheme, in both time and space.

###  5.2 Numerical simulation of coarsening process

The nonlinear term in (1.2) associated with the Lennard-Jones-type potential
gives preference to rotationally invariant patterns. In this subsection, we
perform a two-dimensional numerical simulation showing the coarsening process.
The computational domain is set as $\Omega=(0,L)^{2}$, with $L=12.8$, and the
interface width parameter is taken as $\varepsilon=0.02$. The initial data are
given by

$\phi^{0}_{i,j}=2+0.1\cdot(2r_{i,j}-1),\quad\mbox{$r_{i,j}$ are uniformly
distributed random numbers in $[0,1]$}.$ (5.4)

The second order numerical scheme (4.1) is implemented for this simulation.
For the temporal step size $\Delta t$, we use increasing values of $\Delta t$
in the time evolution: $\Delta t=0.001$ on the time interval $[0,100]$,
$\Delta t=0.004$ on the time interval $[100,500]$, $\Delta t=0.008$ on the
time interval $[500,2000]$ and $\Delta t=0.02$ on the time interval
$[2000,6000]$. Whenever a new time step size is applied, we initiate the two-
step numerical scheme by taking $\phi^{-1}=\phi^{0}$, with the initial data
$\phi^{0}$ given by the final time output of the last time period. The time
snapshots of the evolution by using $\varepsilon=0.02$ are presented in Figure
3, with significant coarsening observed in the system. At early times many
small hills (yellow) are present, with flat base (blue). At the final time,
$t=6000$, a single hill structure emerges, and further coarsening is not
possible.

$t=6,20$

$t=40,60$

$t=100,200$

$t=300,400$

$t=500,900$

$t=2000,6000$

Figure 3: (Color online.) Snapshots of the computed height function $\phi$ at
the indicated times for the parameters $L=12.8$, $\varepsilon=0.02$. Finally,
there is a single droplet structure.

The long time characteristics of the solution, especially the energy decay
rate, are of interest to material scientists. Figure 4 presents the log-log
plot for the energy versus time, with the given physical parameters. Recall
that, at the space-discrete level, the energy, $F_{h}$ is defined via (3.2).
The detailed scaling “exponent” is obtained using least squares fits of the
computed data up to time $t=100$. A clear observation of the $a_{e}t^{b_{e}}$
scaling law can be made, with $a_{e}=59.8167$, $b_{e}=-0.1951$. For the Cahn-
Hilliard flow with polynomial approximation of the double-well energy
potential, various numerical experiments have indicated an approximately
$t^{-1/3}$ energy dissipation law [9, 11]. For the droplet liquid film
equation (1.2), with Lennard-Jones-type energy potential included in (1.1),
this numerical evidence has implied a different energy dissipation scaling
index, with $b_{e}=-0.1951$, in comparison with an approximate $t^{-1/3}$
scaling law for the standard Cahn-Hilliard model.

Figure 4: Log-log plot of the temporal evolution the energy $F_{h}$ for
$\varepsilon=0.02$. The energy decreases like $a_{e}t^{b_{e}}$ until
saturation. The red lines represent the energy plot obtained by the
simulations, while the straight lines are obtained by least squares
approximations to the energy data. The least squares fit is only taken for the
linear part of the calculated data, only up to about time $t=100$. The fitted
line has the form $a_{e}t^{b_{e}}$, with $a_{e}=59.8167$, $b_{e}=-0.1951$.

######  Remark 5.2.

The coarsening dynamics problem usually is a long time process. To improve the
computational efficiency, some adaptive time stepping strategies have been
extensively applied in such a long-time simulation effort; our simulation has
also used variable time step size method, as outlined above. There have been
some theoretical works for the time-adaptive methods in the computation of
coarsening process, such as the ones for the Cahn-Hilliard model [8] and
epitaxial thin film model [40]. The corresponding analysis for the time
adaptive method applied to the droplet liquid film coarsening dynamics will be
left to the future works.

######  Remark 5.3.

For the physical energy given by (1.1), the final steady state solution
depends on $\varepsilon$ and the the mass average of $\phi$, which is a
conserved quantity. It is well-known that, if the mass average is less than or
equal to 1, the steady state solution turns out to be a trivial constant
profile, due to its minimum energy constrained by the mass average. On the
other hand, if the mass average is grater than 1, the time coarsening process
and the final profile have interesting structures, as demonstrated in our
reported numerical results. Under a circumstance that the initial data contain
certain values near zero at certain region, the singular and convex nature of
the Lennard-Jones-type energy potential (1.1) will push the solution away from
the singular value of 0; in this case, if the mass average is greater than 1,
one could still observe the interesting coarsening process, similar to our
numerical simulation result.

##  6 Concluding remarks

Two finite difference numerical schemes, including both the first and second
order accurate algorithms, are proposed and analyzed for the droplet liquid
film model, in which a singular Leonard-Jones energy potential is involved. In
the first order scheme, the convex potential and the surface diffusion terms
are implicitly, while the concave potential term is updated explicitly. A
theoretical justification of the unique solvability and positivity-preserving
property is provided, so that a singularity is avoided in the scheme. Such an
analysis is based on the subtle fact that the numerical solution is equivalent
to a minimization of a convex numerical energy, and the singular nature of the
Leonard-Jones potential term around the value of 0 prevents the numerical
solution reaching such singular value. The energy stability of the numerical
schemes comes from the convex-concave decomposition of the physical energy. In
the second order scheme, an alternate convex-concave decomposition is
formulated so that the concave part corresponds to a quadratic energy. In
turn, the BDF temporal stencil is applied, the combined Leonard-Jones
potential term is treated implicitly, the concave expansive part the is
approximated by a second order Adams-Bashforth explicit extrapolation, and an
artificial Douglas-Dupont regularization term is added to ensure the energy
stability. The unique solvability and the positivity-preserving property for
the second order scheme could be similarly established, in which the convexity
of the $\phi^{-8}$ term has played an important role. The preconditioned
steepest descent (PSD) iteration solver has been applied to implement the
nonlinear numerical schemes, due to the fact the implicit part turns out to be
the gradient of a strictly convex energy functional. A few numerical results
are also presented in this article, which demonstrates the robustness of the
proposed numerical scheme.

## Acknowledgements

This work is supported in part by the National Science Foundation (USA) grants
NSF DMS-2012669 (C. Wang), and NSF DMS-1719854, DMS-2012634 (S. Wise), the
National Natural Science Foundation of China (NSFC) 11871105, 11571045 and
Science Challenge Project TZ2018002 (Z. Zhang).

## References

  * [1] O. Alexander and S.G. Bankoff.  Dewetting of a heated surface by an evaporating liquid film under conjoining/disjoining pressures.  J. Colloid Interface Sci., 218:152–166, 1999.
  * [2] O. Alexander, S.H. Davis, and S.G. Bankoff.  Long-scale evolution of thin liquid films.  Rev. Mod. Phys., 69:931–980, 1997.
  * [3] K.S. Anna, Z.Y. Zhang, and J.F. Wendelken.  Island diffusion and coarsening on metal (100) surfaces.  Phys. Rev. Lett., 79:3210–3213, 1997.
  * [4] J. Becker, G. Grun, R. Seeman, H. Mantz, K. Jacobs, K. R. Mecke, and R. Blossey.  Complex dewetting scenarios captured by thin-film models.  Nature Materials, 2:59–63, 2003.
  * [5] A.L. Bertozzi, G. Grun, and T. P. Witelski.  Dewetting films: bifurcations and concentrations, institute of physics publishing.  Institute of physics publishing, 14:1569–1592, 2001.
  * [6] W. Chen, C. Wang, S. Wang, X. Wang, and S. Wise.  Energy stable numerical schemes for a ternary Cahn-Hilliard system.  J. Sci. Comput., 84:27, 2020.
  * [7] W. Chen, C. Wang, X. Wang, and S.M. Wise.  Positivity-preserving, energy stable numerical schemes for the Cahn-Hilliard equation with logarithmic potential.  J. Comput. Phys.: X, 3:100031, 2019.
  * [8] W. Chen, X. Wang, Y. Yan, and Z. Zhang.  A second order bdf numerical scheme with variable steps for the Cahn-Hilliard equation.  SIAM J. Numer. Anal., 57:495–525, 2019.
  * [9] K. Cheng, W. Feng, C. Wang, and S.M. Wise.  An energy stable fourth order finite difference scheme for the Cahn-Hilliard equation.  J. Comput. Appl. Math., 362:574–595, 2019.
  * [10] K. Cheng, C. Wang, and S.M. Wise.  An energy stable Fourier pseudo-spectral numerical scheme for the square phase field crystal equation.  Commun. Comput. Phys., 26:1335–1364, 2019.
  * [11] K. Cheng, C. Wang, S.M. Wise, and X. Yue.  A second-order, weakly energy-stable pseudo-spectral scheme for the Cahn-Hilliard equation and its solution by the homogeneous linear iteration method.  J. Sci. Comput., 69:1083–1114, 2016.
  * [12] L. Cherfils, A. Miranville, and S. Zelik.  The Cahn-Hilliard equation with logarithmic potentials.  Milan J. Math., 79:561–596, 2011.
  * [13] 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.  Commun. Math. Sci., 17:921–939, 2019.
  * [14] 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.  Commun. Comput. Phys., 28:967–998, 2020.
  * [15] Q. Du, L. Ju, X. Li, and Z. Qiao.  Maximum principle preserving exponential time differencing schemes for the nonlocal Allen-Cahn equation.  SIAM J. Numer. Anal., 57:875–898, 2019.
  * [16] Q. Du, L. Ju, X. Li, and Z. Qiao.  Maximum bound principles for a class pf semi-linear parabolic equations and exponential time differencing schemes.  SIAM Rev., 2020.  Accepted and in press: https://arxiv.org/abs/2015.114.65.
  * [17] W. Feng, Z. Guan, J.S. Lowengrub, C. Wang, S.M. Wise, and Y. Chen.  A uniquely solvable, energy stable numerical scheme for the functionalized Cahn-Hilliard equation and its convergence analysis.  J. Sci. Comput., 76(3):1938–1967, 2018.
  * [18] W. Feng, A.J. Salgado, C. Wang, and S.M. Wise.  Preconditioned steepest descent methods for some nonlinear elliptic equations involving p-Laplacian terms.  J. Comput. Phys., 334:45–67, 2017.
  * [19] W. Feng, C. Wang, S.M. Wise, and Z. Zhang.  A second-order energy stable Backward Differentiation Formula method for the epitaxial thin film equation with slope selection.  Numer. Methods Partial Differ. Equ., 34(6):1975–2007, 2018.
  * [20] A. Giorgini, M. Grasselli, and A. Miranville.  The Cahn-Hiliard-Oono equation with singular potential.  Math. Models Methods Appl. Sci., 27(13):2485–2510, 2017.
  * [21] K.B. Glasner.  Ostwald ripening in thin film equations.  SIAM J. Appl. Math., 69:473–493, 2008.
  * [22] K.B. Glasner and T. P. Witelski.  Coarsening dynamics of dewetting films.  Phys. Rev. E., 67:1–12, 2003.
  * [23] J. Guo, C. Wang, S.M. Wise, and X. Yue.  An $H^{2}$ convergence of a second-order convex-splitting, finite difference scheme for the three-dimensional Cahn-Hilliard equation.  Commun. Math. Sci., 14:489–515, 2016.
  * [24] Y. Hao, Q. Huang, and C. Wang.  A third order BDF energy stable linear scheme for the no-slope-selection thin film model.  Commun. Comput. Phys., 2020.  Accepted and in press.
  * [25] J.N. Israelachvili.  Intermolecular and Surface Forces.  Academic Press, New York, NY, 1992.
  * [26] J. Kim.  Numerical simulations of phase separation dynamics in a water-oil-surfactant system.  J. Colloid Interface Sci., 303:272–279, 2006.
  * [27] R.V. Kohn and F. Otto.  Upper bound on coarsening rate.  Commun. Math. Phys., 229:375–395, 2002.
  * [28] R.V. Kohn and X.D. Yan.  Coarsening rates for models of multicomponent phase separation.  Interfaces and Free Boundaries, 6:135–149, 2004.
  * [29] D. Li and Z. Qiao.  On second order semi-implicit Fourier spectral methods for 2D Cahn-Hilliard equations.  J. Sci. Comput., 70:301–341, 2017.
  * [30] D. Li and Z. Qiao.  On the stabilization size of semi-implicit Fourier-spectral methods for 3D Cahn-Hilliard equations.  Commun. Math. Sci., 15:1489–1506, 2017.
  * [31] W. Li, W. Chen, C. Wang, Y. Yan, and R. He.  A second order energy stable linear scheme for a thin film model without slope selection.  J. Sci. Comput., 76(3):1905–1937, 2018.
  * [32] R. Limary and P. F. Green.  Late-stage coarsening of an unstable structured liquid film.  Phys. Rev. E., 66:1–6, 2002.
  * [33] R. Limary and P. F. Green.  Dynamics of droplets on the surface of a structured fluid film: Late-stage coarsening.  Langmuir, 19:2419–2424, 2003.
  * [34] X. Meng, Z. Qiao, C. Wang, and Z. Zhang.  Artificial regularization parameter analysis for the no-slope-selection epitaxial thin film model.  CSIAM Trans. Appl. Math., 1:441–462, 2020.
  * [35] A. Miranville.  On a phase-field model with a logarithmic nonlinearity.  Appl. Math., 57:215–229, 2012.
  * [36] A. Oron, S. H. Davis, and S. G. Bankoff.  Long-scale evolution of thin liquid films.  Rev. Mod. Phys., 69:932–977, 1997.
  * [37] F. Otto, T. Rump, and D. Slepcev.  Coarsening rates for a droplet model: Rigorous upper bounds.  SIAM J. Math. Anal., 38(2):503, 2006.
  * [38] R.L. Pbgo and A.M. Stuart.  Front migration in the nonlinear Cahn-Hilliard equation.  Proc. R. Soc. Lond., 422:261–278, 1989.
  * [39] L.M. Pismen.  Spinodal dewetting in a volatile liquid film.  Phys. Rev. E., 70:1–9, 2004.
  * [40] Z. Qiao, Z. Zhang, and T. Tang.  An adaptive time-stepping strategy for the molecular beam epitaxy models.  SIAM J. Sci. Comput., 33:1395–1414, 2011.
  * [41] S. Ravati and B.D. Favis.  Interfacial coarsening of ternary polymer blends with partial and complete wetting structures.  Polymer, 54:6739–6751, 2013.
  * [42] U. Thiele, M. G. Velarde, and K. NeuKer.  Dewetting: film rupture by nucleation in the spinodal regime.  Phys. Rev. Lett., 87:1–4, 2001.
  * [43] C. Wang and S.M. Wise.  An energy stable and convergent finite-difference scheme for the modified phase field crystal equation.  SIAM J. Numer. Anal., 49:945–969, 2011.
  * [44] S.M. Wise.  Unconditionally stable finite difference, nonlinear multigrid simulation of the Cahn-Hilliard-Hele-Shaw system of equations.  J. Sci. Comput., 44:38–68, 2010.
  * [45] S.M. Wise, C. Wang, and J. Lowengrub.  An energy stable and convergent finite-difference scheme for the phase field crystal equation.  SIAM J. Numer. Anal., 47:2269–2288, 2009.
  * [46] Y. Yan, W. Chen, C. Wang, and S.M. Wise.  A second-order energy stable BDF numerical scheme for the Cahn-Hilliard equation.  Commun. Comput. Phys., 23:572–602, 2018.