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 $\textrm {D}$ and fluid $\textrm {A}$. We consider an open time interval $(0, t_{fin}) \subseteq \mathbb{R}_{>0}$ and a spatial domain corresponding to a simply connected time-independent subset $\varOmega \subseteq \mathbb{R}^d$ ($d=2,3$). We make use of a NSCH type model that describes the evolution of a so-called order parameter $\varphi \in [-1,1]$ representing pure species $\textrm {D}$ and $\textrm {A}$ when $\varphi =1$ and $\varphi =-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 (Abels et al. Reference Abels, Garcke and Grün2012)

(2.1a, 2.1b, 2.1c, and 2.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 $\varPsi$ are given as

(2.2a)$$\begin{gather} {} \boldsymbol{J} := m \frac{ \rho_{\textrm{A}} - \rho_{\textrm{D}} } { 2 } \boldsymbol{\nabla} \mu , \end{gather}$$

(2.2b)$$\begin{gather}\boldsymbol{\tau} := \eta ( \boldsymbol{\nabla} \boldsymbol{u} + ( \boldsymbol{\nabla} \boldsymbol{u} )^{\rm T} ) , \end{gather}$$

(2.2c)$$\begin{gather}\boldsymbol{\zeta} := - \sigma \varepsilon \boldsymbol{\nabla} \varphi \otimes \boldsymbol{\nabla} \varphi + \boldsymbol{I} \left( \frac{ \sigma \varepsilon } { 2 } | \boldsymbol{\nabla} \varphi |^2 + \frac{ \sigma } { \varepsilon } \varPsi \right) , \end{gather}$$

(2.2d)$$\begin{gather}\varPsi \left( \varphi \right) := \tfrac{1}{4} (\varphi^2 -1)^2 . \end{gather}$$

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 _\textrm {DA}$ according to $2\sqrt {2}\sigma =3\sigma _\textrm {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]$ (Grün, Guillén-González & Metzger Reference Grün, Guillén-González and Metzger2016). We include a density extension that ensures positive densities even for the non-physical scenario $\varphi \notin [-1,1]$ (Bonart, Kahle & Repke Reference Bonart, Kahle and Repke2019)

(2.3) \begin{equation} \rho(\varphi) = \begin{cases} \frac{ 1 } { 4 } \rho_{\textrm{A}}, & \varphi \leqslant - 1 - 2 \lambda ,\\ \frac{ 1 } { 4 } \rho_{\textrm{A}} + \frac{ 1 } { 4 } \rho_{\textrm{A}} \lambda^{{-}2} ( 1 + 2 \lambda + \varphi )^2, & \varphi \in ( - 1 - 2 \lambda , - 1 - \lambda) ,\\ \dfrac{ 1 + \varphi } { 2 } \rho_{\textrm{D}} + \dfrac{ 1 - \varphi } { 2 } \rho_{\textrm{A}}, & \varphi \in [ - 1 - \lambda, 1 + \lambda ] ,\\ \rho_{\textrm{D}} + \frac{ 3 } { 4 } \rho_{\textrm{A}} - \frac{ 1 } { 4 } \rho_{\textrm{A}} \lambda^{{-}2} ( 1 + 2 \lambda - \varphi )^2, & \varphi \in ( 1 + \lambda, 1 + 2 \lambda ) ,\\ \rho_{\textrm{D}} + \frac{ 3 } { 4 } \rho_{\textrm{A}}, & \varphi \geqslant 1 + 2 \lambda , \end{cases} \end{equation}

where $\lambda = \rho _{\textrm {A}} / ( \rho _{\textrm {D}} - \rho _{\textrm {A}} )$. For the viscosity interpolation, we apply the Arrhenius mixture-viscosity model (Arrhenius Reference Arrhenius1887)

(2.4)\begin{equation} \log \eta( \varphi ) = \frac{ \left( 1 + \varphi \right) \varLambda \log \eta_{\textrm{D}} + \left( 1 - \varphi \right) \log \eta_{\textrm{A}} } { \left( 1 + \varphi \right) \varLambda + \left( 1 - \varphi \right) } \,, \end{equation}

where $\varLambda = { \rho _{\textrm {D}} M_{\textrm {A}} }/{\rho _{\textrm {A}} M_{\textrm {D}} }$ is the intrinsic volume ratio between the two fluids (with $M_{\textrm {A}}$ and $M_{\textrm {D}}$ their respective molar masses).

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:

(2.5a)$$\begin{gather} \rho_{i}\partial_t \boldsymbol{u}_i + \rho_{i} \left( \boldsymbol{u}_i \boldsymbol{\cdot} \boldsymbol{\nabla} \right) \boldsymbol{u}_i - \eta_{i} \Delta \boldsymbol{u}_i + \boldsymbol{\nabla} p_i = 0\quad \textrm{in}\ \varOmega_{i} = \varOmega_{i}(t) , \end{gather}$$

(2.5b)$$\begin{gather}\boldsymbol{\nabla}\boldsymbol{\cdot} \boldsymbol{u}_i = 0 \quad \textrm{in } \varOmega_{i} , \end{gather}$$

(2.5c)$$\begin{gather}\boldsymbol{u}_{i}\boldsymbol{\cdot}\boldsymbol{n} = \mathcal{V}\quad \textrm{on } \varGamma = \varGamma(t) , \end{gather}$$

(2.5d)$$\begin{gather} {[}\![ \boldsymbol{u}\boldsymbol{\cdot} \boldsymbol{t}_j ]\!] = 0 \quad \textrm{on } \varGamma \text{, for }j=1,\ldots,d-1 , \end{gather}$$

(2.5e)$$\begin{gather} {[}\![ -( \boldsymbol{\nabla} \boldsymbol{u} + ( \boldsymbol{\nabla} \boldsymbol{u} )^{\rm T} ) \boldsymbol{n} + p \boldsymbol{n} ]\!] = \sigma_{\textrm{D}\textrm{A}} \kappa \boldsymbol{n} \quad \textrm{on } \varGamma , \end{gather}$$

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

As opposed to the NSCH model, the sharp-interface model (2.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 $\varOmega _i(t)$ and $\varGamma (t)$. There is an intrinsic coupling between the velocity field and the evolution of $\varOmega _i$ and $\varGamma$, according to (2.5c). Because this equation holds on both sides of the interface and $\mathcal {V}$ is single valued, (2.5c) implies $[\![ \boldsymbol {u} \boldsymbol {\cdot } \boldsymbol {n} ]\!] = 0$.

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}_{\geqslant 0}$ and $\theta \in [0,2{\rm \pi} )$ and the coordinate transformation $\boldsymbol {x}=(x_1,x_2)=r(\cos \theta,\sin \theta )$. We regard perturbations of the interface $\varGamma _0=\partial \varOmega _{\textrm {D},0}$ corresponding to the following parametrization:

(3.3)\begin{equation} \varGamma_{\delta}(t)=\{\boldsymbol{x}\in\mathbb{R}^2: \boldsymbol{x}=R_{\delta}(\theta,t)\,(\cos(\theta),\sin(\theta)),\theta\in[0,2{\rm \pi})\} , \end{equation}

where

(3.4) \begin{align} R_{\delta} ( \theta , t ) & = R_0 + R_0\, \delta( \beta \cos ( k \theta ) + \sqrt{1-\beta^2 } \, \sin ( k \theta ) ) \cos(\nu t) {\rm e}^{- \alpha t} \nonumber\\ & = {\rm Re} ( R_0 + R_0\, \delta ( \beta \cos ( k \theta ) + \sqrt{1-\beta^2 } \, \sin ( k \theta ) ) {\rm e}^{- \gamma t} ) . \end{align}

The interface configuration (3.3)–(3.4) 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\leqslant \beta \leqslant 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}(\varOmega _{\textrm {D}})=\operatorname {meas}(\varOmega _{\textrm {D},0})+{O}(\delta ^2)$ and the barycentre of $\varOmega _{\textrm {D}}$ is located at the origin. This implies that $k \in \mathbb {N}_{\geqslant 2}$. The damping rate $\alpha \geqslant 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 (3.4) provides a representation of $R_{\delta }$ as the real part of a complex-valued function, with $\gamma := \alpha - {\rm i} \nu$. This form enables us to condense some of the expressions that appear in the sequel.

In conjunction with the interface configuration (3.3)–(3.4), we consider linear asymptotic solutions of the sharp-interface problem of the form

(3.5)\begin{equation} \left(\boldsymbol{u}_{i},p_i,\boldsymbol{n}, \boldsymbol{t},\kappa,\mathcal{V}\right)= \left(\boldsymbol{u}_{i},p_i,\boldsymbol{n}, \boldsymbol{t},\kappa,\mathcal{V}\right)_0 +\delta\left(\boldsymbol{u}_{i},p_i,\boldsymbol{n}, \boldsymbol{t},\kappa,\mathcal{V}\right)_1, \end{equation}

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

(3.6a)$$\begin{gather} \rho\partial_t \boldsymbol{u}_{i,1} - \eta \Delta \boldsymbol{u}_{i,1} + \boldsymbol{\nabla} p_{i,1} = 0 \quad \textrm{in } \varOmega_{i,0} , \end{gather}$$

(3.6b)$$\begin{gather}\boldsymbol{\nabla}\boldsymbol{\cdot} \boldsymbol{u}_{i,1} = 0 \quad \textrm{in } \varOmega_{i,0} , \end{gather}$$

(3.6c)$$\begin{gather}\boldsymbol{u}_{i,1} \boldsymbol{\cdot} \boldsymbol{n}_0 = \mathcal{V}_1 \quad \textrm{on } \varGamma_0 , \end{gather}$$

(3.6d)$$\begin{gather} {[}\![ \boldsymbol{u}_{1} \boldsymbol{\cdot} \boldsymbol{t}_0 ]\!] = 0 \quad \textrm{on } \varGamma_0 , \end{gather}$$

(3.6e)$$\begin{gather} {[}\![ - \eta \left(\boldsymbol{\nabla} \boldsymbol{u}_{1} + \left(\boldsymbol{\nabla} \boldsymbol{u}_{1} \right)^{\rm T} \right) \boldsymbol{n}_0 + p_{1} \boldsymbol{n}_0 ]\!] = \sigma_{\textrm{D}\textrm{A}} \kappa_1 \boldsymbol{n}_0 \quad \textrm{on } \varGamma_0 , \end{gather}$$

for $i \in \{\textrm {D},\textrm {A}\}$. The nonlinear advective term has dropped since the only nonlinear 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 (3.6), because its multiplication with $\boldsymbol {u}_{i,1}$ is a second-order term, its dot product with $\boldsymbol {u}_{i,0}$ vanishes and $[\![ p_0 \boldsymbol {n}_1 ]\!] = \sigma _{\textrm {D}\textrm {A}} R^{-1}_0 \boldsymbol {n}_1$ cancels with the right-hand side $\sigma _{\textrm {D}\textrm {A}} \kappa _0 \boldsymbol {n}_1$.

The remainder of this section is devoted to finding general solutions to (3.6) for different perturbation wavenumbers $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 (3.6a)–(3.6b) for both the droplet $\textrm {D}$ and ambient $\textrm {A}$ domains. We proceed by deriving the vorticity equation corresponding to (3.6a), 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 (3.6a) and by using the identity $\boldsymbol {\nabla } \times \boldsymbol {\nabla }({\cdot }) = 0$. Introducing the vorticity $\omega =\boldsymbol {\nabla } \times \boldsymbol {u}$, we infer from (3.6a) that

(3.7)\begin{equation} \rho \partial_t \left( \boldsymbol{\nabla} \times \boldsymbol{u} \right) - \eta \varDelta \left( \boldsymbol{\nabla} \times \boldsymbol{u} \right) + \boldsymbol{\nabla} \times\boldsymbol{\nabla}p= \rho \partial_t \omega - \eta \Delta \omega = 0. \end{equation}

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 (3.7), we assume the following separation of variables form:

(3.8)\begin{equation} \omega ( r,\theta , t ) := P(r) \varTheta(\theta) T(t). \end{equation}

Substitution in the polar coordinate representation of the diffusion equation gives

(3.9)\begin{equation} P(r) \varTheta(\theta) T'(t) = \frac{\eta}{\rho r^2}P(r) \varTheta''(\theta) T(t)+ \frac{\eta}{\rho r}P'(r) \varTheta(\theta) T(t) + \frac{\eta}{\rho}P''(r) \varTheta(\theta) T(t) , \end{equation}

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

(3.10a)$$\begin{gather} T'(t) ={-}\frac{\eta}{ \rho} m^2 T(t) \quad m \in \mathbb{C}, \end{gather}$$

(3.10b)$$\begin{gather}\varTheta''(\theta) ={-}n^2 \varTheta(\theta) \quad n \in \mathbb{C} , \end{gather}$$

(3.10c)$$\begin{gather}r^2 P''(r) + r P'(r) + (m^2 r^2 - n^2) P(r) = 0 , \end{gather}$$

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

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

(3.11a)$$\begin{gather} T (t) = c\, {\rm e}^{- ({\eta}/{ \rho}) m^2 t}, \end{gather}$$

(3.11b)$$\begin{gather}\varTheta (\theta) = c_1\, {\rm e}^{ i n t} + c_2 {\rm e}^{ {-}i n t}. \end{gather}$$

From the periodicity of the droplet perturbations in the angular dependence, conforming to (3.3), we infer that $n\in \mathbb {Z}_{\geqslant {}0}$. The arbitrary constants $c_1$ and $c_2$ can then be selected such that (3.11b) reduces to the sum of two real-valued trigonometric functions

(3.12)\begin{equation} \varTheta (\theta) = C\cos( n \theta) + D\sin( n \theta)\quad (n \in \mathbb{Z}_{\geqslant0}) , \end{equation}

where $C,D \in \mathbb {R}$ are coefficients that determine the angular orientation of the solution. Regarding (3.10c), we note that this equation corresponds to Bessel's equation with a complex-valued scaling $m\in \mathbb {C}$. Solutions of (3.10c) 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 (3.10c) consists of a linear combination of two Bessel functions of order $n\in \mathbb {Z}_{\geqslant 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, ${\rm J}_n$, and a Hankel function of the second kind, $H^{(2)}_n$,

(3.13)\begin{equation} P(r) = A m^2 {\rm J}_n ( m r ) + B m^2 H^{(2)}_n ( m r ) , \end{equation}

with $A,B \in \mathbb {C}$. Substitution of (3.11a), (3.12) and (3.13) into (3.8) gives the general rotationally periodic solution of the vorticity equation (3.7)

(3.14)\begin{equation} \omega(r,\theta,t) = [A m^2 {\rm J}_n ( m r ) + B m^2 H^{(2)}_n ( m r )] \left[ C\cos( n \theta) + D\sin( n \theta) \right] {\rm e}^{- ({\eta}/{ \rho}) m^2 t} , \end{equation}

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 (3.14), we introduce a streamfunction $\chi$ according to

(3.15)\begin{equation} \Delta \chi ={-} \omega . \end{equation}

The velocity can be retrieved from this streamfunction as

(3.16)\begin{equation} \boldsymbol{u} = \frac{1}{r} \frac{\partial}{\partial\theta} \chi\, \boldsymbol{e}_r - \frac{\partial}{\partial r} \chi \, \boldsymbol{e}_\theta. \end{equation}

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

(3.17)\begin{equation} \chi (r,\theta,t) = \varUpsilon ( r ) \left[ C\cos( n \theta) + D\sin( n \theta) \right] {\rm e}^{- ({\eta}/{ \rho}) m^2 t} . \end{equation}

Substitution of (3.17) into (3.15) leads to the following ordinary differential equation for $\varUpsilon$:

(3.18)\begin{equation} \varUpsilon'' ( r ) + \frac{1}{r} \varUpsilon' ( r ) - \frac{n^2}{r^2} \varUpsilon ( r ) ={-} A m^2 {\rm J}_n ( m r ) - B m^2 H^{(2)}_n ( m r ) . \end{equation}

The general solution to the non-homogeneous ordinary differential equation (3.18) is given by

(3.19)\begin{equation} \varUpsilon(r) = E r^{n} + F r^{{-}n} + A {\rm J}_n ( m r) + B H^{(2)}_{n} ( m r ) . \end{equation}

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

(3.20a)$$\begin{gather} u_r = \frac{1}{r} \frac{\partial}{\partial\theta} \chi = \frac{n}{r} \varUpsilon ( r ) \left[ D \cos( n \theta) - C \sin( n \theta) \right] {\rm e}^{- ({\eta}/{\rho}) m^2 t} \,, \end{gather}$$

(3.20b)$$\begin{gather}u_\theta ={-} \frac{\partial}{\partial r} \chi ={-} \varUpsilon' ( r ) \left[ C \cos( n \theta) + D \sin( n \theta) \right] {\rm e}^{- ({\eta}/{ \rho}) m^2 t} . \end{gather}$$

The velocity solutions (3.20) 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 (3.20) to extract general solutions that are bounded in the droplet domain $\varOmega _{\textrm {D},0}$ (respectively ambient domain $\mathbb {R}^2\setminus \overline {\varOmega _{\textrm {D},0}}$) as $r\to {}0$ (respectively $r\to \infty$) and as $t\to \infty$. To ensure boundedness of the solutions in the limit $t\to \infty$, we insist that ${\rm Re}(m^2)\geqslant {}0$. To assess the boundedness of the solutions (3.20) in the spatial dependence, we note that the Bessel function ${\rm J}_n(mr)$ is singular at $r\to \infty$ for all $m\in \mathbb {C}$ such that $|{\rm 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 ${\rm Im}(m)>0$. In addition, in relation to (3.2), we note that $H^{(2)}_{n}(mr)$ vanishes in the limit $r\to \infty$ if ${\rm Im}(m)\leqslant {}0$. Moreover, $r^n$ (respectively $r^{-n}$) is singular in the limit $r\to \infty$ (respectively $r\to {}0$). On account of their singularity at the origin, $H^{(2)}_{n}$ and $r^{-n}$ are inadmissible in the droplet domain $\varOmega _{\textrm {D},0}\supset \{0\}$. Hence, in the droplet domain, it must hold that $B=0$ and $F=0$. Conversely, ${\rm 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:

(3.21a)\begin{align} u_{\textrm{D},r} &= \exp\left({- \frac{\eta_\textrm{D}}{ \rho_\textrm{D}} m_\textrm{D}^2 t}\right) \left[ D_\textrm{D} \cos( n_\textrm{D} \theta) - C_\textrm{D} \sin( n_\textrm{D} \theta) \right]\nonumber\\ &\quad\times\frac{n_\textrm{D}}{r} [ A {\rm J}_{n_\textrm{D}} ( m_\textrm{D} r ) + E r^{n_\textrm{D}}], \end{align}

(3.21b)\begin{align} u_{\textrm{D},\theta}&= \exp\left({- \frac{\eta_\textrm{D}}{ \rho_\textrm{D}} m_\textrm{D}^2 t}\right) \left[ C_\textrm{D} \cos( n_\textrm{D} \theta) + D_\textrm{D} \sin( n_\textrm{D} \theta) \right]\nonumber\\ &\quad \times \left[- A \left( m_\textrm{D} {\rm J}_{n_\textrm{D}-1} ( m_\textrm{D} r ) - \frac{n_\textrm{D}}{r} {\rm J}_{n_\textrm{D}} ( m_\textrm{D} r ) \right) - E n_\textrm{D} r^{n_\textrm{D}-1} \right], \end{align}

(3.21c)\begin{align} u_{\textrm{A},r}&= \exp\left({- \frac{\eta_\textrm{A}}{ \rho_\textrm{A}} m_\textrm{A}^2 t}\right) \left[ D_\textrm{A} \cos( n_\textrm{A} \theta) - C_\textrm{A} \sin( n_\textrm{A} \theta) \right]\nonumber\\ &\quad \times\frac{n_\textrm{A}}{r}[B H^{(2)}_{n_\textrm{A}} ( m_\textrm{A} r ) + F r^{{-}n_\textrm{A}}], \end{align}

(3.21d)\begin{align} u_{\textrm{A},\theta}&= \exp\left({- \frac{\eta_\textrm{A}}{ \rho_\textrm{A}} m_\textrm{A}^2 t}\right) \left[ C_\textrm{A} \cos( n_\textrm{A} \theta) + D_\textrm{A} \sin( n_\textrm{A} \theta) \right]\nonumber\\ &\quad \times \left[ - B \left( m_\textrm{A} H^{(2)}_{n_\textrm{A}-1} ( m_\textrm{A} r ) - \frac{n_\textrm{A}}{r} H^{(2)}_{n_\textrm{A}} ( m_\textrm{A} r ) \right) + F n_\textrm{A} r^{{-}n_\textrm{A}-1} \right] , \end{align}

subject to ${\rm Re}(m_i^2)\geqslant {}0$ $(i\in \{\textrm {D},\textrm {A}\})$ and ${\rm Im}(m_{\textrm {A}})<0$.

3.2.3. Pressure solutions

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

(3.22)\begin{equation} p_{\textrm{D}}(r,\theta,t)= \exp\left({- \frac{\eta_\textrm{D}}{ \rho_\textrm{D}} m_\textrm{D}^2 t}\right) r^{\tilde{n}_{\textrm{D}}}(\tilde{D}_{\textrm{D}}\cos(\tilde{n}_{\textrm{D}}\theta)+\tilde{C}_{\textrm{D}}\sin(\tilde{n}_{\textrm{D}}\theta)), \end{equation}

with $\tilde {n}_{\textrm {D}}\in \mathbb {Z}_{\geqslant 0}$ and $\tilde {C}_{\textrm {D}},\tilde {D}_{\textrm {D}}\in \mathbb {C}$ arbitrary constants. By substituting (3.21) and (3.22) into (3.6a), we deduce that compatibility between (3.22) and (3.21) imposes that the index $\tilde {n}_{\textrm {D}}\in \mathbb {Z}_{\geqslant 0}$ and coefficients $\tilde {C}_{\textrm {D}},\tilde {D}_{\textrm {D}}\in \mathbb {C}$ satisfy

(3.23a–c)\begin{equation} \tilde{D}_{\textrm{D}}=D_{\textrm{D}}E\eta_{\textrm{D}}m_{\textrm{D}}^2,\quad \tilde{C}_{\textrm{D}}={-}C_{\textrm{D}}E\eta_{\textrm{D}}m_{\textrm{D}}^2,\quad \tilde{n}_{\textrm{D}}=n_{\textrm{D}}. \end{equation}

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

(3.24a)$$\begin{gather} p_{\textrm{D}} = E \eta_\textrm{D} m_\textrm{D}^2\exp\left({- \frac{\eta_\textrm{D}}{ \rho_\textrm{D}} m_\textrm{D}^2 t}\right) r^{n_\textrm{D}} \left[ D_\textrm{D} \cos(n_\textrm{D} \theta) - C_\textrm{D} \sin(n_\textrm{D} \theta) \right] , \end{gather}$$

(3.24b)$$\begin{gather}p_{\textrm{A}} ={-} F \eta_\textrm{A} m_\textrm{A}^2\exp\left({- \frac{\eta_\textrm{A}}{ \rho_\textrm{A}} m_\textrm{A}^2 t}\right) r^{{-}n_\textrm{A}} \left[ D_\textrm{A} \cos(n_\textrm{A} \theta) - C_\textrm{A} \sin(n_\textrm{A} \theta) \right] . \end{gather}$$

Remark 3.2 It is noteworthy that the Bessel and Hankel functions that appear in the velocity solutions (3.21) are absent from the pressure solutions (3.24). 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 (3.21) and (3.24), involve twelve unknown coefficients: $A$, $B$, $C_\textrm {D}$, $C_\textrm {A}$, $D_\textrm {D}$, $D_\textrm {A}$, $E$, $F$, $m_\textrm {D}$, $m_\textrm {A}$, $n_\textrm {D}$ and $n_\textrm {A}$. We next extract from the general solutions the subspace that complies with the kinematic interface conditions (3.6c) and (3.6d), and the dynamic condition (3.6e), by introducing auxiliary conditions on the coefficients.

3.3.1. Kinematic compatibility conditions

The kinematic conditions (3.6c)–(3.6d) can be equivalently reformulated as

(3.25a)$$\begin{gather} \boldsymbol{u}_\textrm{D} \boldsymbol{\cdot} \boldsymbol{n}_0 = \boldsymbol{u}_\textrm{A} \boldsymbol{\cdot} \boldsymbol{n}_0= \mathcal{V}_1 \quad \textrm{on } \varGamma_0 , \end{gather}$$

(3.25b)$$\begin{gather}\boldsymbol{u}_\textrm{D} \boldsymbol{\cdot} \boldsymbol{t}_0 = \boldsymbol{u}_\textrm{A} \boldsymbol{\cdot} \boldsymbol{t}_0 \quad \textrm{on } \varGamma_0 . \end{gather}$$

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 (3.25a) (respectively (3.25b)) thus pertains to the radial (respectively angular) components of $\boldsymbol {u}$ in (3.21a) and (3.21c) (respectively (3.21b) and (3.21d)). To impose (3.25a), we require the infinitesimal interface velocity $\mathcal {V}_1$ corresponding to (3.4)

(3.26)\begin{equation} \mathcal{V}(\theta,t) = \partial_t{}R_{\delta}(\theta,t) ={-}\delta R_0 \gamma\, {\rm e}^{- \gamma t}( \beta \cos ( k \theta ) + \sqrt{1-\beta^2 } \sin ( k \theta ) ) , \end{equation}

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 (3.25a) that

(3.27) \begin{align} &\exp\left({- \frac{\eta_\textrm{D}}{ \rho_\textrm{D}} m_\textrm{D}^2 t}\right)\left[ D_\textrm{D} \cos( n_\textrm{D} \theta) - C_\textrm{D} \sin( n_\textrm{D} \theta) \right] [ R^{{-}1}_0 A n_\textrm{D} {\rm J}_{n_\textrm{D}} ( m_\textrm{D} R_0 ) + E n_\textrm{D} R_0^{n_\textrm{D}-1} ]\nonumber\\ &\quad= \exp\left({- \frac{\eta_\textrm{A}}{ \rho_\textrm{A}} m_\textrm{A}^2 t}\right)\left[ D_\textrm{A} \cos( n_\textrm{A} \theta) - C_\textrm{A} \sin( n_\textrm{A} \theta) \right]\nonumber\\ &\qquad \times [ R^{{-}1}_0 B n_\textrm{A} H^{(2)}_{n_\textrm{A}} ( m_\textrm{A} R_0 ) + F n_\textrm{A} R_0^{{-}n_\textrm{A}-1}]\nonumber\\ &\quad={-}R_0 \gamma\, {\rm e}^{- \gamma t} (\beta \cos ( k \theta ) + \sqrt{1-\beta^2 }\sin ( k \theta )) . \end{align}

The equalities in (3.27) must hold for all $\theta \in [0,2{\rm \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,

(3.28a)$$\begin{gather} \frac{\eta_\textrm{D}}{\rho_\textrm{D}} m_\textrm{D}^2 = \frac{\eta_\textrm{A}}{\rho_\textrm{A}} m_\textrm{A}^2 = \gamma , \end{gather}$$

(3.28b)$$\begin{gather}n_\textrm{D} = n_\textrm{A} = k , \end{gather}$$

(3.28c)$$\begin{gather}D_\textrm{D} = D_\textrm{A} = \beta , \end{gather}$$

(3.28d)$$\begin{gather}C_\textrm{D} = C_\textrm{A} ={-}\sqrt{1-\beta^2} . \end{gather}$$

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

(3.29a)$$\begin{gather} A k R^{{-}1}_0 {\rm J}_k (m_\textrm{D} R_0) + E k R_0^{k-1} ={-} \gamma R_0 , \end{gather}$$

(3.29b)$$\begin{gather}B k R_0^{{-}1} H^{(2)}_{k} ( m_\textrm{A} R_0 ) + F k R_0^{{-}k-1} ={-} \gamma R_0 , \end{gather}$$

(3.29c) \begin{align} &- A ( m_\textrm{D} {\rm J}_{k-1} ( m_\textrm{D} R_0) - k R_0^{{-}1} {\rm J}_k ( m_\textrm{D} R_0 ) ) -E k R_0^{k-1}\nonumber\\ &\quad + B ( m_\textrm{A} H_{k-1}^{(2)}(m_\textrm{A} R_0) - k R_0^{{-}1} H_k ^{(2)}(m_\textrm{A} R_0 )) - F k R_0^{{-}k-1} = 0 . \end{align}

3.3.2. Dynamic compatibility conditions

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

(3.30a)\begin{align} &- \eta_{\textrm{D}} \boldsymbol{n}_0 (\boldsymbol{\nabla} \boldsymbol{u}_{\textrm{D}} + (\boldsymbol{\nabla} \boldsymbol{u}_{\textrm{D}} )^{\rm T} ) \boldsymbol{n}_0 + p_{\textrm{D}} \nonumber\\ &\quad + \eta_{\textrm{A}} \boldsymbol{n}_0 (\boldsymbol{\nabla} \boldsymbol{u}_{\textrm{A}} + (\boldsymbol{\nabla} \boldsymbol{u}_{\textrm{A}} )^{\rm T} ) \boldsymbol{n}_0 - p_{\textrm{A}} = \sigma_{\textrm{D}\textrm{A}} \kappa_1 \quad \textrm{on } \varGamma , \end{align}

(3.30b)\begin{align} &- \eta_{\textrm{D}} \boldsymbol{t}_0 (\boldsymbol{\nabla} \boldsymbol{u}_{\textrm{D}} + (\boldsymbol{\nabla} \boldsymbol{u}_{\textrm{D}} )^{\rm T} ) \boldsymbol{n}_0 + \eta_{\textrm{A}} \boldsymbol{t}_0 (\boldsymbol{\nabla} \boldsymbol{u}_{\textrm{A}} + (\boldsymbol{\nabla} \boldsymbol{u}_{\textrm{A}} )^{\rm T} ) \boldsymbol{n}_0 = 0 \quad \textrm{on } \varGamma. \end{align}

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

(3.31)\begin{align} \kappa(\theta,t) & = \kappa_0+\delta\kappa_1(\theta,t)+{O}(\delta^2)\nonumber\\ & =R_0^{{-}1} + \delta R^{{-}1}_0 e^{-\gamma t} ( k^2 - 1 ) ( \beta \cos ( k \theta ) - \sqrt{1-\beta^2} \sin ( k \theta)) + {O} (\delta^2 ) . \end{align}

From the expressions for the velocity (3.21) and pressure (3.24), the relation between the coefficients in (3.28), and the dynamic conditions (3.30), it then follows that

(3.32a) \begin{align} &- A \eta_\textrm{D} [ 2 m_\textrm{D} k {R_0}^{{-}1} {\rm J}_{k-1}(m_\textrm{D} {R_0}) - 2 (k+1) k {R_0}^{{-}2} {\rm J}_k (m_\textrm{D} {R_0})]\nonumber\\ &\qquad+ B \eta_\textrm{A} [ 2 m_\textrm{A} k {R_0}^{{-}1} H_{k-1}^{(2)}(m_\textrm{A} {R_0}) - 2 (k+1) k {R_0}^{{-}2} H_k ^{(2)}(m_\textrm{A} {R_0}) ]\nonumber\\ &\qquad+ E\eta_\textrm{D} [ m_\textrm{D}^2 {R_0}^k - 2 (k-1) k {R_0}^{k-2} ]+ F \eta_\textrm{A} [ m_\textrm{A}^2 {R_0}^{{-}k} - 2 (k+1) k {R_0}^{{-}k-2} ]\nonumber\\ &\quad= \sigma_{\textrm{D}\textrm{A}} {R_0}^{{-}1} (k^2 - 1) , \end{align}

(3.32b) \begin{align} &-A \eta_\textrm{D} [ 2 m_\textrm{D} {R_0}^{{-}1} {\rm J}_{k-1}(m_\textrm{D} {R_0}) + ( m_{\textrm{D}}^2 - 2 (k+1) k {R_0}^{{-}2} ) {\rm J}_k (m_\textrm{D} {R_0}) ]\nonumber\\ &\qquad+ B \eta_\textrm{A} [ 2 m_\textrm{A} {R_0}^{{-}1} H_{k-1}^{(2)}(m_\textrm{A} {R_0}) + ( m_\textrm{A}^2 - 2 (k+1) k {R_0}^{{-}2} ) H_k ^{(2)}(m_\textrm{A} {R_0}) ]\nonumber\\ &\qquad+ 2 E \eta_\textrm{D} (k-1) k {R_0}^{k-2} + 2 F \eta_\textrm{A} (k+1) k {R_0}^{{-}k-2}\nonumber\\ &\quad= 0 . \end{align}

3.4. Dispersion relation

For each wavenumber $k\in \mathbb {N}_{\geqslant {}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 (3.6a)–(3.6b), subject to the limit condition (3.2), are given by (3.21) and (3.24). These general solutions contain twelve coefficients. Eight of these coefficients are determined by the kinematic interface condition (3.6c), in accordance with (3.28). The kinematic conditions (3.6c) and (3.6d) imply that the remaining four coefficients, $A,B,E$ and $F$ must satisfy the three identities in (3.29). The dynamic condition (3.6e) demands that, in addition, these four coefficients satisfy the two identities in (3.32). The remaining five conditions on the coefficients can be cast in the form

(3.33)\begin{equation} \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{\mathsf{A}}(k,\gamma)} \begin{pmatrix} A \\ B \\ E \\ F \end{pmatrix} = \underbrace{ \begin{pmatrix} -\gamma R_0 \\ -\gamma R_0 \\ 0 \\ \sigma_{\textrm{D}\textrm{A}} R^{{-}1}_0 (k^2 - 1) \\ 0 \end{pmatrix}}_{{\boldsymbol{b}}(k,\gamma)} \,, \end{equation}

in such a manner that the first three equations in (3.33) represent (3.29) and the latter two represent (3.32). Noting the dependence of (3.29) and (3.32) on the wavenumber $k$, and recalling the dependence of $m_{\textrm {D}}$ and $m_{\textrm {A}}$ in these equation on the temporal response coefficient $\gamma$ via (3.28a), we infer that the entries of ${\boldsymbol{\mathsf{A}}}$ depend on $k$ and $\gamma$. With five constraints and four unknowns, the system of equations (3.33) 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{\mathsf{A}}}(k,\gamma )$. The relation between the existence of a solution and the condition

(3.34)\begin{equation} {\boldsymbol{b}}(k,\gamma)\in\operatorname{span}({\rm{col}}({\boldsymbol{\mathsf{A}}}(k,\gamma))) , \end{equation}

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

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

(3.35)\begin{equation} \det\left((\boldsymbol{\mathsf{A}}\mid\boldsymbol{b})(k,\gamma)\right)=0 , \end{equation}

where $(\boldsymbol{\mathsf{A}}\mid \boldsymbol {b})$ corresponds to $\boldsymbol{\mathsf{A}}$ augmented by $\boldsymbol {b}$. The equivalence between (3.34) and (3.35) follows from the fact that the column vectors of $\boldsymbol{\mathsf{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{\mathsf{A}}$. Moreover, by virtue of the linear independence of the columns of $\boldsymbol{\mathsf{A}}$, if (3.35) holds, then (3.33) has a unique solution. This solution corresponds to the coefficients $(A,B,E,F)$ that, in combination with (3.28), define the droplet and ambient velocity-pressure pairs corresponding to $(k,\gamma )$ according to (3.21) and (3.24).

To facilitate and generalize the root finding of the determinant in (3.35), we non-dimensionalize the matrix entries of the augmented matrix based on the droplet density, $\rho _{\textrm {D}}$, droplet viscosity, $\eta _{\textrm {D}}$, and droplet radius $R_0$. The non-dimensionalized parameters are indicated with a tilde diacritic. Additionally, we introduce the following condensed notation:

(3.36)\begin{equation} \left.\begin{gathered} \mathcal{J} = {\rm J}_k ( \tilde{m}_\textrm{D} (\tilde{\gamma}) ),\quad \mathcal{H} = H^{(2)}_k ( \tilde{m}_\textrm{A} (\tilde{\gamma}) ) ,\quad \zeta = 2(k-1)k , \\ \hat{\mathcal{J}} = \tilde{m}_\textrm{D} (\tilde{\gamma}) {\rm J}_{k-1} ( \tilde{m}_\textrm{D} (\tilde{\gamma}) ) ,\quad \hat{\mathcal{H}} = \tilde{m}_\textrm{A} (\tilde{\gamma}) H^{(2)}_{k-1} ( \tilde{m}_\textrm{A} (\tilde{\gamma}) ) ,\quad \xi = 2 (k+1) k . \end{gathered}\right\} \end{equation}

The non-dimensionalized augmented matrix can then be expressed as

(3.37) \begin{align} & (\tilde{\boldsymbol{\mathsf{A}}}|\tilde{\boldsymbol{b}})(k,\tilde{\gamma})\nonumber\\ &\quad = \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 \\ -2 k \hat{\mathcal{J}} + \xi \mathcal{J} & \tilde{\eta}_\textrm{A} \big[ 2 k \hat{\mathcal{H}} - \xi \mathcal{H} \big] & \tilde{m}^2_\textrm{D} - \zeta & \tilde{\eta}_\textrm{A} \big [ \tilde{m}^2_\textrm{A} - \xi \big] & \tilde{\sigma}_{\textrm{D}\textrm{A}}\left(k^2 - 1\right) \\ -2\hat{\mathcal{J}} - \left(\tilde{m}^2_\textrm{D} - \xi \right) \mathcal{J} & \tilde{\eta}_\textrm{A} \big[ 2 \hat{\mathcal{H}} + \left(\tilde{m}^2_\textrm{A} - \xi \right) \mathcal{H}\big] & \zeta & -\tilde{\eta}_\textrm{A} \xi & 0 \end{pmatrix}. \end{align}

It is not generally feasible to determine the roots of $\det ((\tilde {\boldsymbol{\mathsf{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 (3.37) 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 3.3 The roots of $\det ((\tilde {\boldsymbol{\mathsf{A}}}|\tilde {\boldsymbol {b}})(k,\tilde {\gamma }))$ are not unique: one can infer that

(3.38) \begin{equation} [(\tilde{\boldsymbol{\mathsf{A}}}|\tilde{\boldsymbol{b}})(k,\tilde{\gamma}{}^*)] = [(\tilde{\boldsymbol{\mathsf{A}}}|\tilde{\boldsymbol{b}})(k,\tilde{\gamma})]^*, \end{equation}

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{\mathsf{A}}}|\tilde {\boldsymbol {b}})(k,\tilde {\gamma }))$, then so is $\tilde {\gamma }{}^*$. Since, in addition, it must hold that $\textrm {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{\mathsf{A}}}|\tilde {\boldsymbol {b}})(k,\tilde {\gamma }))$ form a totally disconnected set and, accordingly, for each root there exists a neighbourhood in which that root is unique. A detailed investigation of the uniqueness of the roots of $\det ((\tilde {\boldsymbol{\mathsf{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 (3.4), we obtain the real-valued velocity and pressure fields by taking the real parts of (3.21) and (3.24) after substitution of the relations (3.28), 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 picolitre-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{\mathsf{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 wavenumber 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}\ \mathrm {\mu } \textrm {m}$ suspended in air.

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.

4.1. Set-up and discretization

The oscillating-droplet test cases that we consider pertain to doubly symmetric modes. The set-up of the test cases is similar to that in Demont et al. (Reference Demont, van Zwieten, Diddens and van Brummelen2022, § 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 $\varOmega =(0,50)^2\ \mathrm {\mu }\textrm {m}^2$ and prescribe symmetry conditions on $\varGamma _{{sym}}:=\{(x_1,x_2)\in \partial \varOmega :\{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 (3.1a), while the diffuse-interface model exhibits a moving interface, we prescribe auxiliary conditions in accordance with an initially circular droplet. Specifically, with reference to (3.4), we select $t_0$ such that $-\textrm {Im}(\gamma )t_0={\rm \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, $\varGamma _{{ext}}:=\partial \varOmega \setminus \varGamma _{{sym}}$, is furnished with Dirichlet conditions for velocity and homogeneous Neumann conditions for the order parameter and the chemical potential

(4.1) \begin{equation} \left.\begin{array}{c@{}} \boldsymbol{u}({\cdot},t)=\delta\boldsymbol{u}_{\textrm{A}}({\cdot},t_0+t)\\ \partial_n\varphi=0\\ \partial_n\mu=0 \end{array}\right\}\quad\text{on}\ \varGamma_{{ext}}, \text{ for }t\in[0,T), \end{equation}

where $\delta \boldsymbol {u}_{\textrm {A}}$ corresponds to the ambient velocity solution (3.21) 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.

(4.2)\begin{equation} \varphi(\boldsymbol{x},0)=\varphi_0(\boldsymbol{x}):=\tanh\left(\frac{d_{{\pm}}(\boldsymbol{x},\varGamma_0)}{\sqrt{2}\varepsilon}\right), \end{equation}

where $d_{\pm }(\boldsymbol {x},\varGamma _0)$ represents the signed distance from $\boldsymbol {x}$ to $\varGamma _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 (4.2) is meta-stable if $\varepsilon$ is sufficiently small compared with the radius of curvature of $\varGamma _0$. In conjunction with (4.2), we impose the following initial condition for velocity:

(4.3)\begin{equation} \boldsymbol{u}(\boldsymbol{x},0)= \begin{cases} \delta\boldsymbol{u}_{\textrm{D}}(\boldsymbol{x},t_0) & \text{if }\boldsymbol{x}\in\varOmega_{\textrm{D},0},\\ \delta\boldsymbol{u}_{\textrm{A}}(\boldsymbol{x},t_0) & \text{if }\boldsymbol{x}\in\varOmega_{\textrm{A},0}, \end{cases} \end{equation}

where the data in the right member of (4.3) correspond to the velocity solutions according to (3.21) 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 with 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 Demont et al. (Reference Demont, van Zwieten, Diddens and van Brummelen2022). For coherence, we present a concise overview of the numerical methodology. The weak form of the NSCH equations (2.1) is discretized with respect to the spatial dependence with $P3-P2$ (Taylor–Hood) $\mathcal {C}^0$ truncated hierarchical B-splines (see Hughes, Cottrell & Bazilevs Reference Hughes, Cottrell and Bazilevs2005; Cottrell, Hughes & Bazilevs Reference Cottrell, Hughes and Bazilevs2009; Giannelli, Jüttler & Speleers Reference Giannelli, Jüttler and Speleers2012) for the velocity and pressure fields, and $P3$ $\mathcal {C}^0$ truncated hierarchical B-splines for the order parameter and chemical potential; see van Brummelen et al. (Reference van Brummelen, Demont and van Zwieten2021, § 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 (Dörfler Reference Dörfler1996; Bertoluzza et al. Reference Bertoluzza, Nochetto, Quarteroni, Siebert and Veeser2012). 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 Demont et al. (Reference Demont, van Zwieten, Diddens and van Brummelen2022) and van Brummelen et al. (Reference van Brummelen, Demont and van Zwieten2021) 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\ \mathrm {\mu }\textrm {m}$, and we perform $L_{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_{max}$ the original parameter values for $\varepsilon$ and $m$ are used. A skew-symmetric formulation according to Layton (Reference Layton2008) 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 (Wu, van Zwieten & van der Zee Reference Wu, van Zwieten and van der Zee2014) 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. John, Matthies & Rang (Reference John, Matthies and Rang2006). For the temporal discretization, we employ a time-step size $\tau = 2^{-7} \times 10\ \mathrm {\mu }\textrm {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 the generalized minimal residual method (GMRES) with a preconditioner based on a partition of the NSCH system into NS and CH subsystems; see Demont et al. (Reference Demont, van Zwieten, Diddens and van Brummelen2022) for further details.

To illustrate the set-up 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\ \mathrm {\mu } \textrm {m}$, mobility $m=9.5272\times 10^{-13} \textrm {m}^2\textrm {s}\ \textrm {kg}^{-1}$, maximum number of refinement levels $L_{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 (respectively bottom) half of each panel displays the velocity (respectively pressure) field. The right (respectively left) half of each panel depicts the diffuse-interface simulation (respectively 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

(4.4)

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 _{max}$ relative to the baseline interface thickness $\varepsilon _{max} = 2^{-7}\times 10^{2}\ \mathrm {\mu } \textrm {m} = 7.8125\times 10^{-1}\ \mathrm {\mu } \textrm {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_{max}$ with $m_{{max}} = 2.4389632\times 10^{-10}\ \textrm {m}^d \textrm {s}\ \textrm {kg}^{-1}$. The range of mobility parameters has been determined empirically such that $[2^{-12},1]m_{{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 ${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 (4.4), for $k=2$ for a half-period of oscillation.

Table 4.

Deviation between the diffuse-interface solution and the sharp-interface solution according to (4.4), for $k=4$ for a half-period of oscillation.

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$ vs $\varepsilon$. For both wavenumbers, we observe an optimal scaling $m\propto \varepsilon ^{a_{opt}}$ with $a_{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_{opt}}$ larger than the optimal mobility for $k=2$. This suggests that the optimal mobility in fact scales with $(\varepsilon /\ell )^{a_{opt}}$, where $\ell$ represents another characteristic length scale of the interface which, for the considered droplet-oscillation test case, is proportional to the wavelength 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_{opt}}$ in the sharp-interface limit. The results in § 4.2 also convey that the constant of proportionality in the scaling relation $m_{\varepsilon }=\mathscr {C}\varepsilon ^{a_{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 (respectively figure 6) plots for each $\varepsilon \in \mathcal {E}$ the ratio of the deviation $\operatorname {dev}(\varepsilon,m)$ in the columns of able 3 (respectively table 4) to the minimal deviation $\operatorname {dev}(\varepsilon,m_{opt,\varepsilon })$ vs the ratio $m/m_{opt,\varepsilon }$. Noting that the curves in figures 5 and 6 exhibit a vanishing slope near $m/m_{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_{opt,\varepsilon },m_{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_{opt,\varepsilon })$ vs normalized mobility $m/m_{opt,\varepsilon }$ for $k=2$.

Figure 6.

Normalized deviation $\operatorname {dev}(\varepsilon,m)/\operatorname {dev}(\varepsilon,m_{opt,\varepsilon })$ vs normalized mobility $m/m_{opt,\varepsilon }$ for $k=4$.

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

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 _{max})^a \, m_0)$ vs $\varepsilon$ for $a\in \{0,\ldots,3\}$. Herein, $m_0$ corresponds to the optimal sampled mobility for $\varepsilon _{max}$, viz. $m_0=2^{-3}m_{max}$ for $k=2$ and $m_0=2^{-2}m_{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_{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_{opt,\varepsilon })={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 _{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 non-classical sharp-interface Navier–Stokes/Mullins–Sekerka model; see Abels & Garcke (Reference Abels and Garcke2018).

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$.