# Numerical Investigation of the Sharp-Interface Limit of the Navier-
Stokes–Cahn-Hilliard Equations

T.H.B. Demont  S.K.F. Stoter  E.H. van Brummelen  Eindhoven University of
Technology, Department of Mechanical Engineering, P.O. Box 513, 5600 MB
Eindhoven, The Netherlands

###### Abstract

In this article, we study the behavior of the Abels-Garcke-Grün Navier-
Stokes–Cahn-Hilliard diffuse-interface model for binary-fluid flows, as the
diffuse-interface thickness passes to zero. We consider this so-called sharp-
interface limit in the setting of the classical oscillating-droplet problem.
To provide reference limit solutions, we derive new analytical expressions for
small-amplitude oscillations of a viscous droplet in a viscous ambient fluid
in two dimensions. We probe the sharp-interface limit of the Navier-
Stokes–Cahn-Hilliard equations by means of an adaptive finite-element method,
in which the refinements are guided by an a-posteriori error-estimation
procedure. The adaptive-refinement procedure enables us to consider diffuse-
interface thicknesses that are significantly smaller than other relevant
length scales in the droplet-oscillation problem, allowing an exploration of
the asymptotic regime. For two distinct modes of oscillation, we determine the
optimal scaling relation between the diffuse-interface thickness parameter and
the mobility parameter. Additionally, we examine the effect of deviations from
the optimal scaling of the mobility parameter on the approach of the diffuse-
interface solution to the sharp-interface solution.

###### keywords:

Navier-Stokes–Cahn-Hilliard equations, Sharp-interface limit, Two-dimensional
oscillating droplet, Analytical solution , Adaptive finite-element methods

##  1 Introduction

Binary-fluid flows in which the two fluid components are separated by a
molecular transition layer are omnipresent in science and engineering.
Examples are inkjet printing and additive manufacturing. Mathematical-physical
models for binary-fluid flows generally fall under one of two categories,
namely sharp-interface or diffuse-interface models. In sharp-interface models,
the surface that separates the two fluid components is represented explicitly
by a manifold of co-dimension one. This manifold carries kinematic and dynamic
interface conditions, which act as boundary conditions on the initial
boundary-value problems of the two contiguous fluid components and, in
addition, determine the evolution of the manifold. Sharp-interface models are
therefore of free-boundary type. In diffuse-interface models, the interface
between the two fluid components is represented as a thin-but-finite
transition layer, in which the two components are mixed in a proportion that
varies continuously and monotonously between the two pure species across the
layer. The strength of diffuse-interface models lies in their intrinsic
ability to account for topological changes of the fluid-fluid interface due to
coalescence or break-up of droplets or wetting, i.e. the propagation of the
fluid-fluid front along a (possibly elastic) solid substrate [28, 19, 31, 34].

Diffuse-interface models for two immiscible incompressible fluid species are
generally described by the Navier-Stokes–Cahn-Hilliard (NSCH) equations. The
NSCH equations represent a class of models, of which various renditions have
been proposed over the last half century: by Hohenberg and Halperin in the
late 1970s [16], by Lowengrub and Truskinovsky in the late 1990s [24], by
Shokrpour et al. in 2018 [29] and by Abels, Garcke and Grün in 2012 [4]. In
this article, we focus on the latter model, in view of its thermodynamic
consistency and its consistent reduction to the underlying single-fluid
Navier–Stokes equations in the pure species setting.

NSCH models invariably contain three parameters related to the diffuse
interface, viz. an interface-thickness parameter, $\varepsilon$, a mobility
parameter, $m$, and a surface-tension parameter, $\sigma$. The interface-
thickness parameter represents the transverse length scale of the transition
layer between the two fluid components, and the transition layer collapses
(specifically, is supposed to collapse) onto a manifold of co-dimension one in
the so-called sharp-interface limit $\varepsilon\to+0$. The mobility parameter
is responsible for the rate at which phase-diffusion occurs in the vicinity of
the diffuse interface. In the phase-separated regime in which the NSCH
equations are typically applied as a binary-fluid model, the mobility
parameter is responsible for the rate at which the interface equilibrates. In
the mixture regime, it governs the dynamics of the Ostwald-ripening effect.
The surface-tension parameter controls the excess free energy
$\sigma_{\textsc{da}}$ of the diffuse-interface according to
$2\sqrt{2}\sigma=3\sigma_{\textsc{da}}$. It is to be noted that for the NSCH
equations this proportionality holds independent of $\varepsilon$, as opposed
to the Navier-Stokes–Korteweg equations.

Contemporary understanding of the sharp-interface limit of the NSCH equations
is incomplete. An overview of known results and open questions is provided in
[3, §4.3]. One prominent open question pertains to the appropriate scaling of
the mobility parameter in relation to the interface-thickness parameter, in
the sharp-interface limit. The limit solution of the NSCH equations depends on
the scaling $m:=m_{\varepsilon}$. Abels and Garcke [3, §4.1] establish that,
if $m_{\varepsilon}=\mathcal{O}(1)$ as $\varepsilon\rightarrow+0$, their NSCH
model converges to the nonclassical sharp-interface
Navier–Stokes/Mullins–Sekerka model; see also [18, §4]. If, on the other hand,
$m_{\varepsilon}$ vanishes suitably as $\varepsilon\to+0$, the classical
sharp-interface binary-fluid model is obtained, where the interface is
transported by the fluid velocity [24, 3, 35]. However, the decay of the
mobility cannot be too fast: if $m_{\varepsilon}=o(\varepsilon^{3})$ as
$\varepsilon\to+0$, the resulting limit solution of the NSCH model generally
violates the Young–Laplace condition on the pressure jump across the interface
[5]. These results suggest that $m_{\varepsilon}\propto{}\varepsilon^{a}$ with
$0<a\leq{}3$ as $\varepsilon\to+0$ represents a necessary and sufficient
condition to achieve a classical sharp-interface solution. Still, the details
of the approach of the diffuse-interface solution to the sharp-interface limit
solution for these various admissible scalings of the mobility are not
currently known, and different scalings have been proposed in the literature,
in particular in the context of numerical simulation approaches. In [11], the
scaling $m_{\varepsilon}\propto\varepsilon^{3}$ is considered, based on the
argument that this proportionality fixes the diffusive time scale and, thus,
the equilibration rate of the diffuse interface. This cubic scaling of the
mobility with respect to the interface thickness (in terms of their usual
dimensional forms) is also propounded in [21], supported by numerical
investigations. Based on partial matched-asymptotic-analysis arguments, Ref.
[25] concludes that $m\propto\varepsilon^{2}$ is the appropriate scaling.
However, because the matching procedure in this reference is incomplete, it is
unclear whether this scaling relation in fact represents a necessary or
sufficient condition. On the basis of a consideration of curvature-induced
expansion/contraction modes at the diffuse interface, it is argued in [18]
that $m_{\varepsilon}\propto\varepsilon^{a}$ with $1\leq{}a<2$. It is to be
noted that the aforementioned scalings of the mobility pertain to situations
without moving contact lines and topological changes; see, e.g. [35].

In this article, we address the open questions of the optimal scaling of the
mobility and the approach to the sharp-interface limit solution by
computational investigation of the behavior of the Abels–Garcke–Grün Navier-
Stokes–Cahn-Hilliard model for different interface dynamics and different
mobility parameters as it limits toward a sharp-interface description of a
two-dimensional oscillating droplet. To enable an exploration of the
asymptotic regime, we apply an adaptive finite-element method, in which the
adaptivity is guided by an a-posteriori error estimate; see [31, 11] for
details.

We conduct our analysis of the sharp-interface limit of the NSCH equations in
the context of the prototypical oscillating-droplet problem, in two
dimensions. Despite the fact that the oscillating-droplet problem is
classical, it appears that the two-dimensional setting has not been
extensively investigated, and that solutions of the two-dimensional problem
have not been reported in the literature. The investigation of the
oscillating-droplet problem dates back to Rayleigh, who presented the well-
known frequency of oscillation of an inviscid droplet in vacuo in 1879 [30].
This result was extended by Lamb in the 1930s to include the effect of an
inviscid ambient fluid [22]. In 1960, Reid generalized the theory of
oscillating droplets in vacuo by including the effect of viscosity [27]. A
complete theory, comprising solutions for small oscillations of a viscous
droplet in a viscous ambient fluid, was then finally presented by Miller and
Scriven in 1968 [26]. The aforementioned references however exclusively
consider the three-dimensional case and the results, especially those for the
viscous solutions, do not trivially extend to the two-dimensional case.
Clearly, the three-dimensional case is the practically relevant one, but the
two-dimensional case has raison d’être independently as a means of
verification for mathematical models and numerical methods. Our analysis of
the sharp-interface limit of the NSCH equations requires access to closed form
solutions of the sharp-interface model, on the one hand to provide initial and
boundary data for the NSCH equations, and on the other hand to systematically
determine the deviation of the diffuse-interface solution relative to the
sharp-interface solution. A secondary objective of this work is therefore to
establish closed-form expressions for small-amplitude oscillations of a
viscous droplet in a viscous ambient fluid. Our derivation follows that of
Miller and Scriven, but we deviate from their derivation by a more complete
elaboration of intermediate steps and assumptions and, in particular, an
explicit accounting of the complex-valued nature of the different fields, and
by presenting closed-form expressions of the final results.

The remainder of this article is structured as follows. In Section 2, we lay
out the Abels–Garcke–Grün Navier-Stokes–Cahn-Hilliard model equations, and the
coupled Navier–Stokes free-boundary problem that they should reduce to in the
sharp-interface limit. In Section 3, we derive a closed form expression for
the sharp-interface model corresponding to small-amplitude oscillations of a
viscous droplet in a viscous ambient fluid in two dimensions. We make use of
these expressions in Section 4, where we study the approach of the NSCH
solution to the sharp-interface solution in the limit $\varepsilon\to+0$, by
means of systematic numerical experiments.

##  2 Governing equations

We consider a binary-fluid system, where both fluids are modeled as being
incompressible, isothermal, immiscible, and Newtonian with finite viscosity.
In accordance with the later focus on a submerged droplet, we denote one of
the fluids by d, for droplet, and the other by a, for ambient. Various
modeling frameworks for describing the fluid motion exist. These make use of
either a diffuse interface representation or a sharp-interface representation.
Figure 1 illustrates the different relevant domains and material parameters.
We consider in this work the incompressible NSCH model — specifically, the
model developed by Abels, Garcke, and Grün in [4] — in order to describe the
binary fluid dynamics. Motivation for this choice lies in its thermodynamic
consistency, incompressibility, and consistent reduction to the underlying
single-fluid Navier–Stokes equations in the pure species setting. Recently,
the well-posedness of the Abels–Garcke-Grün model in various settings has been
shown [14, 1, 2].

Figure 1: Schematic of the physical setting of a submerged immiscible fluid,
modeled with either a diffuse interface representation or a sharp-interface
representation.

###  2.1 Diffuse-interface representation

In diffuse-interface models, the two immiscible fluids are separated by a
layer of finite thickness constituted by a mixture of both fluids, reflecting
a gradual transition between fluid d and fluid a. We consider an open time
interval $(0,t_{\textrm{fin}})\subseteq\mathds{R}_{>0}$ and a spatial domain
corresponding to a simply connected time-independent subset
$\Omega\subseteq\mathds{R}^{d}$ ($d=2,3$). We make use of a Navier-
Stokes–Cahn-Hilliard type model that describes the evolution of a so-called
order parameter $\varphi\in[-1,1]$ representing pure species d and a when
$\phi=1$ and $\phi=-1$, respectively, and a mixture of both when
$\varphi\in(-1,1)$, in addition to the velocity and pressure of the mixture.
The NSCH model as presented by Abels, Garcke, and Grün is given by [4]:

$\displaystyle\partial_{t}\left(\rho\boldsymbol{u}\right)+\nabla\cdot\left(\rho\boldsymbol{u}\otimes\boldsymbol{u}\right)+\nabla\cdot\left(\boldsymbol{u}\otimes\boldsymbol{J}\right)+\nabla
p-\nabla\cdot\boldsymbol{\tau}-\nabla\cdot\boldsymbol{\zeta}$
$\displaystyle=0,$ (1a) $\displaystyle\nabla\cdot\boldsymbol{u}$
$\displaystyle=0,$ (1b)
$\displaystyle\partial_{t}\varphi+\nabla\cdot\left(\varphi\boldsymbol{u}\right)-\nabla\cdot\left(m\nabla\mu\right)$
$\displaystyle=0,$ (1c)
$\displaystyle\mu+\sigma\varepsilon\Delta\varphi-\frac{\sigma}{\varepsilon}\Psi^{\prime}$
$\displaystyle=0,$ (1d)

where the volume-averaged velocity $\boldsymbol{u}$, the pressure $p$, the
order parameter $\varphi$ and the chemical potential $\mu$ are the unknown
fields. The closure relations for the relative mass flux $\boldsymbol{J}$, the
viscous stress $\boldsymbol{\tau}$, the capillary stress $\boldsymbol{\zeta}$
and the mixture energy density $\Psi$ are given as:

$\displaystyle\boldsymbol{J}\coloneqq
m\frac{\rho_{\textsc{a}}-\rho_{\textsc{d}}}{2}\nabla\mu\,,$ (2a)
$\displaystyle\boldsymbol{\tau}\coloneqq\eta(\nabla\boldsymbol{u}+(\nabla\boldsymbol{u})^{T})\,,$
(2b)
$\displaystyle\boldsymbol{\zeta}\coloneqq-\sigma\varepsilon\nabla\varphi\otimes\nabla\varphi+\boldsymbol{I}\left(\frac{\sigma\varepsilon}{2}|\nabla\varphi|^{2}+\frac{\sigma}{\varepsilon}\Psi\right)\,,$
(2c)
$\displaystyle\Psi\left(\varphi\right)\coloneqq\frac{1}{4}\left(\varphi^{2}-1\right)^{2}\,.$
(2d)

The remaining parameters are material and model parameters. The model
parameters are the mobility parameter $m>0$ and the interface thickness
parameter $\varepsilon>0$, which affect the time and length scale of the
diffuse interface, respectively. The material parameters are $\sigma$, a
rescaling of the droplet-ambient surface tension $\sigma_{\textsc{da}}$
according to $2\sqrt{2}\sigma=3\sigma_{\textsc{da}}$, the mixture density
$\rho$, and the mixture viscosity $\eta$. The mixture density and viscosity
generally depend on $\varphi$. To ensure existence of a solution to the system
of equations, we must allow $\varphi$ to take on values outside of $[-1,1]$
[15]. We include a density extension that ensures positive densities even for
the nonphysical scenario $\varphi\notin[-1,1]$ [9]:

$\rho(\varphi)=\left\\{\begin{tabular}[]{ll}$\frac{1}{4}\rho_{\textsc{a}},$&$\varphi\leq-1-2\lambda\,,$\\\
$\frac{1}{4}\rho_{\textsc{a}}+\frac{1}{4}\rho_{\textsc{a}}\lambda^{-2}\left(1+2\lambda+\varphi\right)^{2},$&$\varphi\in(-1-2\lambda,-1-\lambda)\,,$\\\
$\frac{1+\varphi}{2}\rho_{\textsc{d}}+\frac{1-\varphi}{2}\rho_{\textsc{a}},$&$\varphi\in[-1-\lambda,1+\lambda]\,,$\\\
$\rho_{\textsc{d}}+\frac{3}{4}\rho_{\textsc{a}}-\frac{1}{4}\rho_{\textsc{a}}\lambda^{-2}\left(1+2\lambda-\varphi\right)^{2},$&$\varphi\in(1+\lambda,1+2\lambda)\,,$\\\
$\rho_{\textsc{d}}+\frac{3}{4}\rho_{\textsc{a}},$&$\varphi\geq
1+2\lambda\,,$\end{tabular}\right.$ (3)

where
$\lambda=\rho_{\textsc{a}}/\left(\rho_{\textsc{d}}-\rho_{\textsc{a}}\right)$.
For the viscosity interpolation, we apply the Arrhenius mixture-viscosity
model [6]:

$\log\eta(\varphi)=\frac{\left(1+\varphi\right)\Lambda\log\eta_{\textsc{d}}+\left(1-\varphi\right)\log\eta_{\textsc{a}}}{\left(1+\varphi\right)\Lambda+\left(1-\varphi\right)}\,,$
(4)

where
$\Lambda=\frac{\rho_{\textsc{d}}M_{\textsc{a}}}{\rho_{\textsc{a}}M_{\textsc{d}}}$
is the intrinsic volume ratio between the two fluids (with $M_{\textsc{a}}$
and $M_{\textsc{d}}$ their respective molar masses).

###### Remark 1

To eliminate the Ostwald-ripening effect in the pure species, a degenerate
dependence of the mobility on the phase field can be introduced, according to
$m(\varphi)\geq 0$ with inequality if and only if $|\varphi|<1$. However, as a
degenerate mobility introduces complications with regard to numerical-
approximation procedures [7], we opt for a constant mobility parameter.

###### Remark 2

In this work, we use a volume-fraction-based Arrhenius relation, i.e.
$\Lambda=1$. Because the denominator in Eq. 4 then reduces to a non-zero
constant, this choice eliminates singularities, and the mixture viscosity is
bounded away from zero in a finite interval including $[-1,1]$. See _Remark 2_
in [31] for further details.

###  2.2 Sharp-interface limit

As $\varepsilon\rightarrow+0$, the width of the diffuse interface in the NSCH
model reduces to zero. As pointed out in the introduction, the particular
model that arises in this limit depends on the scaling relation of the
mobility $m$. If the mobility also tends to zero appropriately, then the
following classical sharp-interface model is obtained:

$\displaystyle\rho_{i}\partial_{t}\boldsymbol{u}_{i}+\rho_{i}\left(\boldsymbol{u}_{i}\cdot\nabla\right)\boldsymbol{u}_{i}-\eta_{i}\Delta\boldsymbol{u}_{i}+\nabla
p_{i}$ $\displaystyle=0$ $\displaystyle\textrm{ in
}\Omega_{i}=\Omega_{i}(t)\,,$ (5a)
$\displaystyle\nabla\cdot\boldsymbol{u}_{i}$ $\displaystyle=0$
$\displaystyle\textrm{ in }\Omega_{i}\,,$ (5b)
$\displaystyle\boldsymbol{u}_{i}\cdot\boldsymbol{n}$
$\displaystyle=\mathcal{V}\quad$ $\displaystyle\textrm{ on
}\Gamma=\Gamma(t)\,,$ (5c) $[\\![\boldsymbol{u}\cdot\boldsymbol{t}_{j}]\\!]$
$\displaystyle=0\quad$ $\displaystyle\textrm{ on }\Gamma\text{, for
}j=1,\cdots,d-1\,,$ (5d)
$[\\![-(\nabla\boldsymbol{u}+(\nabla\boldsymbol{u})^{T})\boldsymbol{n}+p\boldsymbol{n}]\\!]$
$\displaystyle=\sigma_{\textsc{d}\textsc{a}}\kappa\boldsymbol{n}\quad$
$\displaystyle\textrm{ on }\Gamma\,,$ (5e)

for $i\in\\{\textsc{d},\textsc{a}\\}$, and with $\boldsymbol{n}$ the unit
normal vector on $\Gamma$ external to $\Omega_{\textsc{d}}$, $\mathcal{V}$ the
interface normal velocity, $\kappa$ the (additive) curvature of the interface,
and $[\\![\cdot]\\!]$ the interface jump operator
$\mbox{$[\\![g]\\!]$}=g|_{\textsc{a}}-g|_{\textsc{d}}$. We adhere to the
convention that curvature is negative if the center of the osculating circle
in the normal plane is located in the droplet domain.

As opposed to the NSCH model, the sharp-interface model (5) represents a set
of equations for each fluid species separately, complemented by appropriate
coupling conditions at the interface. The sharp-interface model represents a
free-boundary problem. The domains on which the various fields are defined
evolve in time, as reflected by the time-dependence of $\Omega_{i}(t)$ and
$\Gamma(t)$. There is an intrinsic coupling between the velocity field and the
evolution of $\Omega_{i}$ and $\Gamma$, according to (5c). Because this
equation holds on both sides of the interface and $\mathcal{V}$ is single-
valued, Eq. 5c implies
$\mbox{$[\\![\boldsymbol{u}\cdot\boldsymbol{n}]\\!]$}=0$.

##  3 Response of an oscillating droplet

With the objective of providing a reference solution for the sharp-interface
limit of the NSCH equations, we now consider solutions of the free-boundary
problem (5) corresponding to small perturbations of a circular droplet set in
an ambient fluid, in two dimensions. Denoting by $R_{0}$ the radius of the
droplet, one can verify that

$\displaystyle\Omega_{\textsc{d},0}$
$\displaystyle=\\{\boldsymbol{x}\in\mathbb{R}^{2}:|\boldsymbol{x}|<R_{0}\\}$
$\displaystyle\qquad\Omega_{\textsc{a},0}$
$\displaystyle=\mathbb{R}^{2}\setminus\overline{\Omega_{\textsc{d},0}}\,,$
(6a) $\displaystyle\boldsymbol{u}_{\textsc{d},0}$ $\displaystyle=0$
$\displaystyle\qquad\boldsymbol{u}_{\textsc{a},0}$ $\displaystyle=0\,,$ (6b)
$\displaystyle p_{\textsc{d},0}$
$\displaystyle=\sigma_{\textsc{d}\textsc{a}}/R_{0}\ $ $\displaystyle\qquad
p_{\textsc{a},0}$ $\displaystyle=0\,,$ (6c)

represents a stationary solution to (5). We will use (6) as a generating
solution. We consider perturbations of the solution (6) that are suitably
bounded and vanish toward infinity. The latter condition can be expressed as:

$\lim\limits_{|\boldsymbol{x}|\rightarrow\infty}\big{(}\boldsymbol{u},p\big{)}(\boldsymbol{x},t)=0\,.$
(7)

Our derivation of the natural response of such a droplet follows that of
Miller and Scriven [26], except that we provide a more complete elaboration of
intermediate steps and assumptions and, in particular, an explicit accounting
of the complex-valued nature of the different fields. The approach essentially
comprises four steps. First, we linearize the governing equations around the
generating solution (6), perturbed by a small deformation of the interface.
Second, we derive the general solutions corresponding to the natural response
in both domains separately. In the third step, we incorporate the interface
coupling conditions by constraining the free parameters in the general
solutions. Finally, the characteristic temporal response (frequency of
oscillation and rate of damping) of the assumed interface displacement, as
well as the corresponding shapes, follow from a solution-existence condition.

###  3.1 Formal linearization

We consider small-amplitude perturbations of the interface that are sinusoidal
along the circumference of the droplet. To facilitate the presentation, we
introduce polar coordinates $r\in\mathbb{R}_{\geq{}0}$ and $\theta\in[0,2\pi)$
and the coordinate transformation
$\boldsymbol{x}=(x_{1},x_{2})=r(\cos\theta,\sin\theta)$. We regard
perturbations of the interface $\Gamma_{0}=\partial\Omega_{\textsc{d},0}$
corresponding to the following parametrization:

$\Gamma_{\delta}(t)=\big{\\{}\boldsymbol{x}\in\mathbb{R}^{2}:\boldsymbol{x}=R_{\delta}(\theta,t)\,(\cos(\theta),\sin(\theta)),\theta\in[0,2\pi)\big{\\}}\,,$
(8)

where

$\displaystyle R_{\delta}(\theta,t)$
$\displaystyle=R_{0}+R_{0}\,\delta\left(\beta\cos(k\theta)+\sqrt{1-\beta^{2}}\,\sin(k\theta)\right)\cos(\nu
t)e^{-\alpha t}$ (9)
$\displaystyle=\Re\left(R_{0}+R_{0}\,\delta\left(\beta\cos(k\theta)+\sqrt{1-\beta^{2}}\,\sin(k\theta)\right)e^{-\gamma
t}\right)\,.$

The interface configuration (8)-(9) represents a damped oscillation of the
droplet with a mode-shape described by the mode number $k\in\mathbb{N}$, with
angular orientation dependent on $0\leq\beta\leq 1$, and with an amplitude
described by $\delta\ll{}1$ as the fraction of the droplet radius $R_{0}$. Our
interest is restricted to droplet configurations for which
$\operatorname{meas}(\Omega_{\textsc{d}})=\operatorname{meas}(\Omega_{\textsc{d},0})+\mathcal{O}(\delta^{2})$
and the barycenter of $\Omega_{\textsc{d}}$ is located at the origin. This
implies that $k\in\mathbb{N}_{\geq 2}$. The damping rate $\alpha\geq 0$ and
the frequency of oscillation $\nu>0$ implicitly depend on the mode number and
will follow from the subsequent analysis. The second expression in (9)
provides a representation of $R_{\delta}$ as the real part of a complex-valued
function, with $\gamma\coloneqq\alpha-i\nu$. This form enables us to condense
some of the expressions that appear in the sequel.

In conjunction with the interface configuration (8)-(9), we consider linear
asymptotic solutions of the sharp-interface problem of the form

$\big{(}\boldsymbol{u}_{i},p_{i},\boldsymbol{n},\boldsymbol{t},\kappa,\mathcal{V})=\big{(}\boldsymbol{u}_{i},p_{i},\boldsymbol{n},\boldsymbol{t},\kappa,\mathcal{V})_{0}+\delta\big{(}\boldsymbol{u}_{i},p_{i},\boldsymbol{n},\boldsymbol{t},\kappa,\mathcal{V})_{1}\,,$
(10)

i.e. functions conforming to (10) that satisfy (5) modulo terms of $o(\delta)$
as $\delta\to{}0$. Substituting (10) into the sharp-interface equations (5),
collecting terms of distinct orders in $\delta$, and noting that all terms of
$\mathcal{O}(1)$ vanish on account of the fact that the first term in (10)
represents a solution to (5), we obtain the following infinitesimal conditions
on the second term in (10):

$\displaystyle\rho\partial_{t}\boldsymbol{u}_{i,1}-\eta\Delta\boldsymbol{u}_{i,1}+\nabla
p_{i,1}$ $\displaystyle=0$ $\displaystyle\textrm{ in }\Omega_{i,0}\,,$ (11a)
$\displaystyle\nabla\cdot\boldsymbol{u}_{i,1}$ $\displaystyle=0$
$\displaystyle\textrm{ in }\Omega_{i,0}\,,$ (11b)
$\displaystyle\boldsymbol{u}_{i,1}\cdot\boldsymbol{n}_{0}$
$\displaystyle=\mathcal{V}_{1}$ $\displaystyle\textrm{ on }\Gamma_{0}\,,$
(11c) $[\\![\boldsymbol{u}_{1}\cdot\boldsymbol{t}_{0}]\\!]$ $\displaystyle=0$
$\displaystyle\textrm{ on }\Gamma_{0}\,,$ (11d)
$[\\![-\eta\left(\nabla\boldsymbol{u}_{1}+\left(\nabla\boldsymbol{u}_{1}\right)^{T}\right)\boldsymbol{n}_{0}+p_{1}\boldsymbol{n}_{0}]\\!]$
$\displaystyle=\sigma_{\textsc{d}\textsc{a}}\kappa_{1}\boldsymbol{n}_{0}\qquad$
$\displaystyle\textrm{ on }\Gamma_{0}\,,$ (11e)

for $i\in\\{\textsc{d},\textsc{a}\\}$. The non-linear advective term has
dropped since the only non-linear first-order perturbation terms are cross
terms between $\boldsymbol{u}_{i,1}$ and $\boldsymbol{u}_{i,0}=0$. Similarly,
only $\boldsymbol{n}_{0}$ appears in the first-order conditions (11), because
its multiplication with $\boldsymbol{u}_{i,1}$ is a second order term, its
dot-product with $\boldsymbol{u}_{i,0}$ vanishes, and
$\mbox{$[\\![p_{0}\boldsymbol{n}_{1}]\\!]$}=\sigma_{\textsc{d}\textsc{a}}R^{-1}_{0}\boldsymbol{n}_{1}$
cancels with the right-hand-side
$\sigma_{\textsc{d}\textsc{a}}\kappa_{0}\boldsymbol{n}_{1}$.

###### Remark 3

The first-order conditions are set on the stationary generating domains,
$\Omega_{i,0}$, and the generating interface, $\Gamma_{0}$. This is a
universal characteristic of linearizations of free-boundary problems. The
perturbation of the interface according to (9) appears implicitly in (11) in
the interface-velocity perturbation, $\mathcal{V}_{1}$, and the curvature
perturbation, $\kappa_{1}$.

The remainder of this section is devoted to finding general solutions to (11)
for different perturbation wave-numbers $k$. For purposes of readability,
henceforth we suppress the subscripts corresponding to the order of
perturbation.

###  3.2 First-order solutions in the droplet and ambient domain

Next, we derive the general first-order solutions $(\boldsymbol{u},p)_{i,1}$
in accordance with the differential equations (11a)–(11b) for both the droplet
d and ambient a domains. We proceed by deriving the vorticity equation
corresponding to (11a), which we solve by means of separation of variables.
The first-order velocity field is subsequently retrieved from the vorticity
solutions. The corresponding first-order pressure fields are derived as those
that yield balance of linear momentum.

####  3.2.1 Vorticity solution

The pressure may be eliminated from the governing equations by taking the curl
of (11a) and by using the identity $\nabla\times\nabla(\cdot)=0$. Introducing
the vorticity $\omega=\nabla\times\boldsymbol{u}$, we infer from (11a) that

$\displaystyle\rho\partial_{t}\left(\nabla\times\boldsymbol{u}\right)-\eta\Delta\left(\nabla\times\boldsymbol{u}\right)+\nabla\times\nabla{}p=\rho\partial_{t}\omega-\eta\Delta\omega=0\,.$
(12)

Let us note that in a two-dimensional setting, vorticity can be represented as
a scalar-valued field. The evolution equation for this scalar vorticity field
can be recognized as a diffusion equation.

To determine the general solution to (12), we assume the following separation
of variables form:

$\omega(r,\theta,t)\coloneqq\Psi(r)\,\Theta(\theta)\,T(t).$ (13)

Substitution in the polar coordinate representation of the diffusion equation
gives

$\displaystyle\Psi(r)\,\Theta(\theta)\,T^{\prime}(t)=\frac{\eta}{\rho
r^{2}}\Psi(r)\,\Theta^{\prime\prime}(\theta)\,T(t)+\frac{\eta}{\rho
r}\Psi^{\prime}(r)\,\Theta(\theta)\,T(t)+\frac{\eta}{\rho}\Psi^{\prime\prime}(r)\,\Theta(\theta)\,T(t)\,,$
(14)

where primes denote differentiation. The usual separation of variables
argument leads to

$\displaystyle
T^{\prime}(t)=-\frac{\eta}{\rho}m^{2}\,T(t)\quad\,\,\,m\in\mathbb{C}\,,$ (15a)
$\displaystyle\Theta^{\prime\prime}(\theta)=-n^{2}\,\Theta(\theta)\qquad
n\in\mathbb{C}\,,$ (15b) $\displaystyle
r^{2}\Psi^{\prime\prime}(r)+r\Psi^{\prime}(r)+\left(m^{2}r^{2}-n^{2}\right)\Psi(r)=0\,,$
(15c)

with $\mathbb{C}$ the set of complex numbers.

The general solutions of (15a) and (15b) consist of complex-valued exponential
functions, according to

$\displaystyle T(t)=c\,e^{-\frac{\eta}{\rho}m^{2}t}\,,$ (16a)
$\displaystyle\Theta(\theta)=c_{1}\,e^{int}+c_{2}e^{-int}\,.$ (16b)

From the periodicity of the droplet perturbations in the angular dependence,
conforming to (8), we infer that $n\in\mathbb{Z}_{\geq{}0}$. The arbitrary
constants $c_{1}$ and $c_{2}$ can then be selected such that (16b) reduces to
the sum of two real-valued trigonometric functions:

$\Theta(\theta)=C\cos(n\theta)+D\sin(n\theta)\qquad(n\in\mathbb{Z}_{\geq
0})\,,$ (17)

where $C,D\in\mathbb{R}$ are coefficients that determine the angular
orientation of the solution. Regarding (15c), we note that this equation
corresponds to Bessel’s equation with a complex-valued scaling
$m\in\mathbb{C}$. Solutions of (15c) therefore consist of extensions of Bessel
functions to the complex plane. Such extensions of Bessel functions are well
defined, by virtue of the fact that Bessel functions are analytic functions on
$\mathbb{R}$ and can hence be extended to analytic functions on $\mathbb{C}$
via their power-series expansion. The general solution of (15c) consists of a
linear combination of two Bessel functions of order $n\in\mathbb{Z}_{\geq 0}$.
For reasons that will become clear once we consider the boundary conditions,
we choose to work with a Bessel function of the first kind, $J_{n}$, and a
Hankel function of the second kind, $H^{(2)}_{n}$:

$\Psi(r)=Am^{2}J_{n}(mr)+Bm^{2}H^{(2)}_{n}(mr)\,,$ (18)

with $A,B\in\mathbb{C}$. Substitution of (16a), (17) and (18) into (13) gives
the general rotationally periodic solution of the vorticity equation (12):

$\omega(r,\theta,t)=\Big{[}Am^{2}J_{n}(mr)+Bm^{2}H^{(2)}_{n}(mr)\Big{]}\Big{[}C\cos(n\theta)+D\sin(n\theta)\Big{]}\,e^{-\frac{\eta}{\rho}m^{2}t}\,,$
(19)

for arbitrary $A,B\in\mathbb{C}$, $C,D\in\mathbb{R}$ and $m\in\mathbb{C}$.

####  3.2.2 Velocity solutions

To obtain the velocity fields from the general vorticity solution (19), we
introduce a stream function $\chi$ according to:

$\Delta\chi=-\omega\,.$ (20)

The velocity can be retrieved from this stream function as

$\boldsymbol{u}=\frac{1}{r}\frac{\partial}{\partial\theta}\chi\,\boldsymbol{e}_{r}-\frac{\partial}{\partial
r}\chi\,\boldsymbol{e}_{\theta}\,.$ (21)

Based on the expression for $\omega$ in (19), we anticipate that a particular
solution to (20) is of the form:

$\chi(r,\theta,t)=\Upsilon(r)\,\big{[}C\cos(n\theta)+D\sin(n\theta)\big{]}\,e^{-\frac{\eta}{\rho}m^{2}t}\,.$
(22)

Substitution of (22) into (20) leads to the following ordinary differential
equation for $\Upsilon$:

$\Upsilon^{\prime\prime}(r)+\frac{1}{r}\Upsilon^{\prime}(r)-\frac{n^{2}}{r^{2}}\Upsilon(r)=-Am^{2}J_{n}(mr)-Bm^{2}H^{(2)}_{n}(mr)\,.$
(23)

The general solution to the nonhomogeneous ordinary differential equation (23)
is given by:

$\Upsilon(r)=Er^{n}+Fr^{-n}+AJ_{n}(mr)+BH^{(2)}_{n}(mr)\,.$ (24)

The first and second term in (24) constitute the homogeneous part of the
solution. The third and fourth term represent the particular part. From (21)
it then follows that the $r$ and $\theta$ components of the velocity field are
given by:

$\displaystyle
u_{r}=\frac{1}{r}\frac{\partial}{\partial\theta}\chi=\frac{n}{r}\Upsilon(r)\,\big{[}D\cos(n\theta)-C\sin(n\theta)\big{]}\,e^{-\frac{\eta}{\rho}m^{2}t}\,,$
(25a) $\displaystyle u_{\theta}=-\frac{\partial}{\partial
r}\chi=-\Upsilon^{\prime}(r)\,\big{[}C\cos(n\theta)+D\sin(n\theta)\big{]}\,e^{-\frac{\eta}{\rho}m^{2}t}\,.$
(25b)

The velocity solutions (25) are not generally bounded in the limits $r\to{}0$,
$r\to\infty$ and $t\to\infty$. Auxiliary conditions must be imposed on the
coefficients in (25) to extract general solutions that are bounded in the
droplet domain $\Omega_{\textsc{d},0}$ (resp. ambient domain
$\mathbb{R}^{2}\setminus\overline{\Omega_{\textsc{d},0}}$) as $r\to{}0$ (resp.
$r\to\infty$) and as $t\to\infty$. To ensure boundedness of the solutions in
the limit $t\to\infty$, we insist that $\Re(m^{2})\geq{}0$. To assess the
boundedness of the solutions (25) in the spatial dependence, we note that the
Bessel function $J_{n}(mr)$ is singular at $r\to\infty$ for all
$m\in\mathbb{C}$ such that $|\Im(m)|>0$, and the Hankel function
$H^{(2)}_{n}(mr)$ is singular at the origin and at $r\to\infty$ for all
$m\in\mathbb{C}$ with $\Im(m)>0$. In addition, in relation to (7), we note
that $H^{(2)}_{n}(mr)$ vanishes in the limit $r\to\infty$ if $\Im(m)\leq{}0$.
Moreover, $r^{n}$ (resp. $r^{-n}$) is singular in the limit $r\to\infty$
(resp. $r\to{}0$). On account of their singularity at the origin,
$H^{(2)}_{n}$ and $r^{-n}$ are inadmissible in the droplet domain
$\Omega_{\textsc{d},0}\supset\\{0\\}$. Hence, in the droplet domain, it must
hold that $B=0$ and $F=0$. Conversely, $J_{n}(mr)$ and $r^{n}$ are
inadmissible in the ambient domain, in view of their singularity at
$r\to\infty$. Hence, in the ambient domain, it must hold that $A=0$ and $E=0$.

Summarizing, we obtain the following general bounded complex-valued velocity
solutions in the droplet and ambient domains:

$\displaystyle u_{\textsc{d},r}$
$\displaystyle=e^{-\frac{\eta_{\textsc{d}}}{\rho_{\textsc{d}}}m_{\textsc{d}}^{2}t}\,\Big{[}D_{\textsc{d}}\cos(n_{\textsc{d}}\theta)-C_{\textsc{d}}\sin(n_{\textsc{d}}\theta)\Big{]}$
$\displaystyle\hskip
100.0pt\times\frac{n_{\textsc{d}}}{r}\Big{[}AJ_{n_{\textsc{d}}}(m_{\textsc{d}}r)+Er^{n_{\textsc{d}}}\Big{]}\,,$
(26a) $\displaystyle u_{\textsc{d},\theta}$
$\displaystyle=e^{-\frac{\eta_{\textsc{d}}}{\rho_{\textsc{d}}}m_{\textsc{d}}^{2}t}\,\Big{[}C_{\textsc{d}}\cos(n_{\textsc{d}}\theta)+D_{\textsc{d}}\sin(n_{\textsc{d}}\theta)\Big{]}$
$\displaystyle\hskip
100.0pt\times\Big{[}-A\Big{(}m_{\textsc{d}}J_{n_{\textsc{d}}-1}(m_{\textsc{d}}r)-\frac{n_{\textsc{d}}}{r}J_{n_{\textsc{d}}}(m_{\textsc{d}}r)\Big{)}-En_{\textsc{d}}r^{n_{\textsc{d}}-1}\Big{]}\,,$
(26b) $\displaystyle u_{\textsc{a},r}$
$\displaystyle=e^{-\frac{\eta_{\textsc{a}}}{\rho_{\textsc{a}}}m_{\textsc{a}}^{2}t}\,\Big{[}D_{\textsc{a}}\cos(n_{\textsc{a}}\theta)-C_{\textsc{a}}\sin(n_{\textsc{a}}\theta)\Big{]}$
$\displaystyle\hskip
100.0pt\times\frac{n_{\textsc{a}}}{r}\Big{[}BH^{(2)}_{n_{\textsc{a}}}(m_{\textsc{a}}r)+Fr^{-n_{\textsc{a}}}\Big{]}\,,$
(26c) $\displaystyle u_{\textsc{a},\theta}$
$\displaystyle=e^{-\frac{\eta_{\textsc{a}}}{\rho_{\textsc{a}}}m_{\textsc{a}}^{2}t}\,\Big{[}C_{\textsc{a}}\cos(n_{\textsc{a}}\theta)+D_{\textsc{a}}\sin(n_{\textsc{a}}\theta)\Big{]}$
$\displaystyle\hskip
100.0pt\times\Big{[}-B\Big{(}m_{\textsc{a}}H^{(2)}_{n_{\textsc{a}}-1}(m_{\textsc{a}}r)-\frac{n_{\textsc{a}}}{r}H^{(2)}_{n_{\textsc{a}}}(m_{\textsc{a}}r)\Big{)}+Fn_{\textsc{a}}r^{-n_{\textsc{a}}-1}\Big{]}\,,$
(26d)

subject $\Re(m_{i}^{2})\geq{}0$ $(i\in\\{\textsc{a},\textsc{d}\\})$ and
$\Im(m_{\textsc{a}})<0$.

####  3.2.3 Pressure solutions

To facilitate the derivation of the infinitesimal pressure solutions
associated with (26), we first note that by virtue of (11a) and (11b), the
pressure solutions are harmonic functions. Considering functions that are
sinusoidal and periodic in the angular dependence, that conform to (26) in the
temporal dependence, and that are appropriately bounded, we find the following
general expression for $p_{\textsc{d}}$:

$p_{\textsc{d}}(r,\theta,t)=e^{-\frac{\eta_{\textsc{d}}}{\rho_{\textsc{d}}}m_{\textsc{d}}^{2}t}\,r^{\tilde{n}_{\textsc{d}}}\big{(}\tilde{D}_{\textsc{d}}\cos(\tilde{n}_{\textsc{d}}\theta)+\tilde{C}_{\textsc{d}}\sin(\tilde{n}_{\textsc{d}}\theta)\big{)}\,,$
(27)

with $\tilde{n}_{\textsc{d}}\in\mathbb{Z}_{\geq 0}$ and
$\tilde{C}_{\textsc{d}},\tilde{D}_{\textsc{d}}\in\mathbb{C}$ arbitrary
constants. By substituting (26) and (27) into (11a), we deduce that
compatibility between (27) and (26) imposes that the index
$\tilde{n}_{\textsc{d}}\in\mathbb{Z}_{\geq 0}$ and coefficients
$\tilde{C}_{\textsc{d}},\tilde{D}_{\textsc{d}}\in\mathbb{C}$ satisfy

$\tilde{D}_{\textsc{d}}=D_{\textsc{d}}E\eta_{\textsc{d}}m_{\textsc{d}}^{2},\qquad\tilde{C}_{\textsc{d}}=-C_{\textsc{d}}E\eta_{\textsc{d}}m_{\textsc{d}}^{2},\qquad\tilde{n}_{\textsc{d}}=n_{\textsc{d}}\,.$
(28)

The infinitesimal pressure solution in the ambient domain can be determined
similarly. In summary, we obtain:

$\displaystyle p_{\textsc{d}}$
$\displaystyle=E\eta_{\textsc{d}}m_{\textsc{d}}^{2}\,e^{-\frac{\eta_{\textsc{d}}}{\rho_{\textsc{d}}}m_{\textsc{d}}^{2}t}\,r^{n_{\textsc{d}}}\big{[}D_{\textsc{d}}\cos(n_{\textsc{d}}\theta)-C_{\textsc{d}}\sin(n_{\textsc{d}}\theta)\big{]}\,,$
(29a) $\displaystyle p_{\textsc{a}}$
$\displaystyle=-F\eta_{\textsc{a}}m_{\textsc{a}}^{2}\,e^{-\frac{\eta_{\textsc{a}}}{\rho_{\textsc{a}}}m_{\textsc{a}}^{2}t}r^{-n_{\textsc{a}}}\big{[}D_{\textsc{a}}\cos(n_{\textsc{a}}\theta)-C_{\textsc{a}}\sin(n_{\textsc{a}}\theta)\big{]}\,.$
(29b)

###### Remark 4

It is noteworthy that the Bessel and Hankel functions that appear in the
velocity solutions (26) are absent in the pressure solutions (29). This can be
rationalized by noting that these special functions originated directly from
the (pressure free) vorticity equation and thus satisfy the momentum equation
for a uniform pressure field.

###  3.3 Interface conditions

The general solutions for the pressure and velocity fields in the droplet and
ambient domains according to (26) and (29), involve twelve unknown
coefficients: $A$, $B$, $C_{\textsc{d}}$, $C_{\textsc{a}}$, $D_{\textsc{d}}$,
$D_{\textsc{a}}$, $E$, $F$, $m_{\textsc{d}}$, $m_{\textsc{a}}$,
$n_{\textsc{d}}$ and $n_{\textsc{a}}$. We next extract from the general
solutions the subspace that complies with the kinematic interface conditions
(11c) and (11d), and the dynamic condition (11e), by introducing auxiliary
conditions on the coefficients.

####  3.3.1 Kinematic compatibility conditions

The kinematic conditions (11c)-(11d) can be equivalently reformulated as:

$\displaystyle\boldsymbol{u}_{\textsc{d}}\cdot\boldsymbol{n}_{0}=\boldsymbol{u}_{\textsc{a}}\cdot\boldsymbol{n}_{0}=\mathcal{V}_{1}\qquad$
$\displaystyle\textrm{ on }\Gamma_{0}\,,$ (30a)
$\displaystyle\boldsymbol{u}_{\textsc{d}}\cdot\boldsymbol{t}_{0}=\boldsymbol{u}_{\textsc{a}}\cdot\boldsymbol{t}_{0}$
$\displaystyle\textrm{ on }\Gamma_{0}\,.$ (30b)

The generating solutions of the interface normal and tangent vectors
corresponding to the circular droplet, $\boldsymbol{n}_{0}$ and
$\boldsymbol{t}_{0}$, simply coincide with the radial and angular basis
vectors, respectively. The kinematic condition (30a) (resp. (30b)) thus
pertains to the radial (resp. angular) components of $\boldsymbol{u}$ in (26a)
and (26c) (resp. (26b) and (26d)). To impose (30a), we require the
infinitesimal interface velocity $\mathcal{V}_{1}$ corresponding to (9):

$\mathcal{V}(\theta,t)=\partial_{t}{}R_{\delta}(\theta,t)=-\delta
R_{0}\gamma\,e^{-\gamma
t}\,\big{(}\beta\cos(k\theta)+\sqrt{1-\beta^{2}}\sin(k\theta)\big{)}\,,$ (31)

where we retain the entire complex form, to facilitate the exposition. Noting
that $\mathcal{V}_{0}$ vanishes and, hence,
$\mathcal{V}=\delta\mathcal{V}_{1}$, we infer from (30a) that

$\displaystyle
e^{-\frac{\eta_{\textsc{d}}}{\rho_{\textsc{d}}}m_{\textsc{d}}^{2}t}\,\Big{[}D_{\textsc{d}}\cos(n_{\textsc{d}}\theta)-C_{\textsc{d}}\sin(n_{\textsc{d}}\theta)\Big{]}\,\Big{[}R^{-1}_{0}An_{\textsc{d}}J_{n_{\textsc{d}}}(m_{\textsc{d}}R_{0})+En_{\textsc{d}}R_{0}^{n_{\textsc{d}}-1}\Big{]}$
(32)
$\displaystyle\qquad=e^{-\frac{\eta_{\textsc{a}}}{\rho_{\textsc{a}}}m_{\textsc{a}}^{2}t}\,\Big{[}D_{\textsc{a}}\cos(n_{\textsc{a}}\theta)-C_{\textsc{a}}\sin(n_{\textsc{a}}\theta)\Big{]}\,\Big{[}R^{-1}_{0}Bn_{\textsc{a}}H^{(2)}_{n_{\textsc{a}}}(m_{\textsc{a}}R_{0})+Fn_{\textsc{a}}R_{0}^{-n_{\textsc{a}}-1}\Big{]}$
$\displaystyle\qquad\qquad=-R_{0}\gamma\,e^{-\gamma
t}\,\big{(}\beta\cos(k\theta)+\sqrt{1-\beta^{2}}\sin(k\theta)\big{)}\,.$

The equalities in (32) must hold for all $\theta\in[0,2\pi)$ and all
$t\in\mathbb{R}_{>0}$. Keeping $\theta$ fixed and varying $t$, one can infer
that the temporal exponents must coincide. Subsequently, by fixing $t$ and
varying $\theta$, it follows that all parameters that characterize the
trigonometric terms must be the same. Hence,

$\displaystyle\frac{\eta_{\textsc{d}}}{\rho_{\textsc{d}}}m_{\textsc{d}}^{2}=\frac{\eta_{\textsc{a}}}{\rho_{\textsc{a}}}m_{\textsc{a}}^{2}=\gamma\,,$
(33a) $\displaystyle n_{\textsc{d}}=n_{\textsc{a}}=k\,,$ (33b) $\displaystyle
D_{\textsc{d}}=D_{\textsc{a}}=\beta\,,$ (33c) $\displaystyle
C_{\textsc{d}}=C_{\textsc{a}}=-\sqrt{1-\beta^{2}}\,.$ (33d)

In the sequel, we continue to use $m_{\textsc{d}}$ and $m_{\textsc{a}}$ in the
arguments of the Bessel and Hankel functions, but we tacitly suppose the
relation to $\gamma$ per (33a). Substitution of (33) in (30) yields the
following three conditions on $A,B,E$ and $F$:

$\displaystyle AkR^{-1}_{0}J_{k}(m_{\textsc{d}}R_{0})+EkR_{0}^{k-1}$
$\displaystyle=-\gamma R_{0}\,,$ (34a) $\displaystyle
BkR_{0}^{-1}H^{(2)}_{k}(m_{\textsc{a}}R_{0})+FkR_{0}^{-k-1}$
$\displaystyle=-\gamma R_{0}\,,$ (34b)
$\displaystyle-A\left(m_{\textsc{d}}J_{k-1}(m_{\textsc{d}}R_{0})-kR_{0}^{-1}J_{k}(m_{\textsc{d}}R_{0})\right)-EkR_{0}^{k-1}\qquad\qquad$
$\displaystyle+B\left(m_{\textsc{a}}H_{k-1}^{(2)}(m_{\textsc{a}}R_{0})-kR_{0}^{-1}H_{k}^{(2)}(m_{\textsc{a}}R_{0})\right)-FkR_{0}^{-k-1}$
$\displaystyle=0\,.$ (34c)

####  3.3.2 Dynamic compatibility conditions

The dynamic condition (11e) can be separated into radial and angular
components to form the following two conditions:

$\displaystyle-\eta_{\textsc{d}}\boldsymbol{n}_{0}\left(\nabla\boldsymbol{u}_{\textsc{d}}+\left(\nabla\boldsymbol{u}_{\textsc{d}}\right)^{T}\right)\boldsymbol{n}_{0}+p_{\textsc{d}}+\eta_{\textsc{a}}\boldsymbol{n}_{0}\left(\nabla\boldsymbol{u}_{\textsc{a}}+\left(\nabla\boldsymbol{u}_{\textsc{a}}\right)^{T}\right)\boldsymbol{n}_{0}-p_{\textsc{a}}$
$\displaystyle=\sigma_{\textsc{d}\textsc{a}}\kappa_{1}$ $\displaystyle\textrm{
on }\Gamma\,,$ (35a)
$\displaystyle-\eta_{\textsc{d}}\boldsymbol{t}_{0}\left(\nabla\boldsymbol{u}_{\textsc{d}}+\left(\nabla\boldsymbol{u}_{\textsc{d}}\right)^{T}\right)\boldsymbol{n}_{0}+\eta_{\textsc{a}}\boldsymbol{t}_{0}\left(\nabla\boldsymbol{u}_{\textsc{a}}+\left(\nabla\boldsymbol{u}_{\textsc{a}}\right)^{T}\right)\boldsymbol{n}_{0}$
$\displaystyle=0$ $\displaystyle\textrm{ on }\Gamma\,.$ (35b)

To elaborate on the dynamic interface condition (35a), we require the first-
order perturbation of the interface curvature. From the postulated interface
displacement (9), the complex-valued form of the curvature can be derived up
to second-order terms:

$\kappa(\theta,t)=\kappa_{0}+\delta\kappa_{1}(\theta,t)+\mathcal{O}(\delta^{2})\\\
=R_{0}^{-1}+\delta R^{-1}_{0}e^{-\gamma
t}\big{(}k^{2}-1\big{)}\big{(}\beta\cos(k\theta)-\sqrt{1-\beta^{2}}\sin(k\theta)\big{)}+\mathcal{O}(\delta^{2})\,.$
(36)

From the expressions for the velocity (26) and pressure (29), the relation
between the coefficients in (33), and the dynamic conditions (35), it then
follows that

$\displaystyle-A\eta_{\textsc{d}}\Big{[}2m_{\textsc{d}}k{R_{0}}^{-1}J_{k-1}(m_{\textsc{d}}{R_{0}})-2(k+1)k{R_{0}}^{-2}J_{k}(m_{\textsc{d}}{R_{0}})\Big{]}$
$\displaystyle\qquad+B\eta_{\textsc{a}}\Big{[}2m_{\textsc{a}}k{R_{0}}^{-1}H_{k-1}^{(2)}(m_{\textsc{a}}{R_{0}})-2(k+1)k{R_{0}}^{-2}H_{k}^{(2)}(m_{\textsc{a}}{R_{0}})\Big{]}$
$\displaystyle\qquad+E\eta_{\textsc{d}}\Big{[}m_{\textsc{d}}^{2}{R_{0}}^{k}-2(k-1)k{R_{0}}^{k-2}\Big{]}+F\eta_{\textsc{a}}\Big{[}m_{\textsc{a}}^{2}{R_{0}}^{-k}-2(k+1)k{R_{0}}^{-k-2}\Big{]}$
$\displaystyle\qquad\qquad=\sigma_{\textsc{d}\textsc{a}}{R_{0}}^{-1}(k^{2}-1)\,,$
(37a)
$\displaystyle-A\eta_{\textsc{d}}\Big{[}2m_{\textsc{d}}{R_{0}}^{-1}J_{k-1}(m_{\textsc{d}}{R_{0}})+\big{(}m_{\textsc{d}}^{2}-2(k+1)k{R_{0}}^{-2}\big{)}J_{k}(m_{\textsc{d}}{R_{0}})\Big{]}$
$\displaystyle\qquad+B\eta_{\textsc{a}}\Big{[}2m_{\textsc{a}}{R_{0}}^{-1}H_{k-1}^{(2)}(m_{\textsc{a}}{R_{0}})+\left(m_{\textsc{a}}^{2}-2(k+1)k{R_{0}}^{-2}\right)H_{k}^{(2)}(m_{\textsc{a}}{R_{0}})\Big{]}$
$\displaystyle\qquad+2E\eta_{\textsc{d}}(k-1)k{R_{0}}^{k-2}+2F\eta_{\textsc{a}}(k+1)k{R_{0}}^{-k-2}$
$\displaystyle\qquad\qquad=0\,.$ (37b)

###  3.4 Dispersion relation

For each wave number $k\in\mathbb{N}_{\geq{}2}$, the corresponding
characteristic temporal response coefficient of the solution,
$\gamma\in\mathbb{C}$, as well as the mode shapes, encoded in the remaining
free parameters, follow from a solution-existence condition. To elucidate this
condition, we first recall that the general bounded complex-valued velocity
and pressure solutions of the partial-differential equations (11a)-(11b),
subject to the limit condition (7), are given by (26) and (29). These general
solutions contain twelve coefficients. Eight of these coefficients are
determined by the kinematic interface condition (11c), in accordance with
(33). The kinematic conditions (11c) and (11d) imply that the remaining four
coefficients, $A,B,E$ and $F$ must satisfy the three identities in (34). The
dynamic condition (11e) demands that, in addition, these four coefficients
satisfy the two identities in (37). The remaining five conditions on the
coefficients can be cast in the form

$\underbrace{\begin{pmatrix}{a}_{11}&0&{a}_{13}&0\\\ 0&{a}_{22}&0&{a}_{24}\\\
{a}_{31}&{a}_{32}&{a}_{33}&{a}_{34}\\\ {a}_{41}&{a}_{42}&{a}_{43}&{a}_{44}\\\
{a}_{51}&{a}_{52}&{a}_{53}&{a}_{54}\end{pmatrix}}_{\boldsymbol{A}(k,\gamma)}\begin{pmatrix}A\\\
B\\\ E\\\ F\end{pmatrix}=\underbrace{\begin{pmatrix}-\gamma R_{0}\\\ -\gamma
R_{0}\\\ 0\\\ \sigma_{\textsc{d}\textsc{a}}R^{-1}_{0}(k^{2}-1)\\\
0\end{pmatrix}}_{{\boldsymbol{b}}(k,\gamma)}\,,$ (38)

in such a manner that the first three equations in (38) represent (34) and the
latter two represent (37). Noting the dependence of (34) and (37) on the wave
number $k$, and recalling the dependence of $m_{\textsc{d}}$ and
$m_{\textsc{a}}$ in these equation on the temporal response coefficient
$\gamma$ via (33a), we infer that the entries of ${\boldsymbol{A}}$ depend on
$k$ and $\gamma$. With five constraints and four unknowns, the system of
equations (38) is formally over-constrained, and a solution is non-existent
unless the right-hand-side vector ${\boldsymbol{b}}(k,\gamma)$ is in the
column space of ${\boldsymbol{A}}(k,\gamma)$. The relation between the
existence of a solution and the condition

${\boldsymbol{b}}(k,\gamma)\in\operatorname{span}(\operatorname{col}({\boldsymbol{A}}(k,\gamma)))\,,$
(39)

is indicative of the fact that only specific combinations of the wave number
$k\in\mathbb{N}_{\geq{}2}$ and the temporal response coefficient
$\gamma\in\mathbb{C}$ in the postulated interface configuration (8) correspond
to a natural response of the droplet.

To determine the combinations $(k,\gamma)$ for which the existence condition
(39) is fulfilled, we note that (39) is equivalent to

$\det\big{(}(\boldsymbol{A}\mid\boldsymbol{b})(k,\gamma)\big{)}=0\,,$ (40)

where $(\boldsymbol{A}\mid\boldsymbol{b})$ corresponds to $\boldsymbol{A}$
augmented by $\boldsymbol{b}$. The equivalence between (39) and (40) follows
from the fact that the column vectors of $\boldsymbol{A}$ are linearly
independent for all $(k,\gamma)$ and, hence, the augmented matrix is singular
if and only if the vector $\boldsymbol{b}$ resides in the column space of
$\boldsymbol{A}$. Moreover, by virtue of the linear independence of the
columns of $\boldsymbol{A}$, if (40) holds, then (38) has a unique solution.
This solution corresponds to the coefficients $(A,B,E,F)$ that, in combination
with (33), define the droplet and ambient velocity-pressure pairs
corresponding to $(k,\gamma)$ according to (26) and (29).

To facilitate and generalize the root-finding of the determinant in (40), we
non-dimensionalize the matrix entries of the augmented matrix based on the
droplet density, $\rho_{\textsc{d}}$, droplet viscosity, $\eta_{\textsc{d}}$
and droplet radius $R_{0}$. The non-dimensionalized parameters are indicated
with a tilde diacritic. Additionally, we introduce the following condensed
notation:

$\displaystyle\mathcal{J}=J_{k}(\tilde{m}_{\textsc{d}}(\tilde{\gamma}))\,,\qquad$
$\displaystyle\mathcal{H}=H^{(2)}_{k}(\tilde{m}_{\textsc{a}}(\tilde{\gamma}))\,,\qquad$
$\displaystyle\zeta=2(k-1)k\,,$
$\displaystyle\hat{\mathcal{J}}=\tilde{m}_{\textsc{d}}(\tilde{\gamma})J_{k-1}(\tilde{m}_{\textsc{d}}(\tilde{\gamma}))\,,\qquad$
$\displaystyle\hat{\mathcal{H}}=\tilde{m}_{\textsc{a}}(\tilde{\gamma})H^{(2)}_{k-1}(\tilde{m}_{\textsc{a}}(\tilde{\gamma}))\,,\qquad$
$\displaystyle\xi=2(k+1)k\,.$

The non-dimensionalized augmented matrix can then be expressed as

$(\tilde{\boldsymbol{A}}|\tilde{\boldsymbol{b}})(k,\tilde{\gamma})=\\\
\begin{pmatrix}k\mathcal{J}&0&k&0&-\tilde{\gamma}\\\
0&k\mathcal{H}&0&k&-\tilde{\gamma}\\\
-\hat{\mathcal{J}}+k\mathcal{J}&\hat{\mathcal{H}}-k\mathcal{H}&-k&-k&0\\\
-2k\hat{\mathcal{J}}+\xi\mathcal{J}&\tilde{\eta}_{\textsc{a}}\big{[}2k\hat{\mathcal{H}}-\xi\mathcal{H}\big{]}&\tilde{m}^{2}_{\textsc{d}}-\zeta&\tilde{\eta}_{\textsc{a}}\big{[}\tilde{m}^{2}_{\textsc{a}}-\xi\big{]}&\tilde{\sigma}_{\textsc{d}\textsc{a}}\left(k^{2}-1\right)\\\
-2\hat{\mathcal{J}}-\left(\tilde{m}^{2}_{\textsc{d}}-\xi\right)\mathcal{J}&\tilde{\eta}_{\textsc{a}}\big{[}2\hat{\mathcal{H}}+\left(\tilde{m}^{2}_{\textsc{a}}-\xi\right)\mathcal{H}\big{]}&\zeta&-\tilde{\eta}_{\textsc{a}}\xi&0\end{pmatrix}.$
(41)

It is not generally feasible to determine the roots of
$\det((\tilde{\boldsymbol{A}}|\tilde{\boldsymbol{b}})(k,\tilde{\gamma}))$ with
respect to $\tilde{\gamma}$ in closed form and, in practice, it is necessary
to revert to a numerical root-finding algorithm. Once a root has been
determined, one can extract the kernel of the augmented matrix (41) and scale
the corresponding vector such that its fifth entry is minus one, to obtain the
coefficients $\tilde{A},\tilde{B},\tilde{E},\tilde{F}$.

###### Remark 5

The roots of
$\det((\tilde{\boldsymbol{A}}|\tilde{\boldsymbol{b}})(k,\tilde{\gamma}))$ are
not unique: one can infer that

$\big{[}(\tilde{\boldsymbol{A}}|\tilde{\boldsymbol{b}})(k,\tilde{\gamma}{}^{*})\big{]}=\big{[}(\tilde{\boldsymbol{A}}|\tilde{\boldsymbol{b}})(k,\tilde{\gamma})\big{]}^{*}\,,$
(42)

where $(\cdot)^{*}$ denotes complex conjugation. Because the eigenvalues of
the complex conjugate of a matrix are the complex conjugates of the original
eigenvalues, it follows that if $\tilde{\gamma}$ is a root of
$\det((\tilde{\boldsymbol{A}}|\tilde{\boldsymbol{b}})(k,\tilde{\gamma}))$,
then so is $\tilde{\gamma}{}^{*}$. Since, in addition, it must hold that
$\Re(\tilde{\gamma})>0$, it suffices to consider roots in the fourth quadrant
of the complex plane. Noting that the entries of the augmented matrix are
analytic functions, one can infer that so is its determinant. This implies
that the roots of
$\det((\tilde{\boldsymbol{A}}|\tilde{\boldsymbol{b}})(k,\tilde{\gamma}))$ form
a totally disconnected set and, accordingly, for each root there exists a
neighborhood in which that root is unique. A detailed investigation of the
uniqueness of the roots of
$\det((\tilde{\boldsymbol{A}}|\tilde{\boldsymbol{b}})(k,\tilde{\gamma}))$ in
the fourth quadrant is beyond the scope of this work. In our numerical root-
finding procedure, we have verified that there are no other roots in a region
around the found root.

By virtue of the complex representation of the interface parametrization (9),
we obtain the real-valued velocity and pressure fields by taking the real
parts of (26) and (29) after substitution of the relations (33), and the
temporal response coefficient $\gamma$ and the corresponding coefficients
$A,B,E$, and $F$. Table 1 provides computed parameter values for the physical
setting outlined in Table 2, representing a water-in-air picoliter-sized
droplet. For completeness, we mention that we have applied _Mathematica_ ’s
root-finder to determine $\tilde{\gamma}:=\tilde{\gamma}_{k}$ in the fourth
quadrant of the complex plane such that
$\det((\tilde{\boldsymbol{A}}|\tilde{\boldsymbol{b}})(k,\tilde{\gamma}_{k}))=0$.
The dimensions of parameters $E$ and $F$ depend on the mode number $k$. As a
result, as $k$ increases, the values of $E$ and $F$ grow rapidly, conveying
that these parameters are ill-conditioned in terms of $k$. Furthermore, as the
mode number $k$ increases, the frequency and damping rate of the corresponding
oscillation, both encoded in $\gamma$, increase. This implies that if a
droplet sustains an initial perturbation that is characterized by multiple
modes, the higher wave-number modes decay quickly, and low order modes
dominate the long-term dynamics of viscous-in-viscous oscillating droplets.

Table 1: Modal solution parameter values for a water droplet of radius $R=\sqrt{2}\\!\times\\!10^{1}\,\textrm{\textmu m}$ suspended in air. $k$ | $\gamma\,\,[\textrm{ s}^{-1}]$ | $A\,\,[10^{-11}\textrm{ m}^{2}\textrm{s}^{-1}]$ | $B\,\,[10^{-4}\textrm{ m}^{2}\textrm{s}^{-1}]$ | $E\,\,[10^{5k-7}\textrm{ m}^{2-k}\textrm{s}^{-1}]$ | $F\,\,[10^{-5k-4}\textrm{ m}^{2+k}\textrm{s}^{-1}]$
---|---|---|---|---|---
2 | $18788.18393$ | $4664.160935$ | $1.142101737$ | $2.618366811$ | $1.904600443$
| $-390396.1271\,i$ | $+137.5420287\,i$ | $+4.485396518\,i$ | $+194.0864949\,i$ | $+1.883632636\,i$
3 | $53722.95262$ | $376.1242266$ | $7.390256000$ | $3.392996002$ | $4.168582183$
| $-777097.6733\,i$ | $-953.0487542\,i$ | $+6.296111161\,i$ | $+182.1390754\,i$ | $+3.618449771\,i$
4 | $104333.5005$ | $-123.9009766$ | $14.90136892$ | $3.548390076$ | $7.942104534$
| $-1223261.977\,i$ | $-195.0741019\,i$ | $+1.040690443\,i$ | $+152.0972759\,i$ | $+6.209497220\,i$
5 | $170139.1932$ | $-53.99114858$ | $16.35189970$ | $3.366783037$ | $14.10551086$
| $-1722774.521\,i$ | $-11.92064586\,i$ | $-10.42851063\,i$ | $+121.2100680\,i$ | $+10.15645338\,i$
6 | $250754.7829$ | $-12.54654878$ | $8.074025870$ | $3.017014114$ | $23.99270873$
| $-2270070.007\,i$ | $+6.299980519\,i$ | $-20.94948730\,i$ | $+94.14089637\,i$ | $+16.16915562\,i$

##  4 Numerical experiments

The free-boundary problem (5) formally represents the sharp-interface limit of
the Abels–Garcke–Grün NSCH model (1), provided that the mobility is
appropriately scaled in the limit $\varepsilon\to{}+0$. For sufficiently small
$\delta$, the oscillating-droplet solutions derived in Section 3 can therefore
serve to investigate the approach of the diffuse-interface solution to the
sharp-interface limit solution. In this section, we investigate this sharp-
interface limit numerically, by means of an adaptive finite-element method.
Specifically, we focus on the scaling of the mobility parameter
$m:=m_{\varepsilon}$ in the limit $\varepsilon\to{}+0$, and investigate the
deviation of the diffuse-interface solution from the sharp-interface solution
in relation to $m$. As reference solutions, we consider the lowest mode of
oscillation ($k=2$), as well as the next higher doubly symmetric mode ($k=4$),
for a viscous droplet in a viscous ambient with parameter values according to
Table 2. The corresponding coefficients of the velocity solution (26) and
pressure solution (29) fields are presented in Table 1.

Table 2: Physical and numerical parameter values of the considered numerical experiments. Entries marked with the symbol $\ast$ indicate a range of values, which will be specified in the text. Droplet | Ambient | Interface | Numerical approximation
---|---|---|---
$\rho_{\textsc{d}}$ | $\eta_{\textsc{d}}$ | $\rho_{\textsc{a}}$ | $\eta_{\textsc{a}}$ | $\sigma_{\textsc{da}}$ | $\varepsilon$ | $m$ | $\tau$ | $h_{0}$ | $L_{\textrm{max}}$ | $K$
$\frac{\textrm{kg}}{\textrm{m}^{d}}$ | $\frac{\textrm{kg}\,\textrm{m}^{2-d}}{\textrm{s}}$ | $\frac{\textrm{kg}}{\textrm{m}^{d}}$ | $\frac{\textrm{kg}\,\textrm{m}^{2-d}}{\textrm{s}}$ | $\frac{\textrm{kg}\,\textrm{m}^{3-d}}{\textrm{s}^{2}}$ | m | $\frac{\textrm{m}^{d}\,\textrm{s}}{\textrm{kg}}$ |  $10$ µs | µm | — | —
$10^{3}$ | $10^{-3}$ | $1$ | $1.813\\!\times\\!10^{-5}$ | $7.28\\!\times\\!10^{-2}$ | $\ast$ | $\ast$ | $2^{-7}$ | $5$ | $\ast$ | $\ast$

###  4.1 Setup and discretization

The oscillating-droplet test cases that we consider pertain to doubly
symmetric modes. The setup of the test cases is similar to that in [11,
Sec.5]. To reduce computational expense, we exploit the symmetry of the
configurations and consider only one quarter of the droplet-ambient domain. We
regard a domain $\Omega=(0,50)^{2}\,\mu\text{m}^{2}$ and prescribe symmetry
conditions on
$\Gamma_{\text{sym}}:=\\{(x_{1},x_{2})\in\partial\Omega:\\{x_{1}=0\\}\cup\\{x_{2}=0\\}\\}$.
Because the linear sharp-interface solution is in fact defined on the
generating circular droplet domain and the corresponding ambient domain
according to (6a), while the diffuse-interface model exhibits a moving
interface, we prescribe auxiliary conditions in accordance with an initially
circular droplet. Specifically, with reference to (9), we select $t_{0}$ such
that $-\Im(\gamma)t_{0}=\pi/2$ and, hence, $R_{\delta}(\theta,t_{0})=R_{0}$,
and prescribe initial data corresponding to the reference solution at $t_{0}$
and boundary data corresponding to $t+t_{0}$. The complementary part of the
boundary, $\Gamma_{\text{ext}}:=\partial\Omega\setminus\Gamma_{\text{sym}}$,
is furnished with Dirichlet conditions for velocity and homogeneous Neumann
conditions for the order parameter and the chemical potential:

$\left.\begin{aligned}
\boldsymbol{u}(\cdot,t)&=\delta\boldsymbol{u}_{\textsc{a}}(\cdot,t_{0}+t)\\\
\partial_{n}\varphi&=0\\\
\partial_{n}\mu&=0\end{aligned}\right\\}\quad\text{on
}\Gamma_{\text{ext}},\text{ for }t\in[0,T)\,,$ (43)

where $\delta\boldsymbol{u}_{\textsc{a}}$ corresponds to the ambient velocity
solution (26) with appropriate coefficients and scaling $\delta$, and $[0,T)$
denotes the time-interval under consideration. For the aforementioned
combination of boundary conditions, the pressure variable $p$ is only
determined up to a constant. We impose the auxiliary condition that $p$
vanishes on average.

We impose an initial condition for the order parameter corresponding to a
circular interface, in accordance with the initial configuration of the sharp-
interface reference solution, viz.

$\varphi(\boldsymbol{x},0)=\varphi_{0}(\boldsymbol{x}):=\tanh\bigg{(}\frac{d_{\pm}(\boldsymbol{x},\Gamma_{0})}{\sqrt{2}\varepsilon}\bigg{)}\,,$
(44)

where $d_{\pm}(\boldsymbol{x},\Gamma_{0})$ represents the signed distance from
$\boldsymbol{x}$ to $\Gamma_{0}$. The function
$s\mapsto\tanh(s/\sqrt{2}\varepsilon)$ corresponds to an equilibrium solution
of the Cahn–Hilliard equations for the phase field in one spatial dimension
and, accordingly, the phase field (44) is meta-stable if $\varepsilon$ is
sufficiently small compared to the radius of curvature of $\Gamma_{0}$. In
conjunction with (44), we impose the following initial condition for velocity:

$\boldsymbol{u}(\boldsymbol{x},0)=\begin{cases}\delta\boldsymbol{u}_{\textsc{d}}(\boldsymbol{x},t_{0})&\quad\text{if
}\boldsymbol{x}\in\Omega_{\textsc{d},0}\,,\\\
\delta\boldsymbol{u}_{\textsc{a}}(\boldsymbol{x},t_{0})&\quad\text{if
}\boldsymbol{x}\in\Omega_{\textsc{a},0}\,,\end{cases}$ (45)

where the data in the right member of (45) corresponds to the the velocity
solutions according to (26) in the droplet and ambient domains. We select the
perturbation magnitude $\delta=10^{-2}$, after verifying that this choice
renders the linearization error negligible in comparison to the deviation
between the diffuse-interface and the sharp-interface solutions, for the
$\varepsilon$ considered below. Hence, the selected value of $\delta$ is
suitable for our investigation of the sharp-interface limit. The
characteristic parameters pertaining to the droplet and ambient fluids, and to
the interface are reported in Table 2.

To perform the numerical simulations, we make use of the adaptive finite-
element approximation method presented in [11]. For coherence, we present a
concise overview of the numerical methodology. The weak form of the NSCH
equations (1) is discretized with respect to the spatial dependence with
$P3-P2$ (Taylor-Hood) $\mathcal{C}^{0}$ truncated hierarchical B-splines (see
[17, 10, 13]) for the velocity and pressure fields, and $P3$ $\mathcal{C}^{0}$
truncated hierarchical B-splines for the order parameter and chemical
potential; see [31, §3.1] for further details. The adaptive-refinement
procedure is guided by a two-level hierarchical a-posteriori error estimate,
and follows the standard SEMR (Solve $\rightarrow$ Estimate $\rightarrow$ Mark
$\rightarrow$ Refine) process [8, 12]. To improve the robustness of the
solution procedure on the coarse meshes that occur in the sequence of adaptive
refinements within each time step, an $\varepsilon$-continuation process is
introduced, in which the thickness parameter $\varepsilon$ (and, in
conjunction, the mobility $m\propto\varepsilon^{3}$) is enlarged for the first
$K$ iterations of the adaptive refinement process; see [11, 31] for details.
In each time step, the fluid domain is initially covered with a uniform mesh
comprising $10\times 10$ elements, corresponding to an initial mesh width
$h_{0}=5\,\mu\text{m}$, and we perform $L_{\textrm{max}}$ refinement steps.
Refinement steps $L=0,1,\ldots,K-1$ make use of the $\varepsilon$\- and
$m$-continuation process, while in refinement steps
$L=K,\ldots,L_{\textrm{max}}$ the original parameter values for $\varepsilon$
and $m$ are used. A skew-symmetric formulation according to [23] is used for
the convective term in the Navier–Stokes equations, enhancing the stability of
the discrete approximation by eliminating potential artificial energy
production due to deviations from solenoidality in pure-species regions. On
the coarse meshes, a first order Backward Euler scheme with second order
contractive-expansive splitting of the double-well potential with
stabilization [33] is employed. On the finest mesh, a second order
Crank–Nicolson scheme is applied with implicit treatment of the double-well
potential. The second order Crank–Nicolson scheme provides significant better
accuracy than the Backward Euler scheme; cf. e.g. [20]. For the temporal
discretization, we employ a time-step size $\tau=2^{-7}\times 10\,\mu$s. The
parameter setting of the numerical procedure is also summarized in Table 2.
The nonlinear algebraic systems corresponding to the discretized NSCH
equations, are solved with a Newton procedure, in which the linear tangent
problems are solved with GMRES with a preconditioner based on a partition of
the NSCH system into NS and CH subsystems; see [11] for further details.

To illustrate the setup of the numerical experiments, and the resemblance
between the analytic sharp-interface solution and the numerical approximation
of the diffuse-interface solution for sufficiently small $\varepsilon$, we
conduct numerical experiments with interface-thickness parameter
$\varepsilon=2^{-10}\times{}10^{2}\,\text{\textmu m}$, mobility
$m=9.5272\times{}10^{-13}\,\text{m}^{2}\text{s}/\text{kg}$, maximum number of
refinement levels $L_{\textrm{max}}=7$, and number of continuation levels
$K=5$. Figures 2 and 3 display snapshots of the velocity field and pressure
field at six time instants. The top (resp. bottom) half of each panel displays
the velocity (resp. pressure) field. The right (resp. left) half of each panel
depicts the diffuse-interface simulation (resp. sharp-interface solution). The
figures convey that the sharp-interface solutions and the diffuse-interface
solutions are visually indistinguishable.

Figure 2: Snapshots of the magnitude of the velocity field $|\boldsymbol{u}|$
(top) and pressure field $p$ (bottom) throughout half a droplet oscillation of
mode $k=2$. The left half of each panel displays the analytical sharp-
interface solution, while the right half displays the numerical diffuse-
interface solution. Figure 3: Snapshots of the magnitude of the velocity field
$|\boldsymbol{u}|$ (top) and pressure field $p$ (bottom) throughout half a
droplet oscillation of mode $k=4$. The left half of each panel displays the
analytical sharp-interface solution, while the right half displays the
numerical diffuse-interface solution.

###  4.2 Optimal mobility scaling

To elucidate the dependence of the diffuse-interface solution in the sharp-
interface limit $\varepsilon\to+0$ on the scaling of the mobility
$m:=m_{\varepsilon}$, we conduct numerical experiments for a range of
combinations of $\varepsilon$ and $m$. For each combination of $\varepsilon$
and $m$, we determine the deviation relative to the sharp-interface solution
according to:

$\operatorname{dev}(\varepsilon,m)=\frac{{\left|\kern-1.07639pt\left|\kern-1.07639pt\left|\boldsymbol{u}^{\varepsilon,m}-\boldsymbol{u}\right|\kern-1.07639pt\right|\kern-1.07639pt\right|}_{\Omega\times(0,T)}}{{\left|\kern-1.07639pt\left|\kern-1.07639pt\left|\boldsymbol{u}\right|\kern-1.07639pt\right|\kern-1.07639pt\right|}_{\Omega\times(0,T)}}\quad\text{with}\quad{\left|\kern-1.07639pt\left|\kern-1.07639pt\left|\boldsymbol{u}\right|\kern-1.07639pt\right|\kern-1.07639pt\right|}_{\Omega\times(0,T)}=\frac{1}{T}\int\limits_{0}^{T}|\\!|\boldsymbol{u}(\cdot,t)|\\!|_{L^{2}(\Omega)}\,\text{d}t\,,$
(46)

where the considered length of the time interval, $T$, corresponds to half a
period of oscillation. We regard a set of decreasing interface thickness
parameters
$\varepsilon\in\mathscr{E}:=\\{2^{0},\ldots,2^{-3}\\}\varepsilon_{\text{max}}$
relative to the baseline interface thickness
$\varepsilon_{\text{max}}=2^{-7}\\!\times\\!10^{2}\,\textrm{\textmu
m}=7.8125\\!\times\\!10^{-1}\,\textrm{\textmu m}$. The baseline interface
thickness parameter corresponds to approximately 5% of the droplet radius. For
each $\varepsilon$, we consider mobility parameters in (a relevant subset of)
the set $m\in\mathscr{M}:=\\{2^{0},2^{-1},\ldots,2^{-12}\\}m_{\text{max}}$
with $m_{\text{max}}=2.4389632\times
10^{-10}\,\textrm{m}^{d}\textrm{s}/\textrm{kg}$. The range of mobility
parameters has been determined empirically such that
$[2^{-12},1]m_{\text{max}}$ includes the optimal mobility, i.e. the one for
which $\operatorname{dev}(\varepsilon,m)$ is minimal, for all
$\varepsilon\in\mathscr{E}$. It is to be noted that
$\mathcal{E}\times\mathcal{M}$ contains various monomial scalings of the
mobility with respect to the interface thickness, viz.
$m\propto\varepsilon^{l}$ with $l\in\\{0,1,2,3\\}$.

Tables 3 and 4 present the deviations $\operatorname{dev}(\varepsilon,m)$ for
the two modes of oscillation, $k=2$ and $k=4$, respectively. For each
$\varepsilon\in\mathcal{E}$, the entry corresponding to the mobility
$m\in\mathcal{M}$ that yields the smallest deviation, is highlighted. One can
observe that, indeed, the mobility corresponding to the minimal deviation
decreases as $\varepsilon$ decreases. More precisely, for both modes, the
optimal scaling of the mobility parameter with the interface-thickness
parameter appears to lie between $m\propto{}\varepsilon$ and
$m\propto{}\varepsilon^{2}$. One may moreover note that the entries
corresponding to the optimal mobility decrease by a factor of approximately
two if $\varepsilon$ is halved, which indicates that for the considered
droplet-oscillation case, the diffuse-interface solution approaches the sharp-
interface solution at rate $\mathcal{O}(\varepsilon)$, provided that the
mobility in the diffuse-interface model is appropriately scaled.

Table 3: Deviation between the diffuse-interface solution and the sharp-interface solution according to (46), for $k=2$ for one half period of oscillation. $m/m_{\textrm{max}}$ |  $2^{-3}$ $\varepsilon_{\textrm{max}}$ |  $2^{-2}$ $\varepsilon_{\textrm{max}}$ |  $2^{-1}$ $\varepsilon_{\textrm{max}}$ | $\varepsilon_{\textrm{max}}$
---|---|---|---|---
$2^{-12}$ | $6.8295\\!\times\\!10^{-2}$ |  |  |
$2^{-11}$ |  |  |  |
$2^{-10}$ |  |  |  |
$2^{-9\phantom{0}}$ | $1.2634\\!\times\\!10^{-2}$ | $7.3709\\!\times\\!10^{-2}$ |  |
$2^{-8\phantom{0}}$ | $8.3534\\!\times\\!10^{-3}$ | $4.4566\\!\times\\!10^{-2}$ |  |
$2^{-7\phantom{0}}$ | $1.2941\\!\times\\!10^{-2}$ | $2.1360\\!\times\\!10^{-2}$ | $1.0190\\!\times\\!10^{-1}$ | $1.3056\\!\times\\!10^{-1}$
$2^{-6\phantom{0}}$ | $2.3777\\!\times\\!10^{-2}$ | $1.8713\\!\times\\!10^{-2}$ | $6.6942\\!\times\\!10^{-2}$ | $1.2485\\!\times\\!10^{-1}$
$2^{-5\phantom{0}}$ |  | $4.2245\\!\times\\!10^{-2}$ | $3.5126\\!\times\\!10^{-2}$ | $1.0979\\!\times\\!10^{-1}$
$2^{-4\phantom{0}}$ |  |  | $6.2628\\!\times\\!10^{-2}$ | $7.9845\\!\times\\!10^{-2}$
$2^{-3\phantom{0}}$ | $1.5565\\!\times\\!10^{-1}$ | $1.6940\\!\times\\!10^{-1}$ | $1.5228\\!\times\\!10^{-1}$ | $6.1487\\!\times\\!10^{-2}$
$2^{-2\phantom{0}}$ |  |  |  | $1.7640\\!\times\\!10^{-1}$
$2^{-1\phantom{0}}$ |  |  |  | $4.2783\\!\times\\!10^{-1}$
$1\phantom{{}^{0-1}}$ |  |  |  | $8.2210\\!\times\\!10^{-1}$
Table 4: Deviation between the diffuse-interface solution and the sharp-interface solution according to (46), for $k=4$ for one half period of oscillation. $m/m_{\textrm{max}}$ |  $2^{-3}$ $\varepsilon_{\text{max}}$ |  $2^{-2}$ $\varepsilon_{\text{max}}$ |  $2^{-1}$ $\varepsilon_{\text{max}}$ | $\varepsilon_{\text{max}}$
---|---|---|---|---
$2^{-11}$ | $1.2922\\!\times\\!10^{-1}$ |  |  |
$2^{-10}$ |  |  |  |
$2^{-9\phantom{0}}$ |  |  |  |
$2^{-8\phantom{0}}$ | $2.7543\\!\times\\!10^{-2}$ | $1.4419\\!\times\\!10^{-1}$ |  |
$2^{-7\phantom{0}}$ | $2.3628\\!\times\\!10^{-2}$ | $9.3099\\!\times\\!10^{-2}$ |  |
$2^{-6\phantom{0}}$ | $3.2089\\!\times\\!10^{-2}$ | $5.5083\\!\times\\!10^{-2}$ |  |
$2^{-5\phantom{0}}$ | $5.2591\\!\times\\!10^{-2}$ | $5.2646\\!\times\\!10^{-2}$ | $1.5428\\!\times\\!10^{-1}$ |
$2^{-4\phantom{0}}$ |  | $9.4735\\!\times\\!10^{-2}$ | $1.0710\\!\times\\!10^{-1}$ | $2.9576\\!\times\\!10^{-1}$
$2^{-3\phantom{0}}$ |  |  | $1.2054\\!\times\\!10^{-1}$ | $2.7358\\!\times\\!10^{-1}$
$2^{-2\phantom{0}}$ | $3.3162\\!\times\\!10^{-1}$ | $3.5898\\!\times\\!10^{-1}$ | $2.7331\\!\times\\!10^{-1}$ | $2.0927\\!\times\\!10^{-1}$
$2^{-1\phantom{0}}$ |  |  |  | $2.5501\\!\times\\!10^{-1}$
$1\phantom{{}^{0-1}}$ |  |  |  | $6.2529\\!\times\\!10^{-1}$

To provide a more precise assessment of the optimal scaling relation
$m:=m_{\varepsilon}$, we determine for each $\varepsilon$ the optimal value of
$m$ based on a quadratic log-log interpolation around the minimal values in
Tables 3 and 4. Figure 4 plots the optimal value of $m$ versus $\varepsilon$.
For both wave numbers, we observe an optimal scaling
$m\propto\varepsilon^{a_{\textrm{opt}}}$ with $a_{\textrm{opt}}\approx 1.7$.
It is noteworthy that the constant of proportionality in the scaling relation
is different for the two modes, and that the graphs are offset in the
$\varepsilon$-dependence by a factor of approximately two, i.e. the optimal
mobility for $k=4$ is approximately $2^{a_{\textrm{opt}}}$ larger than the
optimal mobility for $k=2$. This suggests that the optimal mobility in fact
scales with $(\varepsilon/\ell)^{a_{\textrm{opt}}}$, where $\ell$ represents
another characteristic length scale of the interface which, for the considered
droplet-oscillation test case, is proportional to the wave length of the
perturbation. This observation calls for further investigation, but we
consider a detailed analysis of this aspect beyond the scope of the present
work.

Figure 4: Optimal mobility $m$ obtained from quadratic interpolation around
the minima in Tables 3 and 4.

###  4.3 Sensitivity to the proportionality constant

In the previous section, we established optimal values of the mobility
parameter and inferred an optimal scaling
$m\propto\varepsilon^{a_{\textrm{opt}}}$ in the sharp-interface limit. The
results in Section 4.2 also convey that the constant of proportionality in the
scaling relation $m_{\varepsilon}=\mathscr{C}\varepsilon^{a_{\textrm{opt}}}$
depends on the configuration and dynamics of the interface. This raises the
question how sensitive the solution is to suboptimality of the proportionality
constant in the scaling relation.

To elucidate the sensitivity of the deviation of the diffuse-interface
solution to the sharp-interface solution with respect to the mobility in the
limit $\varepsilon\to{}+0$, Figure 5 (resp. Figure 6) plots for each
$\varepsilon\in\mathcal{E}$ the ratio of the deviation
$\operatorname{dev}(\varepsilon,m)$ in the columns of Table 3 (resp. Table 4)
to the minimal deviation
$\operatorname{dev}(\varepsilon,m_{\textrm{opt},\varepsilon})$ versus the
ratio $m/m_{\textrm{opt},\varepsilon}$. Noting that the curves in Figures 5
and 6 exhibit a vanishing slope near $m/m_{\textrm{opt},\varepsilon}=1$, one
can conclude that in the vicinity of the optimal mobility, the relative
deviation is essentially independent of the mobility. However, for larger
departures from the optimal mobility and sufficiently small $\varepsilon$, the
relative deviation increases almost linearly in
$\max(m/m_{\textrm{opt},\varepsilon},m_{\textrm{opt},\varepsilon}/m)$. For the
largest $\varepsilon\in\mathcal{E}$, the relative deviation appears to be less
sensitive to underestimation than to overestimation of the mobility. However,
the results plotted with solid markers in Figures 5 and 6 indicate that in the
sharp-interface limit, the relative deviation is equally sensitive to under-
and overestimation of the mobility.

Figure 5: Normalized deviation
$\operatorname{dev}(\varepsilon,m)/\operatorname{dev}(\varepsilon,m_{\textrm{opt},\varepsilon})$
versus normalized mobility $m/m_{\textrm{opt},\varepsilon}$ for $k=2$. Figure
6: Normalized deviation
$\operatorname{dev}(\varepsilon,m)/\operatorname{dev}(\varepsilon,m_{\textrm{opt},\varepsilon})$
versus normalized mobility $m/m_{\textrm{opt},\varepsilon}$ for $k=4$.

###  4.4 Suboptimal mobility scaling and convergence in the sharp-interface
limit

Figure 7: Error convergence rates for optimal and suboptimal scalings of
mobility $m$, for mode of oscillation $k=2$. Figure 8: Error convergence rates
for optimal and suboptimal scalings of mobility $m$, for mode of oscillation
$k=4$.

To illustrate the effect of the scaling of the mobility on the approach to the
sharp-interface limit solution, Figures 7 and 8 plot the deviation
$\operatorname{dev}(\varepsilon,(\varepsilon/\varepsilon_{\textrm{max}})^{a}\,m_{0})$
versus $\varepsilon$ for $a\in\\{0,\ldots,3\\}$. Herein, $m_{0}$ corresponds
to the optimal sampled mobility for $\varepsilon_{\textrm{max}}$, viz.
$m_{0}=2^{-3}m_{\textrm{max}}$ for $k=2$ and $m_{0}=2^{-2}m_{\textrm{max}}$
for $k=4$; cf. Tables 3 and 4. For reference, the figures also contain the
estimated minimal deviation obtained by minimization of the quadratic
interpolation, $\operatorname{dev}(\varepsilon,m_{\textrm{opt},\varepsilon})$.
The figures convey that for the optimal scaling of the mobility, the diffuse-
interface solution approaches the sharp-interface solution essentially at
order $\varepsilon$, i.e.
$\operatorname{dev}(\varepsilon,m_{\textrm{opt},\varepsilon})=\mathcal{O}(\varepsilon)$
as $\varepsilon\to{}+0$. For the linear and quadratic scaling of the mobility
with the interface thickness, $m\propto\varepsilon^{a}$ with $a\in\\{1,2\\}$,
i.e. the integer scalings of the mobility closest to the optimum, we observe
convergence to the sharp-interface solution, but at a suboptimal (sublinear)
rate. For the constant and cubic scalings of the mobility,
$m\propto\varepsilon^{a}$ with $a\in\\{0,3\\}$, the deviation
$\operatorname{dev}(\varepsilon,m_{0}(\varepsilon/\varepsilon_{\textrm{max}})^{a})$
does not vanish as $\varepsilon\to{}+0$, i.e. the diffuse-interface solution
does not convergence to the sharp-interface solution. For
$m\propto\varepsilon^{0}$, this confirms the known result that for constant
mobility, the Abels–Garcke–Grün NSCH model converges to the nonclassical
sharp-interface Navier–Stokes/Mullins–Sekerka model; see [3].

###### Remark 6

It is noteworthy that the subcubic scaling of the mobility,
$m\propto\varepsilon^{a}$ with $0<{}a<3$, which is necessary to converge to
the classical sharp-interface solution in the limit $\varepsilon\to{}+0$,
implies that the characteristic diffusive time scale
$T_{\mathrm{diff}}:=\varepsilon^{3}/\sigma{}m$ associated with the diffuse
interface, approaches zero in the sharp-interface limit. Consequently, for any
$\varepsilon$-independent characteristic time scale $T_{*}$ in the problem
under consideration, e.g. the period of oscillation of a droplet, it holds
that the ratio $T_{\mathrm{diff}}/T_{*}\to{}+0$ as $\varepsilon\to{}+0$. For
numerical time-integration methods for the NSCH equations, it is therefore
essential that such methods are robust in the limit
$T_{\mathrm{diff}}/\tau\to+0$ (with $\tau$ denoting the time-step size), to
avoid excessive computational complexity in the sharp-interface limit.

##  5 Conclusions

Diffuse-interface binary-fluid models bear significant potential for
describing complex phenomena in fluid mechanics, such as topological changes
of the fluid-fluid interface and dynamic wetting, by virtue of their implicit
representation of the interface. In the absence of topological changes of the
interface, diffuse-interface models should reduce to corresponding classical
sharp-interface models in the so-called sharp-interface limit, viz. if the
interface-thickness parameter, $\varepsilon$, passes to zero. Contemporary
understanding of the sharp-interface limit is however incomplete and, in
particular, the scaling of the mobility parameter, $m$, with $\varepsilon$ as
$\varepsilon\to{}+0$ is incompletely understood. In this article, we
investigated the limit behavior of the Abels–Garcke–Grün Navier-Stokes–Cahn-
Hilliard model for the classical case of an oscillating droplet in two
dimensions, by means of an adaptive finite-element methodology.

To provide reference sharp-interface solutions, we derived new two-dimensional
analytical expressions for the velocity and pressure fields for small-
amplitude oscillations of a viscous droplet in a viscous ambient fluid with
different densities and viscosities.

For mode numbers $k=2,4$, we compared the solutions of the Navier-Stokes–Cahn-
Hilliard model to the corresponding analytical solutions for a decreasing
sequence of interface-thickness parameters and a suitably chosen sequence of
mobility parameters. Based on an analysis of the deviation between the
diffuse-interface solution and the sharp-interface solution, we deduced that
$m_{\textrm{opt},\varepsilon}\propto\varepsilon^{a_{\textrm{opt}}}$ with
${a_{\textrm{opt}}}\approx 1.7$ corresponds to the optimal scaling of the
mobility in the sharp-interface limit. We found that this optimal scaling is
universal for $k=2$ and $k=4$. However, we also observed that the factor of
proportionality differs by a factor of approximately $2^{a_{\textrm{opt}}}$,
indicating that the optimal mobility in fact scales with the ratio of
$\varepsilon$ to another characteristic length scale, in this case,
proportional to the wave length of the perturbation. The observed dependence
of the optimal mobility on the configuration and motion of the interface
warrants further investigation.

For the optimal scaling of the mobility parameter, we observed that the
deviation between the diffuse-interface solution and its sharp-interface limit
decreases according to $\mathcal{O}(\varepsilon)$ in the limit
$\varepsilon\to{}+0$. Our investigation of suboptimal integer scalings of the
mobility conveyed that the sharp-interface limit is also attained for a linear
and quadratic scaling of the mobility, but not for a constant or cubic
scaling. For the linear and quadratic scaling, the approach of the diffuse-
interface solution to the sharp-interface solution occurs at a suboptimal
(sublinear) rate. The fact that the scaling of the mobility with $\varepsilon$
must be subcubic, implies that the characteristic diffusive time scale
$\varepsilon^{3}/\sigma{}m$ (with $\sigma$ as surface tension) passes to zero
in the sharp-interface limit.

## Acknowledgments

This research was partly conducted within the Industrial Partnership Program
Fundamental Fluid Dynamics Challenges in Inkjet Printing (FIP), a joint
research program of Canon Production Printing, Eindhoven University of
Technology, University of Twente, and the Netherlands Organization for
Scientific Research (NWO). T.H.B. Demont and S.K.F. Stoter gratefully
acknowledge financial support through the FIP program. All simulations have
been performed using the open source software package Nutils (www.nutils.org)
[32].

The authors thank Dr. S.W. Rienstra of the Department of Mathematics and
Computer Science at Eindhoven University of Technology for insightful
discussions with respect to the interpretation of matched-asymptotic-expansion
results for sharp-interface limits of the NSCH equations in the literature.

## References

  * [1] H. Abels, D. Depner, and H. Garcke.  Existence of weak solutions for a diffuse interface model for two-phase flows of incompressible fluids with different densities.  J. Math. Fluid Mech., 15(3):453–480, 2013.
  * [2] H. Abels, D. Depner, and H. Garcke.  On an incompressible Navier–Stokes/Cahn–Hilliard system with degenerate mobility.  Ann. I. H. Poincaré – AN, 30(6):1175–1190, 2013.
  * [3] H. Abels and H. Garcke.  Weak Solutions and Diffuse Interface Models for Incompressible Two-Phase Flows, pages 1267–1327.  Springer International Publishing, Cham, 2018.
  * [4] H. Abels, H. Garcke, and G. Grün.  Thermodynamically consistent, frame indifferent diffuse interface models for incompressible two-phase flows with different densities.  Math. Mod. Meth. Appl. Sci., 22:1150013, 2012.
  * [5] H. Abels and D. Lengeler.  On sharp interface limits for diffuse interface models for two-phase flows.  Interfaces Free Boundaries, 395–418, 2014.
  * [6] S. Arrhenius.  Über die innere reibung verdünnter wässeriger lösungen.  Z. Physik. Chem., 1U(1):285–298, 1887.
  * [7] J. Barrett, J. Blowey, and H. Garcke.  Finite element approximation of the Cahn–Hilliard equation with degenerate mobility.  SIAM J. Numer. Anal., 37(1):286–318, 2016/11/21 1999.
  * [8] S. Bertoluzza, R.H. Nochetto, A. Quarteroni, K.G. Siebert, and A. Veeser.  Multiscale and Adaptivity: Modeling, Numerics and Applications.  CIME, 2012.
  * [9] H. Bonart, C. Kahle, and J.-U. Repke.  Comparison of energy stable simulation of moving contact line problems using a thermodynamically consistent Cahn–Hilliard Navier–Stokes model.  J. Comput. Phys., 399:108959, 2019.
  * [10] J.A. Cottrell, T.J.R. Hughes, and Y. Bazilevs.  Isogeometric Analysis: Toward Integration of CAD and FEA.  Wiley, Chichester, 2009.
  * [11] T.H.B. Demont, G.J. van Zwieten, C. Diddens, and E.H. van Brummelen.  A robust and accurate adaptive approximation method for a diffuse-interface model of binary-fluid flows.  Comput. Methods Appl. Mech. Engrg., 400:115563, 2022.
  * [12] W. Dörfler.  A convergent adaptive algorithm for Poisson’s equation.  SIAM J. Numer. Anal., 33:1106–1124, 1996.
  * [13] C. Giannelli, B. Jüttler, and H. Speleers.  THB-splines: The truncated basis for hierarchical splines.  Comput. Aided Geom. Des., 29:485–498, 2012.
  * [14] A. Giorgini.  Well-posedness of the two-dimensional Abels–Garcke–Grün model for two-phase flows with unmatched densities.  Calc. Var. Partial Differ. Equ., 60(3):100, 2021.
  * [15] G. Grün, F. Guillén-González, and S. Metzger.  On fully decoupled, convergent schemes for diffuse interface models for two-phase flow with general mass densities.  Commun. Comput. Phys., 19(5):1473–1502, 2016.
  * [16] P.C. Hohenberg and B.I. Halperin.  Theory of dynamic critical phenomena.  Rev. Mod. Phys., 49:435–479, 1977.
  * [17] T.J.R. Hughes, J.A. Cottrell, and Y. Bazilevs.  Isogeometric analysis: CAD, finite elements, NURBS, exact geometry and mesh refinement.  Comput. Methods Appl. Mech. Engrg., 194:4135–4195, 2005.
  * [18] D. Jacqmin.  Calculation of two-phase Navier-Stokes flows using phase-field modeling.  J. Comput. Phys., 155:96–127, 10 1999.
  * [19] D. Jacqmin.  Contact-line dynamics of a diffuse fluid interface.  J. Fluid Mech., 402:57–88, 2000.
  * [20] V. John, G. Matthies, and J. Rang.  A comparison of time-discretization/linearization approaches for the incompressible Navier–Stokes equations.  Comput. Methods Appl. Mech. Eng., 195(44):5995–6010, 2006.
  * [21] V.V. Khatavkar, P.D. Anderson, and H.E.H. Meijer.  On scaling of diffuse–interface models.  Chem. Eng. Sci., 61(8):2364–2378, 2006.
  * [22] H. Lamb.  Hydrodynamics.  The University Press, New York, 6th edition, 1932.
  * [23] W. Layton.  Introduction to the Numerical Analysis of Incompressible Flows.  The Society for Industrial and Applied Mathematics, Philadelphia, PA, USA, 2008.
  * [24] J. Lowengrub and L. Truskinovsky.  Quasi-incompressible Cahn-Hilliard fluids and topological transitions.  Proceedings of the Royal Society of London A: Mathematical, Physical and Engineering Sciences, 454:2617–2654, 1998.
  * [25] F. Magaletti, F. Picano, M. Chinappi, L. Marino, and C.M. Casciola.  The sharp-interface limit of the Cahn–Hilliard/Navier–Stokes model for binary fluids.  J. Fluid Mech., 714:95–126, 2013.
  * [26] C.A. Miller and L.E. Scriven.  The oscillations of a fluid droplet immersed in another fluid.  J. Fluid Mech., 32(3):417–435, 1968.
  * [27] W.H. Reid.  The oscillations of a viscous liquid drop.  Q. Appl. Math., 18(1):86–89, 1960.
  * [28] P. Seppecher.  Moving contact lines in the Cahn-Hilliard theory.  Int. J. Engng. Sci., 34:977–992, 7 1996.
  * [29] M. Shokrpour Roudbari, G. Şimşek, E.H. van Brummelen, and K.G. van der Zee.  Diffuse-interface two-phase flow models with different densities: A new quasi-incompressible form and a linear energy-stable method.  Math. Models Methods Appl. Sci., 28:733–770, 2018.
  * [30] J.W. Strutt.  VI. on the capillary phenomena of jets.  Proceedings of the Royal Society of London, 29(196-199):71–97, 2021/11/18 1879.
  * [31] E.H. van Brummelen, T.H.B. Demont, and G.J. van Zwieten.  An adaptive isogeometric analysis approach to elasto-capillary fluid-solid interaction.  Int. J. Numer. Meth. Engng., 122(19):5331–5352, 2021/09/30 2021\.
  * [32] G.J. van Zwieten, J. van Zwieten, C.V. Verhoosel, E. Fonn, T.M. van Opstal, and W. Hoitinga.  Nutils (ver. 5.0).  2019\.
  * [33] 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.  Int. J. Numer. Method. Biomed. Eng., 30(2):180–203, 2014.
  * [34] P. Yue and J.J. Feng.  Wall energy relaxation in the Cahn–Hilliard model for moving contact lines.  Phys. Fluids, 23:012106, 2011.
  * [35] P. Yue, C. Zhou, and J.J. Feng.  Sharp-interface limit of the Cahn–Hilliard model for moving contact lines.  J. Fluid Mech., 645:279–294, 2010.