1 Introduction
The dynamics of oscillating drops is of interest to researchers in pattern formation and dynamical systems, as well as having practical applications over a wide variety of scales in areas as diverse as astroseismology, containerless material processing for high purity crystal growth and drug delivery and mixing in microfluidic devices.
Surface tension is responsible for the spherical shape of a drop. In the absence of external forces, if the drop is slightly perturbed, it will recover its spherical shape through decaying oscillations. This problem was first considered by Kelvin (Reference Kelvin1863) and Rayleigh (Reference Rayleigh1879), who described natural oscillations of drops of inviscid fluids. Rayleigh (Reference Rayleigh1879) and Lamb (Reference Lamb1932) derived the, now classic, resonance mode frequency resulting from the restoring force of surface tension:
$$\begin{eqnarray}\unicode[STIX]{x1D714}^{2}=\frac{\unicode[STIX]{x1D70E}}{\unicode[STIX]{x1D70C}}\frac{\ell (\ell -1)(\ell +2)}{R^{3}},\end{eqnarray}$$
where
$\unicode[STIX]{x1D714}$
is the frequency,
$\unicode[STIX]{x1D70E}$
and
$\unicode[STIX]{x1D70C}$
the surface tension and density,
$R$
is the radius and
$\ell$
is the degree of the spherical harmonic
$$\begin{eqnarray}Y_{\ell }^{m}=P_{\ell }^{m}(\cos \unicode[STIX]{x1D703})\text{e}^{\text{i}m\unicode[STIX]{x1D719}}\end{eqnarray}$$
describing the perturbation. Linear analyses including viscosity were carried out by Reid (Reference Reid1960), Chandrasekhar (Reference Chandrasekhar1961) and Miller & Scriven (Reference Miller and Scriven1968). These authors demonstrated the equivalence of this problem to that of a fluid globe oscillating under the influence of self-gravitation, generalizing the previous conclusion of Lamb. Chandrasekhar showed that the return to a spherical shape could take place via monotonic decay as well as via damped oscillations. The problem was further investigated by Prosperetti (Reference Prosperetti1980) using an initial-value code. Weakly nonlinear effects in inviscid fluid drops were investigated by Tsamopoulos & Brown (Reference Tsamopoulos and Brown1983) using a Poincaré–Lindstedt expansion technique.
The laboratory realization of any configuration with only spherically symmetric radially directed forces is difficult. Indeed such experiments have been sent into space, e.g. (Wang, Anilkumar & Lee Reference Wang, Anilkumar and Lee1996; Futterer et al. Reference Futterer, Krebs, Plesa, Zaussinger, Hollerbach, Breuer and Egbers2013) and in parabolic flight (Falcón et al. Reference Falcón, Falcon, Bortolozzo and Fauve2009) in order to eliminate or reduce the perturbing influence of the gravitational field of the Earth. Wang et al. (Reference Wang, Anilkumar and Lee1996) were able to confirm the decrease in frequency with increasing oscillation amplitude predicted by Tsamopoulos & Brown (Reference Tsamopoulos and Brown1983). Wang et al. (Reference Wang, Anilkumar and Lee1996) mention, however, that the treatment of viscosity is not exact. Falcón et al. (Reference Falcón, Falcon, Bortolozzo and Fauve2009) produced spherical capillary wave turbulence and compared its spectrum with theoretical predictions.
In the laboratory, drops have been levitated by using acoustic or magnetic forces and excited by periodic electric modulation (Shen, Xie & Wei Reference Shen, Xie and Wei2010); drops and puddles weakly pinned on a vibrating substrate (Brunet & Snoeijer Reference Brunet and Snoeijer2011) have produced star-like patterns. One of the purposes of such experiments is to provide a measurement of the surface tension. Trinh, Zwern & Wang (Reference Trinh, Zwern and Wang1982) visualized the shapes and internal flow of vibrating drops and compared the frequencies to those of Lamb (Reference Lamb1932) and the damping coefficients to those derived by Marston (Reference Marston1980). These experimental procedures cannot produce a perfectly spherical base state, and indeed, Trinh & Wang (Reference Trinh and Wang1982) and Cummings & Blackburn (Reference Cummings and Blackburn1991) discuss differences between oscillating oblate and prolate drops, and the resulting deviations from (1.1). Because the experimentally observed frequencies remain close to (1.1), it seems likely that the results of our stability analysis are also only mildly affected by a departure from perfect spherical symmetry.
Here, and in a companion paper, we consider a viscous drop under the influence of a time-periodic radial bulk force and of surface tension. Our investigation relies on a variety of mathematical and computational tools. Here, we solve the linear stability problem by adapting to spherical coordinates the Floquet method of Kumar & Tuckerman (Reference Kumar and Tuckerman1994). At the linear level, the instability depends only on the spherical wavenumber
$\ell$
of (1.2) as illustrated by the Lamb–Rayleigh relation (1.1). Thus, perturbations which are not axisymmetric (
$m\neq 0$
in (1.2)) have the same thresholds as the corresponding axisymmetric (
$m=0$
) perturbations. Indeed, the theoretical and numerical investigations listed above have assumed that the drop shape remains axisymmetric. The fully nonlinear Faraday problem, however, usually leads to patterns which are non-axisymmetric. In our complementary investigation (Ebo Adou et al.
Reference Ebo Adou, Tuckerman, Shin, Chergui and Juric2016), we will describe the results of full three-dimensional simulations which calculate the interface motion and the velocity field inside and outside the parametrically forced drop and interpret them in the context of the theory of pattern formation.
2 Governing equations
2.1 Equations of motion
We consider a drop of viscous, incompressible liquid bounded by a spherical free surface that separates the liquid from the exterior in the presence of an uniform radial oscillatory body force. The fluid motion inside the drop satisfies the Navier–Stokes equations
$$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D70C}\left[\frac{\unicode[STIX]{x2202}}{\unicode[STIX]{x2202}t}+(\boldsymbol{U}\boldsymbol{\cdot }\unicode[STIX]{x1D735})\right]\boldsymbol{U}=-\unicode[STIX]{x1D735}P+\unicode[STIX]{x1D702}\unicode[STIX]{x1D6FB}^{2}\boldsymbol{U}-\unicode[STIX]{x1D70C}G(r,t)\boldsymbol{e}_{r}, & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D735}\boldsymbol{\cdot }\boldsymbol{U}=0, & \displaystyle\end{eqnarray}$$
$\boldsymbol{U}$
is the velocity,
$P$
the pressure,
$\unicode[STIX]{x1D70C}$
the density and
$\unicode[STIX]{x1D702}$
the dynamic viscosity.
$G(r,t)$
is an imposed radial parametric acceleration given by 
$$\begin{eqnarray}G(r,t)=(g-a\cos (\unicode[STIX]{x1D714}t))\frac{r}{R},\end{eqnarray}$$
which is regular at the origin.
The interface is located at
$$\begin{eqnarray}r=R+\unicode[STIX]{x1D701}(\unicode[STIX]{x1D703},\unicode[STIX]{x1D719},t).\end{eqnarray}$$
Conservation of volume leads to the requirement that the integral of
$\unicode[STIX]{x1D701}$
over the sphere must be zero.
The boundary conditions applied at the interface are the kinematic condition, which states that the interface is advected by the fluid
$$\begin{eqnarray}\left[\frac{\unicode[STIX]{x2202}}{\unicode[STIX]{x2202}t}+(\boldsymbol{U}\boldsymbol{\cdot }\unicode[STIX]{x1D735})\right]\unicode[STIX]{x1D701}=U_{r}|_{r=R+\unicode[STIX]{x1D701}}\end{eqnarray}$$
and the interface stress balance equation, which is given by
$$\begin{eqnarray}\boldsymbol{n}\boldsymbol{\cdot }\hat{\unicode[STIX]{x1D72B}}-\boldsymbol{n}\boldsymbol{\cdot }\unicode[STIX]{x1D72B}=\unicode[STIX]{x1D70E}\boldsymbol{n}(\unicode[STIX]{x1D735}\boldsymbol{\cdot }\boldsymbol{n})-\unicode[STIX]{x1D735}\unicode[STIX]{x1D70E},\end{eqnarray}$$
where
$\unicode[STIX]{x1D70E}$
is the surface tension coefficient and
$\boldsymbol{n}$
represents the unit outward normal to the surface, both defined only on the interface.
The tensors
$\unicode[STIX]{x1D72B}$
(drop) and
$\hat{\unicode[STIX]{x1D72B}}$
(medium) denote the stress tensor in each fluid and are defined by
$$\begin{eqnarray}\unicode[STIX]{x1D72B}=-P\boldsymbol{I}_{\boldsymbol{d}}+\unicode[STIX]{x1D702}[\unicode[STIX]{x1D735}\boldsymbol{U}+(\unicode[STIX]{x1D735}\boldsymbol{U})^{\text{T}}].\end{eqnarray}$$
For simplicity, we consider a situation in which the outer medium has no effect on the drop, and so we set the density, pressure and stress tensor
$\hat{\unicode[STIX]{x1D72B}}$
outside the drop to zero. The boundary conditions corresponding to the case of drop forced in a medium are described in the appendix. We assume that the surface tension is uniform, so
$\unicode[STIX]{x1D735}\unicode[STIX]{x1D70E}=0$
. The tangential stress balance equation at the free surface then reduces to
$$\begin{eqnarray}\boldsymbol{n}\boldsymbol{\cdot }\unicode[STIX]{x1D72B}\boldsymbol{\cdot }\boldsymbol{t}=0\end{eqnarray}$$
for both unit tangent vectors
$\boldsymbol{t}$
.
The normal stress jump boundary condition determines the curvature of the deformed interface. The Laplace formula relates the normal stress jump to the divergence of the normal field, which is in turn equal to the mean curvature:
$$\begin{eqnarray}-\boldsymbol{n}\boldsymbol{\cdot }\unicode[STIX]{x1D72B}\boldsymbol{\cdot }\boldsymbol{n}=\unicode[STIX]{x1D70E}\unicode[STIX]{x1D735}\boldsymbol{\cdot }\boldsymbol{n}=\unicode[STIX]{x1D70E}\left(\frac{1}{R_{1}}+\frac{1}{R_{2}}\right),\end{eqnarray}$$
with
$R_{1}$
and
$R_{2}$
the principal radii of curvature at a given point of the surface.
For a sphere,
$R_{1}=R_{2}=R$
and (2.8) becomes
$$\begin{eqnarray}P|_{r=R}=2\frac{\unicode[STIX]{x1D70E}}{R}\end{eqnarray}$$
and the solution to the governing equations (2.1a
) and boundary conditions (2.4), (2.7), (2.8) is the motionless equilibrium spherical state at
$r=R$
with
$$\begin{eqnarray}\displaystyle & \displaystyle \bar{\boldsymbol{U}}=0, & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle \bar{P}(r,t)=2\frac{\unicode[STIX]{x1D70E}}{R}-\int _{r}^{R}\unicode[STIX]{x1D70C}G(r^{\prime },t)\,\text{d}r^{\prime }, & \displaystyle\end{eqnarray}$$
$\bar{P}$
is continuously differentiable at the origin because
$G(0,t)=0$
.
2.2 Linearizing the governing equations
We linearize the Navier–Stokes equations about the unperturbed state (2.10) by decomposing the velocity and the pressure
$$\begin{eqnarray}\displaystyle & \displaystyle \boldsymbol{U}=\bar{\boldsymbol{U}}+\boldsymbol{u}, & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle P=\bar{P}+p, & \displaystyle\end{eqnarray}$$
$\boldsymbol{u}$
and
$p$
$$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D70C}\frac{\unicode[STIX]{x2202}\boldsymbol{u}}{\unicode[STIX]{x2202}t}=-\unicode[STIX]{x1D735}p+\unicode[STIX]{x1D702}\unicode[STIX]{x1D6FB}^{2}\boldsymbol{u}, & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D735}\boldsymbol{\cdot }\boldsymbol{u}=0. & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D735}_{H}\boldsymbol{\cdot }\equiv \frac{1}{r\sin \unicode[STIX]{x1D703}}\frac{\unicode[STIX]{x2202}}{\unicode[STIX]{x2202}\unicode[STIX]{x1D703}}\sin \unicode[STIX]{x1D703}\,\hat{e}_{\unicode[STIX]{x1D703}}\boldsymbol{\cdot }+\frac{1}{r\sin \unicode[STIX]{x1D703}}\frac{\unicode[STIX]{x2202}}{\unicode[STIX]{x2202}\unicode[STIX]{x1D719}}\hat{e}_{\unicode[STIX]{x1D719}}, & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D735}_{H}\equiv \frac{\hat{e}_{\unicode[STIX]{x1D703}}}{r}\frac{\unicode[STIX]{x2202}}{\unicode[STIX]{x2202}\unicode[STIX]{x1D703}}+\frac{\hat{e}_{\unicode[STIX]{x1D719}}}{r\sin \unicode[STIX]{x1D703}}\frac{\unicode[STIX]{x2202}}{\unicode[STIX]{x2202}\unicode[STIX]{x1D719}}, & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D6FB}_{H}^{2}\equiv \frac{1}{r\sin \unicode[STIX]{x1D703}}\frac{\unicode[STIX]{x2202}}{\unicode[STIX]{x2202}\unicode[STIX]{x1D703}}\left(\sin \unicode[STIX]{x1D703}\frac{\unicode[STIX]{x2202}}{\unicode[STIX]{x2202}\unicode[STIX]{x1D703}}\right)+\frac{1}{r\sin \unicode[STIX]{x1D703}}\frac{\unicode[STIX]{x2202}}{\unicode[STIX]{x2202}\unicode[STIX]{x1D719}}. & \displaystyle\end{eqnarray}$$
For a solenoidal field satisfying (2.12b ), definitions (2.13c ) lead to
$$\begin{eqnarray}\displaystyle (\unicode[STIX]{x1D6FB}^{2}\boldsymbol{u})_{r} & = & \displaystyle \left(\frac{1}{r^{2}}\frac{\unicode[STIX]{x2202}}{\unicode[STIX]{x2202}r}r^{2}\frac{\unicode[STIX]{x2202}}{\unicode[STIX]{x2202}r}+\unicode[STIX]{x1D6FB}_{H}^{2}-\frac{2}{r^{2}}\right)u_{r}-\frac{2}{r}\unicode[STIX]{x1D735}_{H}\boldsymbol{\cdot }\boldsymbol{u}_{H}\nonumber\\ \displaystyle & = & \displaystyle \left(\frac{1}{r^{2}}\frac{\unicode[STIX]{x2202}}{\unicode[STIX]{x2202}r}r^{2}\frac{\unicode[STIX]{x2202}}{\unicode[STIX]{x2202}r}-\frac{2}{r^{2}}+\frac{2}{r^{3}}\frac{\unicode[STIX]{x2202}}{\unicode[STIX]{x2202}r}r^{2}+\unicode[STIX]{x1D6FB}_{H}^{2}\right)u_{r}\nonumber\\ \displaystyle & = & \displaystyle \left(\frac{1}{r^{3}}\frac{\unicode[STIX]{x2202}}{\unicode[STIX]{x2202}r}r^{2}\frac{\unicode[STIX]{x2202}}{\unicode[STIX]{x2202}r}r+\unicode[STIX]{x1D6FB}_{H}^{2}\right)u_{r}\equiv \widetilde{\mathscr{L}}^{2}u_{r}.\end{eqnarray}$$
We can then eliminate the horizontal velocity
$\boldsymbol{u}_{H}=(u_{\unicode[STIX]{x1D703}},u_{\unicode[STIX]{x1D719}})$
and the pressure
$p$
from (2.12a
) in the usual way by operating with
$\boldsymbol{e}_{r}\boldsymbol{\cdot }\unicode[STIX]{x1D735}\times \unicode[STIX]{x1D735}\times$
on (2.12a
), leading to:
$$\begin{eqnarray}\widetilde{\mathscr{L}}^{2}\left(\frac{\unicode[STIX]{x2202}}{\unicode[STIX]{x2202}t}-\unicode[STIX]{x1D708}\widetilde{\mathscr{L}}^{2}\right)u_{r}=0.\end{eqnarray}$$
Since we are interested in the linear stability of the interface located at
$r=R+\unicode[STIX]{x1D701}$
with
$\unicode[STIX]{x1D701}\ll R$
, we Taylor expand the fields and their radial derivatives around
$r=R$
and retain only the lowest-order terms, which are evaluated at
$r=R$
. The kinematic condition (2.4) becomes
$$\begin{eqnarray}\frac{\unicode[STIX]{x2202}}{\unicode[STIX]{x2202}t}\unicode[STIX]{x1D701}=u_{r}|_{r=R}.\end{eqnarray}$$
We now wish to apply the stress balance equations at
$r=R$
. The components of the stress tensor which we will need are
$$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D6F1}_{r\unicode[STIX]{x1D703}}=\unicode[STIX]{x1D702}\left(\frac{1}{r}\frac{\unicode[STIX]{x2202}u_{r}}{\unicode[STIX]{x2202}\unicode[STIX]{x1D703}}+r\frac{\unicode[STIX]{x2202}}{\unicode[STIX]{x2202}r}\left(\frac{u_{\unicode[STIX]{x1D703}}}{r}\right)\right), & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D6F1}_{r\unicode[STIX]{x1D719}}=\unicode[STIX]{x1D702}\left(\frac{1}{r\sin \unicode[STIX]{x1D703}}\frac{\unicode[STIX]{x2202}u_{r}}{\unicode[STIX]{x2202}\unicode[STIX]{x1D719}}+r\frac{\unicode[STIX]{x2202}}{\unicode[STIX]{x2202}r}\left(\frac{u_{\unicode[STIX]{x1D719}}}{r}\right)\right), & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D6F1}_{rr}=2\frac{\unicode[STIX]{x2202}u_{r}}{\unicode[STIX]{x2202}r}. & \displaystyle\end{eqnarray}$$
$r=R$
, leading to: 
$$\begin{eqnarray}\unicode[STIX]{x1D6F1}_{r\unicode[STIX]{x1D703}}|_{r=R}=\unicode[STIX]{x1D6F1}_{r\unicode[STIX]{x1D719}}|_{r=R}=0.\end{eqnarray}$$
Taking the horizontal divergence of
$\unicode[STIX]{x1D6F1}_{r\unicode[STIX]{x1D703}}\,\hat{e}_{\unicode[STIX]{x1D703}}+\unicode[STIX]{x1D6F1}_{r\unicode[STIX]{x1D719}}\,\hat{e}_{\unicode[STIX]{x1D719}}$
leads to
$$\begin{eqnarray}\displaystyle 0 & = & \displaystyle [\unicode[STIX]{x1D735}_{H}\boldsymbol{\cdot }(\unicode[STIX]{x1D6F1}_{r\unicode[STIX]{x1D703}}\,\hat{e}_{\unicode[STIX]{x1D703}}+\unicode[STIX]{x1D6F1}_{r\unicode[STIX]{x1D719}}\,\hat{e}_{\unicode[STIX]{x1D719}})]_{r=R}\nonumber\\ \displaystyle & = & \displaystyle \unicode[STIX]{x1D702}\left[\unicode[STIX]{x1D735}_{H}\boldsymbol{\cdot }\left(\frac{1}{r}\frac{\unicode[STIX]{x2202}u_{r}}{\unicode[STIX]{x2202}\unicode[STIX]{x1D703}}\,\hat{e}_{\unicode[STIX]{x1D703}}+\frac{1}{r\sin \unicode[STIX]{x1D703}}\frac{\unicode[STIX]{x2202}u_{r}}{\unicode[STIX]{x2202}\unicode[STIX]{x1D719}}\,\hat{e}_{\unicode[STIX]{x1D719}}\right)+\unicode[STIX]{x1D735}_{H}\boldsymbol{\cdot }\left(r\frac{\unicode[STIX]{x2202}}{\unicode[STIX]{x2202}r}\left(\frac{u_{\unicode[STIX]{x1D703}}}{r}\right)\hat{e}_{\unicode[STIX]{x1D703}}+r\frac{\unicode[STIX]{x2202}}{\unicode[STIX]{x2202}r}\left(\frac{u_{\unicode[STIX]{x1D719}}}{r}\right)\hat{e}_{\unicode[STIX]{x1D719}}\right)\right]_{r=R}\nonumber\\ \displaystyle & = & \displaystyle \unicode[STIX]{x1D702}\left[\unicode[STIX]{x1D735}_{H}\boldsymbol{\cdot }\unicode[STIX]{x1D735}_{H}u_{r}+\frac{\unicode[STIX]{x2202}}{\unicode[STIX]{x2202}r}\unicode[STIX]{x1D735}_{H}\boldsymbol{\cdot }\left(u_{\unicode[STIX]{x1D703}}\hat{e}_{\unicode[STIX]{x1D703}}+u_{\unicode[STIX]{x1D719}}\hat{e}_{\unicode[STIX]{x1D719}}\right)\right]_{r=R}\!=\unicode[STIX]{x1D702}\left[\left(\unicode[STIX]{x1D6FB}_{H}^{2}-\frac{\unicode[STIX]{x2202}}{\unicode[STIX]{x2202}r}\frac{1}{r^{2}}\frac{\unicode[STIX]{x2202}}{\unicode[STIX]{x2202}r}r^{2}\right)u_{r}\right]_{r=R}\!,~\end{eqnarray}$$
which is the form of the tangential stress continuity equation that we impose.
We now wish to linearize the normal stress balance equation:
$$\begin{eqnarray}-[\boldsymbol{n}\boldsymbol{\cdot }\unicode[STIX]{x1D72B}\boldsymbol{\cdot }\boldsymbol{n}]_{r=R+\unicode[STIX]{x1D701}}=[\boldsymbol{n}\boldsymbol{\cdot }(P\boldsymbol{I}_{\boldsymbol{d}}-\unicode[STIX]{x1D702}[\unicode[STIX]{x1D735}\boldsymbol{U}+(\unicode[STIX]{x1D735}\boldsymbol{U})^{\text{T}}])\boldsymbol{\cdot }\boldsymbol{n}]_{r=R+\unicode[STIX]{x1D701}}=\unicode[STIX]{x1D70E}\left(\frac{1}{R_{1}}+\frac{1}{R_{2}}\right).\end{eqnarray}$$
The right-hand side of (2.19) is
$\unicode[STIX]{x1D70E}$
times the curvature of a deformed interface and can be shown (Lamb Reference Lamb1932, § 275) to be, up to first order in
$\unicode[STIX]{x1D701}$
,
$$\begin{eqnarray}\unicode[STIX]{x1D70E}\left(\frac{1}{R_{1}}+\frac{1}{R_{2}}\right)_{r=R+\unicode[STIX]{x1D701}}\approx \frac{2\unicode[STIX]{x1D70E}}{R}-\unicode[STIX]{x1D70E}\left(\frac{2}{R^{2}}+\unicode[STIX]{x1D6FB}_{H}^{2}\right)\unicode[STIX]{x1D701}.\end{eqnarray}$$
For the left-hand side of (2.19), we use (2.10) to expand the pressure as
$$\begin{eqnarray}(\bar{P}+p)_{r=R+\unicode[STIX]{x1D701}}\approx (\bar{P}+p)_{r=R}+\left(\frac{\unicode[STIX]{x2202}\bar{P}}{\unicode[STIX]{x2202}r}\right)_{r=R}\unicode[STIX]{x1D701}=2\frac{\unicode[STIX]{x1D70E}}{R}+p|_{r=R}-\unicode[STIX]{x1D70C}G(R,t)\unicode[STIX]{x1D701}.\end{eqnarray}$$
Adding the term resulting from the viscosity leads to
$$\begin{eqnarray}-\boldsymbol{n}\boldsymbol{\cdot }\unicode[STIX]{x1D72B}\boldsymbol{\cdot }\boldsymbol{n}|_{r=R+\unicode[STIX]{x1D701}}\approx 2\frac{\unicode[STIX]{x1D70E}}{R}+p|_{r=R}-\unicode[STIX]{x1D70C}G(R,t)\unicode[STIX]{x1D701}-2\unicode[STIX]{x1D702}\left(\frac{\unicode[STIX]{x2202}u_{r}}{\unicode[STIX]{x2202}r}\right)_{r=R}.\end{eqnarray}$$
Setting (2.20) equal to (2.22) leads to the desired linearized form:
$$\begin{eqnarray}p|_{r=R}-\unicode[STIX]{x1D70C}G(R,t)\unicode[STIX]{x1D701}-2\unicode[STIX]{x1D702}\left(\frac{\unicode[STIX]{x2202}u_{r}}{\unicode[STIX]{x2202}r}\right)_{r=R}=-\unicode[STIX]{x1D70E}\left(\frac{2}{R^{2}}+\unicode[STIX]{x1D6FB}_{H}^{2}\right)\unicode[STIX]{x1D701}.\end{eqnarray}$$
It will be useful to take the horizontal Laplacian of (2.23):
$$\begin{eqnarray}\unicode[STIX]{x1D6FB}_{H}^{2}p|_{r=R}=2\unicode[STIX]{x1D702}\unicode[STIX]{x1D6FB}_{H}^{2}\frac{\unicode[STIX]{x2202}}{\unicode[STIX]{x2202}r}u_{r}|_{r=R}+\unicode[STIX]{x1D70C}G(R,t)\unicode[STIX]{x1D6FB}_{H}^{2}\unicode[STIX]{x1D701}-\unicode[STIX]{x1D70E}\unicode[STIX]{x1D6FB}_{H}^{2}\left(\frac{2}{R^{2}}+\unicode[STIX]{x1D6FB}_{H}^{2}\right)\unicode[STIX]{x1D701}.\end{eqnarray}$$
We can derive another expression for
$\unicode[STIX]{x1D6FB}_{H}^{2}p|_{r=R}$
by taking the horizontal divergence of (2.12a
):
$$\begin{eqnarray}\unicode[STIX]{x1D6FB}_{H}^{2}p=\frac{1}{r^{2}}\left(\unicode[STIX]{x1D70C}\frac{\unicode[STIX]{x2202}}{\unicode[STIX]{x2202}t}-\unicode[STIX]{x1D702}\unicode[STIX]{x1D6FB}^{2}\right)\frac{\unicode[STIX]{x2202}}{\unicode[STIX]{x2202}r}(r^{2}u_{r}).\end{eqnarray}$$
Setting equal the right-hand sides of equations (2.25) and (2.24), we obtain the pressure jump condition
$$\begin{eqnarray}\left[\frac{1}{r^{2}}\left(\unicode[STIX]{x1D70C}\frac{\unicode[STIX]{x2202}}{\unicode[STIX]{x2202}t}-\unicode[STIX]{x1D702}\unicode[STIX]{x1D6FB}^{2}\right)\frac{\unicode[STIX]{x2202}}{\unicode[STIX]{x2202}r}(r^{2}u_{r})-2\unicode[STIX]{x1D702}\unicode[STIX]{x1D6FB}_{H}^{2}\frac{\unicode[STIX]{x2202}}{\unicode[STIX]{x2202}r}u_{r}\right]_{r=R}=\unicode[STIX]{x1D70C}G(R,t)\unicode[STIX]{x1D6FB}_{H}^{2}\unicode[STIX]{x1D701}-\unicode[STIX]{x1D70E}\unicode[STIX]{x1D6FB}_{H}^{2}\left(\frac{2}{R^{2}}+\unicode[STIX]{x1D6FB}_{H}^{2}\right)\unicode[STIX]{x1D701}.\end{eqnarray}$$
This is the only equation in which the parametrical external forcing appears explicitly; note that only the value
$G(r=R,t)$
on the sphere is relevant. We now have a set of linear equations (2.14), (2.15), (2.18) and (2.26) involving only
$u_{r}(r,\unicode[STIX]{x1D703},\unicode[STIX]{x1D719},t)$
and
$\unicode[STIX]{x1D701}(\unicode[STIX]{x1D703},\unicode[STIX]{x1D719},t)$
, which reduce to those for the planar case (Kumar & Tuckerman Reference Kumar and Tuckerman1994) in the limit of
$R\rightarrow \infty$
.
3 Solution to linear stability problem
3.1 Spherical harmonic decomposition
The equations simplify somewhat when we use the poloidal–toroidal decomposition
$$\begin{eqnarray}\boldsymbol{u}=\unicode[STIX]{x1D735}\times (f_{T}\boldsymbol{e}_{r})+\unicode[STIX]{x1D735}\times \unicode[STIX]{x1D735}\times (f\boldsymbol{e}_{r}).\end{eqnarray}$$
The radial velocity component
$u_{r}$
depends only on the poloidal field
$f$
and is given by
$$\begin{eqnarray}u_{r}(r,\unicode[STIX]{x1D703},\unicode[STIX]{x1D719},t)=-\unicode[STIX]{x1D6FB}_{H}^{2}f(r,\unicode[STIX]{x1D703},\unicode[STIX]{x1D719},t).\end{eqnarray}$$
Using
$$\begin{eqnarray}\underbrace{\left(\frac{1}{r^{3}}\frac{\unicode[STIX]{x2202}}{\unicode[STIX]{x2202}r}r^{2}\frac{\unicode[STIX]{x2202}}{\unicode[STIX]{x2202}r}r+\unicode[STIX]{x1D6FB}_{H}^{2}\right)}_{\widetilde{\mathscr{L}}^{2}}\unicode[STIX]{x1D6FB}_{H}^{2}=\unicode[STIX]{x1D6FB}_{H}^{2}\underbrace{\left(\frac{\unicode[STIX]{x2202}^{2}}{\unicode[STIX]{x2202}r^{2}}+\unicode[STIX]{x1D6FB}_{H}^{2}\right)}_{\mathscr{L}^{2}}\end{eqnarray}$$
we express (2.14) in terms of the poloidal field
$$\begin{eqnarray}\unicode[STIX]{x1D6FB}_{H}^{2}\left(\frac{\unicode[STIX]{x2202}}{\unicode[STIX]{x2202}t}-\unicode[STIX]{x1D708}\mathscr{L}^{2}\right)\mathscr{L}^{2}f=0.\end{eqnarray}$$
Functions are expanded in series of spherical harmonics
$Y_{\ell }^{m}(\unicode[STIX]{x1D703},\unicode[STIX]{x1D719})=P_{\ell }^{m}(\cos \unicode[STIX]{x1D703})\text{e}^{\text{i}m\unicode[STIX]{x1D719}}$
satisfying
$$\begin{eqnarray}\unicode[STIX]{x1D6FB}_{H}^{2}Y_{\ell }^{m}(\unicode[STIX]{x1D703},\unicode[STIX]{x1D719})=-\frac{\ell (\ell +1)}{r^{2}}Y_{\ell }^{m}(\unicode[STIX]{x1D703},\unicode[STIX]{x1D719}).\end{eqnarray}$$
We write the deviation of the interface
$\unicode[STIX]{x1D701}(t,\unicode[STIX]{x1D703},\unicode[STIX]{x1D719})$
and the scalar function
$f$
as
$$\begin{eqnarray}\unicode[STIX]{x1D701}(t,\unicode[STIX]{x1D703},\unicode[STIX]{x1D719})=\mathop{\sum }_{\ell =1}^{\infty }\mathop{\sum }_{m=-\ell }^{\ell }\unicode[STIX]{x1D701}_{\ell }^{m}(t)Y_{\ell }^{m}(\unicode[STIX]{x1D703},\unicode[STIX]{x1D719})\quad \text{and}\quad f(r,\unicode[STIX]{x1D703},\unicode[STIX]{x1D719},t)=\mathop{\sum }_{\ell =1}^{\infty }\mathop{\sum }_{m=-\ell }^{\ell }f_{\ell }^{m}(r,t)Y_{\ell }^{m}(\unicode[STIX]{x1D703},\unicode[STIX]{x1D719}).\end{eqnarray}$$
The equations do not couple different spherical modes
$(\ell ,m)$
, allowing us to consider each mode separately. The term multiplying
$\unicode[STIX]{x1D70E}\unicode[STIX]{x1D6FB}_{H}^{2}\unicode[STIX]{x1D701}$
in the normal stress equation (2.26) becomes
$$\begin{eqnarray}\left(\frac{2}{R^{2}}+\unicode[STIX]{x1D6FB}_{H}^{2}\right)\unicode[STIX]{x1D701}_{\ell }^{m}=\left(\frac{2}{R^{2}}-\frac{\ell (\ell +1)}{R^{2}}\right)\unicode[STIX]{x1D701}_{\ell }^{m}=-\frac{(\ell -1)(\ell +2)}{R^{2}}\unicode[STIX]{x1D701}_{\ell }^{m}.\end{eqnarray}$$
The value
$\ell =0$
corresponds to an overall expansion or contraction of the sphere, which is forbidden by mass conservation and is therefore excluded from (3.6). The value
$\ell =1$
corresponds to a shift of the drop, rather than a deformation of the interface, so that surface tension cannot act as a restoring force; it is included in this study only when the surface tension
$\unicode[STIX]{x1D70E}$
is zero and the constant radial bulk force
$g$
is non-zero.
The complete linear problem given by (3.4), (2.15), (2.18) and (2.26) is expressed in terms of
$f_{\ell }^{m}(r,t)$
and
$\unicode[STIX]{x1D701}_{\ell }^{m}(t)$
as:
$$\begin{eqnarray}\displaystyle & \displaystyle \left(\frac{\unicode[STIX]{x2202}}{\unicode[STIX]{x2202}t}-\unicode[STIX]{x1D708}\mathscr{L}_{\ell }^{2}\right)\mathscr{L}_{\ell }^{2}f_{\ell }^{m}=0, & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle \frac{\unicode[STIX]{x2202}}{\unicode[STIX]{x2202}t}\unicode[STIX]{x1D701}_{\ell }^{m}=\frac{\ell (\ell +1)}{R^{2}}f_{\ell }^{m}|_{r=R}, & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle \left(\mathscr{L}_{\ell }^{2}-\frac{2}{r}\frac{\unicode[STIX]{x2202}}{\unicode[STIX]{x2202}r}\right)f_{\ell }^{m}|_{r=R}=0, & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & & \displaystyle \left[\left(\unicode[STIX]{x1D70C}\frac{\unicode[STIX]{x2202}}{\unicode[STIX]{x2202}t}\frac{\unicode[STIX]{x2202}}{\unicode[STIX]{x2202}r}-\unicode[STIX]{x1D702}\left(\frac{\unicode[STIX]{x2202}^{3}}{\unicode[STIX]{x2202}r^{3}}+\frac{2}{r}\frac{\unicode[STIX]{x2202}^{2}}{\unicode[STIX]{x2202}r^{2}}-\ell (\ell +1)\left(\frac{3}{r^{2}}\frac{\unicode[STIX]{x2202}}{\unicode[STIX]{x2202}r}-\frac{4}{r^{3}}\right)\right)\right)f_{\ell }^{m}\right]_{r=R}\nonumber\\ \displaystyle & & \displaystyle \quad =-\left(\unicode[STIX]{x1D70C}G(R,t)+\unicode[STIX]{x1D70E}\frac{(\ell -1)(\ell +2)}{R^{2}}\right)\unicode[STIX]{x1D701}_{\ell }^{m},\end{eqnarray}$$
where we have used (3.2), (3.5) and
$$\begin{eqnarray}\mathscr{L}_{\ell }^{2}\equiv \frac{\unicode[STIX]{x2202}^{2}}{\unicode[STIX]{x2202}r^{2}}-\frac{\ell (\ell +1)}{r^{2}}\end{eqnarray}$$
and have divided through by
$\ell (\ell +1)/R^{2}$
.
The value of
$m$
does not appear in these equations, in much the same way as the Cartesian linear Faraday problem depends only on the wavenumber
$k$
and not on its orientation.
3.2 Floquet solution
The presence of the time-periodic term in (3.11) means that (3.8), (3.9), (3.10), (3.11) comprise a Floquet problem. To solve it, we follow the procedure of Kumar & Tuckerman (Reference Kumar and Tuckerman1994), whereby
$\unicode[STIX]{x1D701}_{\ell }^{m}$
and
$f_{\ell }^{m}$
are written in the Floquet form:
$$\begin{eqnarray}\unicode[STIX]{x1D701}_{\ell }^{m}(t)=\text{e}^{(\unicode[STIX]{x1D707}+\text{i}\unicode[STIX]{x1D6FC})t}\mathop{\sum }_{n}\unicode[STIX]{x1D701}_{n}\text{e}^{\text{i}n\unicode[STIX]{x1D714}t}\quad \text{and}\quad f_{\ell }^{m}(r,t)=\text{e}^{(\unicode[STIX]{x1D707}+\text{i}\unicode[STIX]{x1D6FC})t}\mathop{\sum }_{n}f_{n}(r)\text{e}^{\text{i}n\unicode[STIX]{x1D714}t},\end{eqnarray}$$
where
$\unicode[STIX]{x1D707}+\text{i}\unicode[STIX]{x1D6FC}$
is the Floquet exponent and we have omitted the indices
$(\ell ,m)$
. Substituting the Floquet expansions (3.13) into (3.8) gives, for each temporal frequency
$n$
$$\begin{eqnarray}(\unicode[STIX]{x1D707}+\text{i}(n\unicode[STIX]{x1D714}+\unicode[STIX]{x1D6FC})-\unicode[STIX]{x1D708}\mathscr{L}_{\ell }^{2})\mathscr{L}_{\ell }^{2}f_{n}=0\end{eqnarray}$$
or
$$\begin{eqnarray}(\mathscr{L}_{\ell }^{2}-q_{n}^{2})\mathscr{L}_{\ell }^{2}f_{n}(r)=0,\end{eqnarray}$$
where
$$\begin{eqnarray}q_{n}^{2}\equiv \frac{\unicode[STIX]{x1D707}+\text{i}(n\unicode[STIX]{x1D714}+\unicode[STIX]{x1D6FC})}{\unicode[STIX]{x1D708}}.\end{eqnarray}$$
In order to solve the fourth-order differential equation (3.15), we first solve the homogeneous second-order equation
$$\begin{eqnarray}(\mathscr{L}_{\ell }^{2}-q_{n}^{2})\tilde{f}_{n}=0,\end{eqnarray}$$
which is a modified spherical Bessel, or Riccati–Bessel, equation (Abramowitz & Stegun Reference Abramowitz and Stegun1965). The general solutions are of the form
$$\begin{eqnarray}\tilde{f}_{n}(r)=\tilde{B}_{n}r^{1/2}J_{\ell +1/2}(iq_{n}r)+\tilde{D}_{n}^{2}r^{1/2}J_{-(\ell +1/2)}(iq_{n}r),\end{eqnarray}$$
where
$J_{\ell +1/2}$
is the spherical Bessel function of half-integer order
$\ell +1/2$
.
The remaining second-order differential equation is
$$\begin{eqnarray}\mathscr{L}_{\ell }^{2}f_{n}=\tilde{f}_{n}\end{eqnarray}$$
whose solutions are the general solutions of (3.15) and are given by
$$\begin{eqnarray}f_{n}(r)=A_{n}r^{\ell +1}+B_{n}r^{1/2}J_{\ell +1/2}(\text{i}q_{n}r)+C_{n}r^{-\ell }+D_{n}r^{1/2}J_{-(\ell +1/2)}(\text{i}q_{n}r).\end{eqnarray}$$
Note that
$u_{r}\sim f_{n}/r^{2}$
satisfies (2.14). Eliminating the solutions in (3.20) which diverge at the centre, we are left with:
$$\begin{eqnarray}f_{n}(r)=A_{n}r^{\ell +1}+B_{n}r^{1/2}J_{\ell +1/2}(\text{i}q_{n}r).\end{eqnarray}$$
The constants
$A_{n}$
and
$B_{n}$
can be related to
$\unicode[STIX]{x1D701}_{n}$
by using the kinematic condition (3.9) and the tangential stress condition (3.10) which, for Floquet mode
$n$
, are
$$\begin{eqnarray}\displaystyle & \displaystyle (\unicode[STIX]{x1D707}+\text{i}(n\unicode[STIX]{x1D714}+\unicode[STIX]{x1D6FC}))\unicode[STIX]{x1D701}_{n}=\frac{\ell (\ell +1)}{R^{2}}f_{n}|_{r=R}, & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle \left(\frac{\unicode[STIX]{x2202}^{2}}{\unicode[STIX]{x2202}r^{2}}-\frac{2}{r}\frac{\unicode[STIX]{x2202}}{\unicode[STIX]{x2202}r}+\frac{\ell (\ell +1)}{r^{2}}\right)f_{n}|_{r=R}=0. & \displaystyle\end{eqnarray}$$
Appendix A gives more details on the determination of these constants and also presents the solution and boundary conditions in the case for which the exterior of the drop is a fluid rather than a vacuum.
There remains the normal stress (pressure jump) condition (3.11), the only equation which couples temporal Floquet modes for different
$n$
. Writing
$$\begin{eqnarray}a\cos (\unicode[STIX]{x1D714}t)\mathop{\sum }_{n}\unicode[STIX]{x1D701}_{n}\text{e}^{\text{i}n\unicode[STIX]{x1D714}t}=a\frac{\text{e}^{\text{i}\unicode[STIX]{x1D714}t}+\text{e}^{-\text{i}\unicode[STIX]{x1D714}t}}{2}\mathop{\sum }_{n}\unicode[STIX]{x1D701}_{n}\text{e}^{\text{i}n\unicode[STIX]{x1D714}t}=\frac{a}{2}\mathop{\sum }_{n}(\unicode[STIX]{x1D701}_{n+1}+\unicode[STIX]{x1D701}_{n-1})\text{e}^{\text{i}n\unicode[STIX]{x1D714}t}\end{eqnarray}$$
and inserting the Floquet decomposition (3.13) into (3.11) leads to
$$\begin{eqnarray}\displaystyle & & \displaystyle \left[\left(\unicode[STIX]{x1D70C}(\unicode[STIX]{x1D707}+\text{i}(n\unicode[STIX]{x1D714}+\unicode[STIX]{x1D6FC}))\frac{\unicode[STIX]{x2202}}{\unicode[STIX]{x2202}r}-\unicode[STIX]{x1D702}\left(\frac{\unicode[STIX]{x2202}^{3}}{\unicode[STIX]{x2202}r^{3}}+\frac{2}{r}\frac{\unicode[STIX]{x2202}^{2}}{\unicode[STIX]{x2202}r^{2}}-\ell (\ell +1)\left(\frac{3}{r^{2}}\frac{\unicode[STIX]{x2202}}{\unicode[STIX]{x2202}r}-\frac{4}{r^{3}}\right)\right)\right)f_{n}\right]_{r=R}\nonumber\\ \displaystyle & & \displaystyle \quad +\left(\unicode[STIX]{x1D70C}g+\unicode[STIX]{x1D70E}\frac{(\ell -1)(\ell +2)}{R^{2}}\right)\unicode[STIX]{x1D701}_{n}=\unicode[STIX]{x1D70C}\frac{a}{2}(\unicode[STIX]{x1D701}_{n+1}+\unicode[STIX]{x1D701}_{n-1}).\end{eqnarray}$$
Using (A 9) and (A 10) to express the partial derivatives of
$f_{n}$
as multiples of
$\unicode[STIX]{x1D701}_{n}$
, (3.25) can be written as an eigenvalue problem with eigenvalues
$a$
and eigenvectors
$\{\unicode[STIX]{x1D701}_{n}\}$
$$\begin{eqnarray}\mathscr{A}\unicode[STIX]{x1D701}=a\mathscr{B}\unicode[STIX]{x1D701}.\end{eqnarray}$$
The usual procedure for a stability analysis is to fix the wavenumber (here the spherical mode
$\ell$
) and the forcing amplitude
$a$
, to calculate the exponents
$\unicode[STIX]{x1D707}+\text{i}\unicode[STIX]{x1D6FC}$
of the growing solutions and to select that whose growth rate
$\unicode[STIX]{x1D707}(\ell ,a)$
is largest. Following Kumar & Tuckerman (Reference Kumar and Tuckerman1994), we instead fix
$\unicode[STIX]{x1D707}=0$
and restrict consideration to two kinds of growing solutions, harmonic with
$\unicode[STIX]{x1D6FC}=0$
and subharmonic with
$\unicode[STIX]{x1D6FC}=\unicode[STIX]{x1D714}/2$
. We then solve the problem (3.26) for the eigenvalues
$a$
and eigenvectors
$\{\unicode[STIX]{x1D701}_{n}\}$
and select the smallest, or several smallest, real positive eigenvalues
$a$
. These give the marginal stability curves
$a(\ell ,\unicode[STIX]{x1D707}=0,\unicode[STIX]{x1D6FC}=\unicode[STIX]{x1D714}/2)$
and
$a(\ell ,\unicode[STIX]{x1D707}=0,\unicode[STIX]{x1D6FC}=0)$
without interpolation. Our computations require no more than 10 Fourier modes in representation (3.13). This method can be used to solve any Floquet problem for which the overall amplitude of the time-periodic terms can be varied. The detailed procedure for the solution of the eigenvalue problem (3.26) is described in Kumar & Tuckerman (Reference Kumar and Tuckerman1994), Kumar (Reference Kumar1996).
4 Ideal fluid case and non-dimensionalization
For an ideal fluid drop (
$\unicode[STIX]{x1D708}=0$
) and for a given spherical harmonic
$Y_{\ell }^{m}$
, system (3.8) reduces to
$$\begin{eqnarray}\frac{\unicode[STIX]{x2202}}{\unicode[STIX]{x2202}t}\mathscr{L}_{\ell }^{2}f(r,t)=0.\end{eqnarray}$$
We make the customary assumption that
$\mathscr{L}_{\ell }^{2}f(r,t)$
is not only constant, as implied by (4.1), but zero. In this case, the solution which does not diverge at the drop centre is of the form
$$\begin{eqnarray}f(r,t)=F(t)r^{\ell +1}.\end{eqnarray}$$
As the tangential stress is purely viscous in origin, only the kinematic condition (3.9) and the pressure jump condition across the interface (3.11) are applied. Using (4.2), these are reduced to:
$$\begin{eqnarray}\displaystyle & \displaystyle \dot{\unicode[STIX]{x1D701}}(t)=\frac{\ell (\ell +1)}{R^{2}}F(t)R^{\ell +1}, & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle (\ell +1){\dot{F}}(t)R^{\ell }=-\left(G(R,t)+\frac{\unicode[STIX]{x1D70E}}{\unicode[STIX]{x1D70C}}\frac{(\ell -1)(\ell +2)}{R^{2}}\right)\unicode[STIX]{x1D701}_{\ell }^{m}(t). & \displaystyle\end{eqnarray}$$
By differentiating (4.3) with respect to time and substituting into (4.4), we arrive at
$$\begin{eqnarray}\ddot{\unicode[STIX]{x1D701}}=-\left(g\frac{\ell }{R}+\frac{\unicode[STIX]{x1D70E}}{\unicode[STIX]{x1D70C}}\frac{\ell (\ell -1)(\ell +2)}{R^{3}}-a\frac{\ell }{R}\cos (\unicode[STIX]{x1D714}t)\right)\unicode[STIX]{x1D701}.\end{eqnarray}$$
Focusing for the moment on the unforced equation, we define:
$$\begin{eqnarray}\unicode[STIX]{x1D714}_{0}^{2}\equiv \left(g\frac{\ell }{R}+\frac{\unicode[STIX]{x1D70E}}{\unicode[STIX]{x1D70C}}\frac{\ell (\ell -1)(\ell +2)}{R^{3}}\right).\end{eqnarray}$$
Equation (4.6) is the spherical analogue of the usual dispersion relation for gravity–capillary waves in a plane layer of infinite depth, with the modifications
$$\begin{eqnarray}\displaystyle & \displaystyle gk\rightarrow g\frac{\ell }{R}, & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle \frac{\unicode[STIX]{x1D70E}}{\unicode[STIX]{x1D70C}}k^{3}\rightarrow \frac{\unicode[STIX]{x1D70E}}{\unicode[STIX]{x1D70C}}\frac{\ell (\ell -1)(\ell +2)}{R^{3}}. & \displaystyle\end{eqnarray}$$
$\ell$
is
$2\unicode[STIX]{x03C0}R/\ell$
. The transformation (4.7b
) must be understood in light of the fact that an unperturbed sphere, unlike a planar surface, has a non-zero curvature term proportional to
$2\unicode[STIX]{x1D70E}/R$
, from which (4.7b
) is derived as a deviation via (3.7). In the absence of a bulk force,
$g=0$
, (4.6) becomes the formula of Rayleigh (Reference Rayleigh1879) or Lamb (Reference Lamb1932) for the eigenfrequencies of capillary oscillation of a spherical drop perturbed by a deformation characterized by spherical wavenumber
$\ell$
.
Using definitions (4.6) and
$$\begin{eqnarray}a_{0}\equiv \frac{R\unicode[STIX]{x1D714}_{0}^{2}}{\ell }\end{eqnarray}$$
and non-dimensionalizing time via
$\hat{t}\equiv t\unicode[STIX]{x1D714}$
, equation (4.5) can be converted to the Mathieu equation
$$\begin{eqnarray}\frac{\text{d}^{2}\unicode[STIX]{x1D701}}{\text{d}\hat{t}^{2}}=-\left(\frac{\unicode[STIX]{x1D714}_{0}}{\unicode[STIX]{x1D714}}\right)^{2}\left(1-\frac{a}{a_{0}}\cos \hat{t}\right)\unicode[STIX]{x1D701}\end{eqnarray}$$
combining the multiple parameters
$g$
,
$R$
,
$\unicode[STIX]{x1D70E}$
,
$\unicode[STIX]{x1D70C}$
,
$a$
,
$\unicode[STIX]{x1D714}$
and
$\ell$
into only two,
$\unicode[STIX]{x1D714}/\unicode[STIX]{x1D714}_{0}$
and
$a/a_{0}$
. The stability regions of (4.9) are bounded by tongues which intersect the
$a=0$
axis at
$$\begin{eqnarray}\unicode[STIX]{x1D714}=\frac{2}{n}\unicode[STIX]{x1D714}_{0}.\end{eqnarray}$$
Thus the inviscid spherical Faraday linear stability problem reduces to the Mathieu equation, as it does in the planar case (Benjamin & Ursell Reference Benjamin and Ursell1954). In a Faraday wave experiment or numerical simulation,
$\ell$
, unlike the other parameters, is not known a priori. Instead, in light of (4.6), equation (4.10) is a cubic equation which determines
$\ell$
given the other parameters. Since
$\unicode[STIX]{x1D714}_{0}$
and
$a_{0}$
are functions of
$\ell$
, both of the variables
$\unicode[STIX]{x1D714}/\unicode[STIX]{x1D714}_{0}$
and
$a/a_{0}$
contain the unknown
$\ell$
and so cannot be interpreted as simple non-dimensionalizations of
$\unicode[STIX]{x1D714}$
and
$a$
. (For the purely gravitational case,
$a_{0}=g$
is independent of
$\ell$
.)
Figure 1. Instability tongues of an inviscid fluid drop due to oscillatory forcing with amplitude
$a$
and angular frequency
$\unicode[STIX]{x1D714}$
. Solid curves bound subharmonic tongues and dashed curves bound harmonic tongues. (a) Tongues corresponding to gravitational instability with spherical wavenumbers
$\ell =1$
, 2, 3, 4 originate at
$\unicode[STIX]{x1D714}/\unicode[STIX]{x1D714}_{g}=2\sqrt{\ell }/n$
. (b) Tongues corresponding to capillary instability with spherical wavenumbers
$\ell =2$
, 3, 4, 5 originate at
$\unicode[STIX]{x1D714}/\unicode[STIX]{x1D714}_{c}=2\sqrt{\ell (\ell -1)(\ell +2)}/n$
, with
$n$
odd (even) for subharmonic (harmonic) tongues.
It is useful to examine (4.6) and (4.8) in the two limits of gravity and capillary waves. We first define non-dimensional angular frequencies and oscillation amplitudes which do not depend on
$\ell$
:
$$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D714}_{g}^{2}\equiv \frac{g}{R},\quad a_{g}\equiv R\unicode[STIX]{x1D714}_{g}^{2}=g, & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D714}_{c}^{2}\equiv \frac{\unicode[STIX]{x1D70E}}{\unicode[STIX]{x1D70C}R^{3}},\quad a_{c}\equiv R\unicode[STIX]{x1D714}_{c}^{2}, & \displaystyle\end{eqnarray}$$
so that
$$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D714}_{0}^{2}=\unicode[STIX]{x1D714}_{g}^{2}\,\ell +\unicode[STIX]{x1D714}_{c}^{2}\,\ell (\ell -1)(\ell +2), & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle a_{0}=a_{g}+a_{c}(\ell -1)(\ell +2). & \displaystyle\end{eqnarray}$$
The Bond number measuring the relative importance of the two forces can be written as:
$$\begin{eqnarray}Bo\equiv \frac{\unicode[STIX]{x1D70C}gR^{2}}{\unicode[STIX]{x1D70E}}=\frac{\unicode[STIX]{x1D714}_{g}^{2}}{\unicode[STIX]{x1D714}_{c}^{2}}.\end{eqnarray}$$
In the gravity-dominated regime (large
$Bo$
), we write (4.10) as
$$\begin{eqnarray}\unicode[STIX]{x1D714}^{2}=\left(\frac{2}{n}\right)^{2}\unicode[STIX]{x1D714}_{g}^{2}\left[\ell +\frac{1}{Bo}\ell (\ell -1)(\ell +2)\right].\end{eqnarray}$$
The stability boundaries for
$1/Bo=0$
are given in figure 1(a). The subharmonic and harmonic tongues originate at
$\unicode[STIX]{x1D714}/\unicode[STIX]{x1D714}_{g}=2\sqrt{\ell }$
and
$\unicode[STIX]{x1D714}/\unicode[STIX]{x1D714}_{g}=\sqrt{\ell }$
.
In the capillary-dominated regime (small
$Bo$
) it is more appropriate to write (4.10) as
$$\begin{eqnarray}\unicode[STIX]{x1D714}^{2}=\left(\frac{2}{n}\right)^{2}\unicode[STIX]{x1D714}_{c}^{2}\left[\ell (\ell -1)(\ell +2)+Bo\,\ell \right].\end{eqnarray}$$
The stability boundaries for
$Bo=0$
are given in figure 1(b) in terms of
$\unicode[STIX]{x1D714}/\unicode[STIX]{x1D714}_{c}$
and
$a/a_{c}$
. These consist of families of tongues, which originate on the
$a=0$
axis at
$\unicode[STIX]{x1D714}/\unicode[STIX]{x1D714}_{c}=2\sqrt{\ell (\ell -1)(\ell +2)}/n$
for
$\ell =2$
, 3, 4, 5 and for
$n=1,2,\ldots$
within which the drop has one of the spatial forms corresponding to the spherical wavenumber
$\ell$
and oscillates like
$\text{e}^{\text{i}n\unicode[STIX]{x1D714}t/2}$
. The solid curves bound the first subharmonic instability tongues, which originate on the
$a=0$
axis at
$\unicode[STIX]{x1D714}/\unicode[STIX]{x1D714}_{c}=2\sqrt{\ell (\ell -1)(\ell +2)}$
for
$\ell =2$
, 3, 4, 5, within which the drop oscillates like
$\text{e}^{\text{i}\unicode[STIX]{x1D714}t/2}$
. The dashed curves bound the first harmonic tongues, originating at
$\unicode[STIX]{x1D714}/\unicode[STIX]{x1D714}_{c}=\sqrt{\ell (\ell -1)(\ell +2)}$
describing a response like
$\text{e}^{\text{i}\unicode[STIX]{x1D714}t}$
. Tongues for higher
$n$
are located at still lower values of
$\unicode[STIX]{x1D714}$
.
The curves in figure 1 for different
$\ell$
,
$g$
,
$\unicode[STIX]{x1D70E}/\unicode[STIX]{x1D70C}$
,
$R$
all collapse onto a single set of tongues when they are plotted in terms of
$\unicode[STIX]{x1D714}/\unicode[STIX]{x1D714}_{0}$
and
$a/a_{0}$
. This is shown in figure 2, in which the various tongues correspond to the temporal harmonic index
$n$
. In order to use figure 2 to determine whether the drop is stable against perturbations with spherical wavenumber
$\ell$
when a radial force with amplitude
$a$
and angular frequency
$\unicode[STIX]{x1D714}$
is applied, the following procedure must be used.
Figure 2. Same data as in figure 1, but scaled by
$\unicode[STIX]{x1D714}_{0}$
and
$a_{0}$
. The tongues for the gravitational and the capillary cases and for all values of
$\ell$
all collapse onto a single set of tongues. Solid curves bound subharmonic tongues and dashed curves bound harmonic tongue.
For each
$\ell$
, formulas (4.6) and (4.8) are used to determine the values of
$(\unicode[STIX]{x1D714}_{0},a_{0})$
. If
$(\unicode[STIX]{x1D714}/\unicode[STIX]{x1D714}_{0},a/a_{0})$
is inside one of the instability tongues (usually, but not always, that corresponding to an
$\unicode[STIX]{x1D714}/2$
response with
$n=1$
) then the drop is unstable to perturbations of that
$\ell$
. The drop is stable if
$(\unicode[STIX]{x1D714}/\unicode[STIX]{x1D714}_{0},a/a_{0})$
lies outside the tongues for all
$\ell$
and all
$n$
. This is the procedure described by Benjamin & Ursell (Reference Benjamin and Ursell1954) and carried out by Batson, Zoueshtiagh & Narayanan (Reference Batson, Zoueshtiagh and Narayanan2013) in a cylindrical geometry.
Because
$(\unicode[STIX]{x1D714}/\unicode[STIX]{x1D714}_{0},a/a_{0})$
depends on the unknown
$\ell$
, an experimental or numerical choice of parameters cannot immediately be assigned to a point in figure 2, rendering its interpretation somewhat more obscure. It is perhaps because of this that investigations of the Faraday instability are often presented in dimensional terms. Figures like 1(a,b), in which the two axes are non-dimensional quantities defined in terms of known parameters, represent a good compromise between universality and ease of use. This treatment can also be applied to non-spherical geometries in which the unperturbed surface is flat and the depth is finite.
5 Viscous fluids and scaling laws
We now return to the viscous case. For a viscous fluid, it is not possible to reduce the governing equations to a Mathieu equation even with the addition of a damping term (Kumar & Tuckerman Reference Kumar and Tuckerman1994). As described in § 3.2, the governing equations and boundary conditions are reduced to an eigenvalue problem, whose solution gives the critical oscillation amplitude
$a$
as an eigenvalue.
Figure 3. Instability tongues due to oscillatory forcing of amplitude
$a$
and angular frequency
$\unicode[STIX]{x1D714}$
for a viscous fluid drop. Solid curves bound subharmonic tongues and dashed curves bound harmonic tongues. (a) Tongues corresponding to gravitational instability with viscosity
$\unicode[STIX]{x1D708}/\sqrt{gR^{3}}=0.08$
. Spherical wavenumbers are
$\ell =1$
(magenta), 2 (blue), 3 (red), 4 (green), 5 (black). (b) Tongues correspond to capillary instability with viscosity
$\unicode[STIX]{x1D708}/\sqrt{\unicode[STIX]{x1D70E}R/\unicode[STIX]{x1D70C}}=0.08$
. Spherical wavenumbers are
$\ell =2$
(blue), 3 (red), 4 (green), 5 (black), 6 (purple).
Figure 3 displays the regions of instability of a viscous drop using the same conventions as we did for the ideal fluid case, i.e. treating capillary and gravitational instability separately and plotting the stability boundaries in units of
$\unicode[STIX]{x1D714}_{c}$
,
$\unicode[STIX]{x1D714}_{g}$
,
$a_{c}$
,
$a_{g}$
. Viscosity smoothes the instability tongues and displaces the critical forcing amplitude towards higher values, with a displacement which increases with
$\ell$
. The viscosity used in figure 3 is
$\unicode[STIX]{x1D708}/\sqrt{\unicode[STIX]{x1D70E}R/\unicode[STIX]{x1D70C}}=0.08$
and
$\unicode[STIX]{x1D708}/\sqrt{gR^{3}}=0.08$
, the Ohnesorge number and the inverse square root of the Galileo number for the capillary and the gravitational cases, respectively. This value is chosen to be high enough to show the important qualitative effect of viscosity, but low enough to permit the first few tongues to be shown on a single diagram. It can be seen that the effect of viscosity is greater on tongues with higher
$\ell$
; this important point will be discussed below. For low enough frequency, i.e.
$\unicode[STIX]{x1D714}/\unicode[STIX]{x1D714}_{c}\lesssim 1.2$
in figure 3(a) and
$\unicode[STIX]{x1D714}/\unicode[STIX]{x1D714}_{g}\lesssim 4$
in figure 3(b), it can be seen that the instability is harmonic rather than subharmonic, as discussed by Kumar (Reference Kumar1996).
Figure 4. Viscosity dependence of
$\unicode[STIX]{x1D714}$
(rows 1 and 2) and
$a$
(rows 3 and 4) at threshold. Gravitational (rows 1 and 3) and capillary (rows 2 and 4) Faraday instability shown for
$\ell =2$
(blue), 4 (green) and 6 (purple). Left column (a,d,g,j)
$\unicode[STIX]{x1D714}$
and
$a$
non-dimensionalized in terms of experimentally imposed quantities
$\unicode[STIX]{x1D714}_{g}$
,
$\unicode[STIX]{x1D714}_{c}$
,
$a_{g}$
and
$a_{c}$
. Non-dimensionalization of middle column (b,e,h,k) uses (5.1) and (5.2). It can be seen that
$a\sim \unicode[STIX]{x1D708}$
for
$\unicode[STIX]{x1D708}$
small. Right column (c,f,i,l) shows analogous quantities for the planar infinite-depth case, calculated using (5.1b
) and (5.2b
); the scaling is seen to be exact in this case.
Figure 4 shows the variation with viscosity of the Faraday threshold, more specifically of the coordinates
$(\unicode[STIX]{x1D714},a)$
of the minimum of the primary subharmonic tongue, for
$\ell =2,4,6$
for the capillary and gravitational Faraday instability. The first column shows this dependence using non-dimensional quantities that are independent of
$\ell$
. (We explained the motivation for such a choice in § 4, namely that
$\ell$
is not known a priori in an experiment or full numerical simulation.) The second column shows the non-dimensionalization that best fits the data for all values of
$\ell$
. The appropriate choice for the non-dimensionalization of viscosity is
$$\begin{eqnarray}\tilde{\unicode[STIX]{x1D708}}\equiv \frac{\unicode[STIX]{x1D708}\ell (\ell +1)}{R^{2}\unicode[STIX]{x1D714}_{0}}\quad \text{or}\quad \tilde{\unicode[STIX]{x1D708}}\equiv \frac{\unicode[STIX]{x1D708}k^{2}}{\unicode[STIX]{x1D714}_{0}}\end{eqnarray}$$
in the spherical or planar geometries, respectively, based on the wavelength. See Bechhoefer et al. (Reference Bechhoefer, Ego, Manneville and Johnson1995), who studied the influence of viscosity on Faraday waves. The choice of (5.1) is guided by comparing the viscous and oscillatory terms in (3.25) and corroborated by the fact that the curves in the second column have only a weak dependence on
$\ell$
and on
$\tilde{\unicode[STIX]{x1D708}}$
. We recall that the ratio of the horizontal (angular) wavelength
$2\unicode[STIX]{x03C0}R/\ell$
to the depth (radius)
$R$
goes to zero as
$\ell$
increases, and the curvature of the sphere is less manifested over a horizontal wavelength. With increasing
$\ell$
, the curves for the spherical case can be seen to approach the corresponding quantities in the planar infinite-depth case, shown in the third column, despite the fact that we are far from the large
$\ell$
limit.
Figure 4(a–f) shows that the frequency
$\unicode[STIX]{x1D714}$
which favours waves (corresponding to the bottom of the tongue) is a non-monotonic function of
$\unicode[STIX]{x1D708}$
, for which an explanation is proposed by Kumar & Tuckerman (Reference Kumar and Tuckerman1994). At lower viscosities, the flow is assumed to be irrotational, as in (4.2), and equation (4.6) is modified merely by subtracting a term proportional to
$\unicode[STIX]{x1D708}^{2}$
. This leads to a decrease in the critical
$\unicode[STIX]{x1D714}$
from
$2\unicode[STIX]{x1D714}_{0}$
. At higher viscosities, it is assumed that the response time
$4\unicode[STIX]{x03C0}/\unicode[STIX]{x1D714}$
approaches the viscous time scale, here
$O(\ell (\ell +1)R^{2}/\unicode[STIX]{x1D708})$
, leading to an increase in
$\unicode[STIX]{x1D714}$
with
$\unicode[STIX]{x1D708}$
when the other parameters are fixed. Experimental values for
$\tilde{\unicode[STIX]{x1D708}}$
are, however, rarely greater than one.
Concerning the oscillation amplitude
$a$
, we find that the appropriate choice for non-dimensionalization is
$$\begin{eqnarray}\tilde{a}\equiv \frac{a\ell }{R\unicode[STIX]{x1D714}_{0}^{2}}\quad \text{or}\quad \tilde{a}\equiv \frac{ak}{\unicode[STIX]{x1D714}_{0}^{2}}.\end{eqnarray}$$
This non-dimensionalization causes the three curves in each of figure 4(g,j) to collapse. For small viscosities,
$\tilde{a}$
increases linearly with
$\tilde{\unicode[STIX]{x1D708}}$
; for this reason we plot
$$\begin{eqnarray}\frac{\tilde{a}}{\tilde{\unicode[STIX]{x1D708}}}=\frac{a\ell }{R\unicode[STIX]{x1D714}_{0}^{2}}\left(\frac{\unicode[STIX]{x1D708}\ell (\ell +1)}{R^{2}\unicode[STIX]{x1D714}_{0}}\right)^{-1}=\frac{aR}{\unicode[STIX]{x1D708}\unicode[STIX]{x1D714}_{0}(\ell +1)}\quad \text{or}\quad \frac{\tilde{a}}{\tilde{\unicode[STIX]{x1D708}}}=\frac{a}{\unicode[STIX]{x1D708}\unicode[STIX]{x1D714}_{0}k}\end{eqnarray}$$
in figure 4(h,k). The form of this curve for higher viscosities shows that
$\tilde{a}$
contains terms of higher order in
$\tilde{\unicode[STIX]{x1D708}}$
, as demonstrated by Vega, Knobloch & Martel (Reference Vega, Knobloch and Martel2001).
The viscosity dependence of
$\unicode[STIX]{x1D714}$
and
$a$
, once they are non-dimensionalized, are practically identical for the capillary and gravitational cases. The difference between these, as well as the dependence on
$\ell$
, is taken into account exclusively via
$\unicode[STIX]{x1D714}_{0}(\ell )$
, just as it is for the inviscid problem.
6 Eigenmodes
Thus far we have not discussed the spatial form of the eigenmodes on the interface, beyond stating that they are spherical harmonics. Visualizations of the spherical harmonics, while common, are inadequate or incomplete for depicting the behaviour of the interface in this problem for a number of reasons. First, as stated in the introduction, the linear stability problem is degenerate: the
$2\ell +1$
spherical harmonics which share the same
$\ell$
all have the same linear growth rates and threshold. Second, for
$\ell \geqslant 4$
, the patterns actually realized in experiments or numerical simulations, which are determined by the nonlinear terms, are not individual spherical harmonics, but particular combinations of them. Finally, in a Floquet problem, the motion of the interface is time dependent.
Figure 5 shows the spherical harmonics for
$\ell =4$
. Spherical harmonics can be classified as zonal (
$m=0$
, independent of
$\unicode[STIX]{x1D719}$
, nodal lines which are circles of constant latitude), sectoral (
$m=\pm \ell$
, independent of
$\unicode[STIX]{x1D703}$
, nodal lines which are circles of constant longitude) or tesseral (
$m\neq 0,\pm \ell$
, checkered). Figure 6 depicts the pattern that is realized in numerical simulations for
$\ell =4$
. The pattern oscillates between a cube and an octahedron and is a combination of
$Y_{4}^{0}$
and
$Y_{4}^{4}$
. This pattern is the result of nonlinear selection; at the linear level, many other patterns could be realized. Our companion paper is devoted to a comprehensive description of the motion of the interface and of the velocity field for
$\ell$
between 1 and 6.
Figure 5. Spherical harmonics for
$\ell =4$
.
Figure 6. Subharmonic
$\ell =4$
standing wave pattern oscillates in time between cubic and octahedral shapes.
7 Discussion
We have considered a configuration similar to the classic problem of freely oscillating liquid drops, that of a viscous drop under the influence of a time-periodic radial bulk force and of surface tension. Here, we have carried out a linear stability analysis, while in a companion paper Ebo Adou et al. (Reference Ebo Adou, Tuckerman, Shin, Chergui and Juric2016), we describe the results of a full numerical simulation. We believe both of these investigations to be the first of their kind.
Our investigation relies on a variety of mathematical and computational tools. We have solved the linear stability problem by adapting to spherical coordinates the Floquet method of Kumar & Tuckerman (Reference Kumar and Tuckerman1994). The solution method uses a spherical harmonic decomposition in the angular directions and a Floquet decomposition in time to reduce the problem to a series of radial equations, whose solutions are monomials and spherical Bessel functions. We find that the equations for the inviscid case reduce exactly to the Mathieu equation, as they do for the planar case (Benjamin & Ursell Reference Benjamin and Ursell1954), with merely a reassignment of the parameter definitions. In contrast, for the viscous case, there are additional terms specific to the spherical geometry.
The forcing parameters for which the spherical drop is unstable are organized into tongues. The effect of viscosity is to raise, to smooth and to distort the instability tongues, both with increasing
$\unicode[STIX]{x1D708}$
and also with increasing spherical wavenumber
$\ell$
. Our computations have demonstrated the appropriate scaling for the critical oscillation frequency and amplitude with viscosity, substantially reducing the large parameter space of the problem.
The nonlinear problem is fully three-dimensional and must be treated numerically. In our companion paper (Ebo Adou et al.
Reference Ebo Adou, Tuckerman, Shin, Chergui and Juric2016), we will describe and analyse the patterns corresponding to various values of
$\ell$
that we have computed by forcing a viscous drop at appropriate frequencies.
Acknowledgements
We thank S. Fauve and D. Quéré for helpful discussions.
Appendix A
A.1 Boundary conditions for the two-fluid case
We describe here the modifications necessary in order to take into account the fluid medium surrounding the drop of radius
$R$
, occupying either a finite sphere of radius
$R_{\mathit{out}}$
or an infinite domain. The inner and outer density, dynamic viscosity and poloidal fields are denoted by
$\unicode[STIX]{x1D70C}_{j}$
,
$\unicode[STIX]{x1D702}_{j}$
and
$f^{(j)}$
for
$j=1,2$
.
$\unicode[STIX]{x0394}\unicode[STIX]{x1D6F9}\equiv [\unicode[STIX]{x1D6F9}^{(2)}-\unicode[STIX]{x1D6F9}^{(1)}]_{r=R}$
denotes the jump of any quantity
$\unicode[STIX]{x1D6F9}$
across the interface and applies to all quantities to its right within a term. For each spherical harmonic wavenumber
$\ell$
and each Floquet mode
$n$
, the poloidal fields
$f_{n}^{(j)}$
given in (3.20) are as follows:
$$\begin{eqnarray}\displaystyle & \displaystyle f_{n}^{(1)}(r)=A_{n}^{(1)}r^{\ell +1}+B_{n}^{(1)}r^{1/2}J_{\ell +1/2}(iq_{n}^{(1)}r),\qquad & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle f_{n}^{(2)}(r)=A_{n}^{(2)}r^{\ell +1}+B_{n}^{(2)}r^{1/2}J_{\ell +1/2}(iq_{n}^{(2)}r)+C_{n}^{(2)}r^{-\ell }+D_{n}^{(2)}r^{1/2}J_{-(\ell +1/2)}(\text{i}q_{n}^{(2)}r),\qquad & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}q_{n}^{(j)}\equiv \left[\frac{\unicode[STIX]{x1D707}+\text{i}(n\unicode[STIX]{x1D714}+\unicode[STIX]{x1D6FC})}{\unicode[STIX]{x1D708}_{j}}\right]^{1/2}.\end{eqnarray}$$
The introduction of four more constants requires four additional conditions. Two of these are provided by the exterior boundary conditions. The relations
$$\begin{eqnarray}\displaystyle & \displaystyle 0=u_{r}=\frac{\ell (\ell +1)}{r^{2}}f, & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle 0=\unicode[STIX]{x1D735}_{H}\boldsymbol{\cdot }\boldsymbol{u}_{H}=\frac{1}{r^{2}}\frac{\unicode[STIX]{x2202}}{\unicode[STIX]{x2202}r}r^{2}u_{r}=\frac{1}{r^{2}}\frac{\unicode[STIX]{x2202}}{\unicode[STIX]{x2202}r}\ell (\ell +1)f & \displaystyle\end{eqnarray}$$
$f_{n}^{(2)}$
and its radial derivative must vanish at
$R_{\mathit{out}}$
: 
$$\begin{eqnarray}0=f_{n}^{(2)}(R_{out})=f_{n}^{(2)\prime }(R_{out}).\end{eqnarray}$$
If the exterior is infinite, then
$A_{n}^{(2)}=B_{n}^{(2)}=0$
; otherwise (A 4) couple the six constants. Continuity of the velocity at
$r=R$
, together with (A 3) provides the remaining two additional conditions:
$$\begin{eqnarray}0=\unicode[STIX]{x0394}f=\unicode[STIX]{x0394}f^{\prime }.\end{eqnarray}$$
The kinematic condition (3.22) remains unchanged:
$$\begin{eqnarray}(\unicode[STIX]{x1D707}+\text{i}(n\unicode[STIX]{x1D714}+\unicode[STIX]{x1D6FC}))\unicode[STIX]{x1D701}_{n}=\frac{\ell (\ell +1)}{R^{2}}f_{n}|_{r=R}\end{eqnarray}$$
since
$f$
is continuous across the interface, while the tangential stress condition (3.23) becomes
$$\begin{eqnarray}0=\unicode[STIX]{x1D6E5}\left[\unicode[STIX]{x1D702}\left(\frac{\unicode[STIX]{x2202}^{2}}{\unicode[STIX]{x2202}r^{2}}-\frac{2}{r}\frac{\unicode[STIX]{x2202}}{\unicode[STIX]{x2202}r}+\frac{\ell (\ell +1)}{r^{2}}\right)f_{n}\right].\end{eqnarray}$$
The pressure jump condition (3.25) becomes
$$\begin{eqnarray}\displaystyle 0 & = & \displaystyle \unicode[STIX]{x1D6E5}\left[\unicode[STIX]{x1D70C}(\unicode[STIX]{x1D707}+i(n\unicode[STIX]{x1D714}+\unicode[STIX]{x1D6FC}))\frac{\unicode[STIX]{x2202}}{\unicode[STIX]{x2202}r}-\unicode[STIX]{x1D702}\left(\frac{\unicode[STIX]{x2202}^{3}}{\unicode[STIX]{x2202}r^{3}}+\frac{2}{r}\frac{\unicode[STIX]{x2202}^{2}}{\unicode[STIX]{x2202}r^{2}}-\ell (\ell +1)\left(\frac{3}{r^{2}}\frac{\unicode[STIX]{x2202}}{\unicode[STIX]{x2202}r}-\frac{4}{r^{3}}\right)\right)f_{n}\right.\nonumber\\ \displaystyle & & \displaystyle +\left.\left(\unicode[STIX]{x1D70C}g+\unicode[STIX]{x1D70E}\frac{(\ell -1)(\ell +2)}{R^{2}}\right)\unicode[STIX]{x1D701}_{n}-\unicode[STIX]{x1D70C}\frac{a}{2}(\unicode[STIX]{x1D701}_{n+1}+\unicode[STIX]{x1D701}_{n-1})\right],\end{eqnarray}$$
where
$\unicode[STIX]{x1D70C}$
,
$\unicode[STIX]{x1D702}$
,
$\unicode[STIX]{x2202}^{2}f/\unicode[STIX]{x2202}r^{2}$
and
$\unicode[STIX]{x2202}^{3}f/\unicode[STIX]{x2202}r^{3}$
are all discontinuous across the interface.
A.2 Differentiation relations
We express the governing equations in terms of the constants
$A_{n}$
,
$B_{n}$
,
$C_{n}$
,
$D_{n}$
via
$$\begin{eqnarray}f_{n}(r)=A_{n}r^{\ell +1}+B_{n}r^{1/2}J_{+}+C_{n}r^{-\ell }+D_{n}r^{1/2}J_{-},\qquad\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle \frac{\unicode[STIX]{x2202}}{\unicode[STIX]{x2202}r}f_{n}(r) & = & \displaystyle A_{n}(l+1)r^{\ell }+B_{n}\left(\frac{1}{2}r^{-1/2}J_{+}+\text{i}q_{n}r^{1/2}J_{+}^{\prime }\right)\nonumber\\ \displaystyle & & \displaystyle -\,C_{n}\ell r^{-\ell -1}+D_{n}\left(\frac{1}{2}r^{-1/2}J_{-}+\text{i}q_{n}r^{1/2}J_{-}^{\prime }\right),\qquad\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle \frac{\unicode[STIX]{x2202}^{2}}{\unicode[STIX]{x2202}r^{2}}f_{n}(r) & = & \displaystyle A_{n}\ell (\ell +1)r^{\ell -1}+B_{n}\left(-\frac{1}{4}r^{-3/2}J_{+}+\text{i}q_{n}r^{-1/2}J_{+}^{\prime }-q_{n}^{2}r^{1/2}J_{+}^{\prime \prime }\right)\nonumber\\ \displaystyle & & \displaystyle +\,C_{n}\ell (\ell +1)r^{-\ell -2}+D_{n}\left(-\frac{1}{4}r^{-3/2}J_{-}+\text{i}q_{n}r^{-1/2}J_{-}^{\prime }-q_{n}^{2}r^{1/2}J_{-}^{\prime \prime }\right),\qquad\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle \frac{\unicode[STIX]{x2202}^{3}}{\unicode[STIX]{x2202}r^{3}}f_{n}(r) & = & \displaystyle A_{n}(\ell +1)\ell (\ell -1)r^{\ell -2}\nonumber\\ \displaystyle & & \displaystyle +\,B_{n}\left(\frac{3}{8}r^{-5/2}J_{+}-\frac{3}{4}\text{i}q_{n}r^{-3/2}J_{+}^{\prime }-\frac{3}{2}q_{n}^{2}r^{-1/2}J_{+}^{\prime \prime }-\text{i}q_{n}^{3}r^{1/2}J_{+}^{\prime \prime \prime }\right)\nonumber\\ \displaystyle & & \displaystyle -\,C_{n}\ell (\ell +1)(\ell +2)r^{-\ell -3}\nonumber\\ \displaystyle & & \displaystyle +\,D_{n}\left(\frac{3}{8}r^{-5/2}J_{-}-\frac{3}{4}\text{i}q_{n}r^{-3/2}J_{-}^{\prime }-\frac{3}{2}q_{n}^{2}r^{-1/2}J_{-}^{\prime \prime }-\text{i}q_{n}^{3}r^{1/2}J_{-}^{\prime \prime \prime }\right),\end{eqnarray}$$
$J_{+}$
and
$J_{-}$
denote
$J_{\ell +1/2}$
and
$J_{-(\ell +1/2)}$
, respectively, to be evaluated at
$\text{i}q_{n}^{(j)}R$
. To evaluate the derivatives of the Bessel functions, we use the recurrence relations: 
$$\begin{eqnarray}\displaystyle & \displaystyle J_{\unicode[STIX]{x1D708}}^{\prime }(z)={\textstyle \frac{1}{2}}(J_{\unicode[STIX]{x1D708}-1}(z)-J_{\unicode[STIX]{x1D708}+1}(z)), & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle J_{\unicode[STIX]{x1D708}}^{\prime \prime }(z)={\textstyle \frac{1}{4}}(J_{\unicode[STIX]{x1D708}-2}(z)-2J_{\unicode[STIX]{x1D708}}(z)+J_{\unicode[STIX]{x1D708}+2}(z)), & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle J_{\unicode[STIX]{x1D708}}^{\prime \prime \prime }(z)={\textstyle \frac{1}{8}}(J_{\unicode[STIX]{x1D708}-3}(z)-3J_{\unicode[STIX]{x1D708}-1}(z)+3J_{\unicode[STIX]{x1D708}+1}(z)-J_{\unicode[STIX]{x1D708}+3}(z)). & \displaystyle\end{eqnarray}$$
$A_{n}^{(j)}$
,
$B_{n}^{(j)}$
,
$C_{n}^{(j)}$
,
$D_{n}^{(j)}$
and
$\unicode[STIX]{x1D701}_{n}$
. For the single-fluid case, we express the three conditions (3.22), (3.23), (3.25) in terms of
$A_{n}$
,
$B_{n}$
and
$\unicode[STIX]{x1D701}_{n}$
. Omitting the pressure jump condition leads to a
$6\times 6$
(finite outer sphere),
$4\times 4$
(infinite outer sphere), or
$2\times 2$
(single-fluid) system which can be inverted to obtain values for all of the constants as multiples of
$\unicode[STIX]{x1D701}_{n}$
. The pressure jump condition is then a Floquet problem in
$\{\unicode[STIX]{x1D701}_{n}\}$
, solved as described in § 3.2.
