Stability of a passive viscous droplet in a confined active nematic liquid crystal

Orientational order in the active nematic layer surrounding the droplet is described by a nematic director field $\mathbf{q}(\mathbf{x},t)$ (with $|\mathbf{q}|=1$ in the sharply aligned limit), which satisfies the Leslie-Ericksen model Ericksen (1961, 1962); Leslie (1968); De Gennes and Prost (1995); Gao et al. (2017)

$\displaystyle\frac{\partial\mathbf{q}}{\partial t}+\mathbf{u}\cdot\mathbf{\nabla q}=(\mathbf{I}-\mathbf{qq})\cdot(\varpi\mathbf{e}-\mathbf{w})\cdot\mathbf{q}+\frac{1}{\Gamma}(\mathbf{I}-\mathbf{qq})\cdot\mathbf{h},$

(1)

where $\mathbf{u}$ is the fluid velocity in the nematic. The rate-of-strain tensor $\mathbf{e}=(\nabla\mathbf{u}+\nabla\mathbf{u}^{T})/2$ and rate-of-rotation tensor $\mathbf{w}=(\nabla\mathbf{u}-\nabla\mathbf{u}^{T})/2$ depend on the velocity gradient defined as $(\nabla\mathbf{u})_{ij}=\partial_{i}u_{j}$, and $\varpi\in[0,1]$ denotes the flow alignment parameter. The last term in Eq. (1) captures nematic elasticity, and involves the rotational viscosity $\Gamma$ and molecular field $\mathbf{h}$, which derives from a free energy density. Under the commonly used one-constant approximation De Gennes and Prost (1995), we express the free energy density $\mathcal{F}_{d}$ as

$\displaystyle\mathcal{F}_{d}=\frac{K}{2}\|\nabla\mathbf{q}\|^{2},$

(2)

from which the molecular field is obtained as

$\displaystyle\mathbf{h}=-\frac{\delta\mathcal{F}_{d}}{\delta\mathbf{q}}=K\nabla^{2}\mathbf{q}.$

(3)

Nematic elasticity thus enters Eq. (1) as a diffusive term, with the ratio $K/\Gamma$ playing the role of a diffusivity. The nematic is assumed to satisfy strong parallel anchoring conditions at both the outer domain boundary and droplet interface. Note that in two dimensions the director field is equivalently described by a single angle $\beta(\mathbf{x},t)$ such that $\mathbf{q}=[\cos\beta,\sin\beta]$; we make further use of this representation below.

Conservation of mass and momentum in the active nematic read

$\displaystyle\nabla\cdot\mathbf{u}=0,\qquad\nabla\cdot(\bm{\sigma}^{h}+\bm{\sigma}^{e}+\bm{\sigma}^{a})=\bm{0},$

(4)

where $\bm{\sigma}^{h}$, $\bm{\sigma}^{a}$ and $\bm{\sigma}^{e}$ are the hydrodynamic, active, and elastic (or Leslie-Ericksen) stress tensors, respectively:

	$\displaystyle\bm{\sigma}^{h}$	$\displaystyle=-p\,\mathbf{I}+2\mu\,\mathbf{e},$		(5)
	$\displaystyle\bm{\sigma}^{a}$	$\displaystyle=-\alpha\mathbf{qq},$		(6)
	$\displaystyle\bm{\sigma}^{e}$	$\displaystyle={\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\smash{\left(\tfrac{1}{2}K\|\nabla\mathbf{q}\|^{2}-\varpi\Gamma\mathbf{qq}\,:\,\mathbf{e}\right)\mathbf{I}-K\nabla\mathbf{q}^{T}\cdot\nabla\mathbf{q}-K\left[(\nabla^{2}\mathbf{q})\mathbf{q}+\|\nabla\mathbf{q}\|^{2}\mathbf{qq}\right]+\varpi\Gamma(\mathbf{qq}\,:\,\mathbf{e})\mathbf{qq}.}}$		(7)

Here, $p$ is the fluid pressure and $\mu$ is the shear viscosity of the fluid. The active stress magnitude is $\alpha=n\varsigma_{a}$, where $n$ is the number density of active dipoles and $\varsigma_{a}$ is the dipole strength or stresslet; $\alpha>0$ for an extensile nematic, and $\alpha<0$ in the contractile case Saintillan and Shelley (2013); Saintillan (2018).

In the passive viscous drop, the only stress is hydrodynamic, and the governing equations are

$\displaystyle\nabla\cdot\overline{\mathbf{u}}=0,\qquad\nabla\cdot\overline{\bm{\sigma}}^{h}=\bm{0},$

(8)

where $\overline{\bm{\sigma}}^{h}=-\overline{p}\,\mathbf{I}+2\overline{\mu}\,\overline{\mathbf{e}}$ and $\overline{\mathbf{e}}=\tfrac{1}{2}(\nabla\overline{\mathbf{u}}+\nabla\overline{\mathbf{u}}^{\mathrm{T}})$. Here and in the following, variables associated with flow inside the drop are decorated with an overbar.

At the outer boundary $r=R$, the no-slip condition $\mathbf{u}=\mathbf{0}$ applies. There, the nematic director is tangent to the boundary: $\mathbf{q}\cdot\hat{\mathbf{r}}=0$.

To parametrize the shape of the drop, we introduce the function $F(r,\theta,t)=r-\eta R-\zeta(\theta,t)$ and note that the interface $\mathcal{C}$ between the two fluids is defined as $F=0$. Consequently, the outward unit normal to the interface is

$\displaystyle\mathbf{n}=\frac{\nabla F}{|\nabla F|}.$

(9)

The kinematic boundary condition describing the evolution of the interface is

$\displaystyle\frac{\partial F}{\partial t}+\mathbf{u}\cdot\nabla F=0,$

(10)

i.e.,

$\displaystyle\frac{\partial\zeta}{\partial t}=u_{r}-\frac{u_{\theta}}{r}\frac{\partial\zeta}{\partial\theta}.$

(11)

Note that the fluid velocity is continuous at the interface: $\mathbf{u}=\overline{\mathbf{u}}$ at $F=0$.

Neglecting Marangoni effects, the dynamic boundary condition is given by the stress balance

$$\mathbf{n}\cdot(\bm{\sigma}^{h}+\bm{\sigma}^{a}+\bm{\sigma}^{e}-\overline{\bm{\sigma}}^{h})=\gamma\kappa\mathbf{n},$$

(12)

where $\gamma$ is the surface tension, assumed to be uniform, and $\kappa=\nabla\cdot\mathbf{n}$ is the local curvature. The parallel anchoring condition for the nematic at the drop surface reads: $\mathbf{q}\cdot\mathbf{n}=0$ at $F=0$.

We scale the equations using length scale $R$, time scale $\mu/|\alpha|$, and velocity scale $|\alpha|R/\mu$. The nematodynamic equation (1) becomes

$\displaystyle\frac{\partial\mathbf{q}}{\partial t}+\mathbf{u}\cdot\mathbf{\nabla q}=(\mathbf{I}-\mathbf{qq})\cdot(\varpi\mathbf{e}-\mathbf{w})\cdot\mathbf{q}+\frac{1}{\mathrm{Pe}}(\mathbf{I}-\mathbf{qq})\cdot\nabla^{2}\mathbf{q},$

(13)

where the active Péclet number $\mathrm{Pe}=\Gamma|\alpha|R^{2}/K\mu$ compares the effects of the fluid flow to nematic relaxation on the dynamics of $\mathbf{q}$.

In the rest of the paper, we will assume $\mathrm{Pe}\to 0$ (strong nematic relaxation), in which case the equation for $\mathbf{q}$ simplifies to

$\displaystyle(\mathbf{I}-\mathbf{q}\mathbf{q})\cdot\nabla^{2}\mathbf{q}=\mathbf{0}\qquad\implies\qquad\nabla^{2}\beta=0.$

(14)

In this limit, the configuration of the nematic is thus purely a boundary value problem governed by the geometry of the nematic layer through the anchoring boundary conditions. The only effect of the flow on $\mathbf{q}$ is through the kinematic boundary condition (11), which governs the shape of the interface.

In the nematic, we use $|\alpha|$ to scale the pressure and stress. The stress tensors inside the active nematic then become

$$\bm{\sigma}^{h}=-p\,\mathbf{I}+2\,\mathbf{e},\qquad\bm{\sigma}^{a}=-S\mathbf{qq},\qquad\bm{\sigma}^{e}=-\xi^{2}\nabla\mathbf{q}^{\mathrm{T}}\cdot\nabla\mathbf{q},$$

(15)

Note that the elastic stress was simplified under that assumption of $\mathrm{Pe}\rightarrow 0$ and by making use of Eq. (14), and that its isotropic part is henceforth omitted as it can absorbed in the fluid pressure. Inside the passive drop, the hydrodynamic stress reads

$$\overline{\bm{\sigma}}^{h}=-\overline{p}\,\mathbf{I}+2\lambda\,\overline{\mathbf{e}}.$$

(16)

where $\lambda=\overline{\mu}/\mu$ is the viscosity ratio. The active stress coefficient in Eq. (15) is $S=\alpha/|\alpha|$, with $S=+1$ for an extensile nematic and $S=-1$ for a contractile nematic. The relative magnitude of elastic over active stresses is captured by the dimensionless parameter $\xi=\sqrt{K/|\alpha|}/R$, or ratio of the active length scale $\sqrt{K/|\alpha|}$ to the system size. Active stresses are negligible when $\xi\gg 1$, whereas they dominate when $\xi\ll 1$, a regime where active turbulence is typically observed. In the present work, we assume $\xi\sim O(1)$.

The kinematic boundary condition (11) remains unchanged under non-dimensionalization, and the stress balance (12) becomes

$$\mathbf{n}\cdot(\bm{\sigma}^{h}+\bm{\sigma}^{a}+\bm{\sigma}^{e}-\overline{\bm{\sigma}}^{h})=\frac{1}{\mathrm{Ca}}\kappa\mathbf{n},$$

(17)

where we have introduced the active capillary number $\mathrm{Ca}=|\alpha|R/\gamma$ comparing the magnitudes of active vs. capillary stresses. Under the assumption of vanishing Péclet number, the system is entirely governed by five parameters: the radius ratio $\eta$, the sign $S$ of active stresses, the capillary number $\mathrm{Ca}$, the dimensionless elastic constant $\xi$, and the viscosity ratio $\lambda$.

In the unperturbed state where the drop is circular and centered at the origin ($\zeta_{0}=0$), the nematic field lines are concentric circles with

$$\beta_{0}=\arctan\frac{y}{x}+\frac{\pi}{2}=\theta+\frac{\pi}{2},$$

(18)

which satisfies Eq. (14) along with the anchoring boundary conditions at $r=\eta$ and $1$. The corresponding nematic director is $\mathbf{q}_{0}=\hat{\bm{\theta}}$. Both fluids are still ($\mathbf{u}_{0}=\overline{\mathbf{u}}_{0}=\mathbf{0}$). The base-state active and elastic stress tensors in the outer nematic layer are given by

$$\bm{\sigma}^{a}_{0}=-S\mathbf{q}_{0}\mathbf{q}_{0}=-S\hat{\bm{\theta}}\hat{\bm{\theta}},\qquad\bm{\sigma}^{e}_{0}=-\xi^{2}\nabla\mathbf{q}_{0}^{\mathrm{T}}\cdot\nabla\mathbf{q}_{0}=-\frac{\xi^{2}}{r^{2}}\hat{\mathbf{r}}\hat{\mathbf{r}}.$$

(19)

They both induce a pressure gradient in the nematic fluid:

$$\nabla p_{0}=\nabla\cdot(\bm{\sigma}^{a}_{0}+\bm{\sigma}^{e}_{0})=\left(\frac{S}{r}+\frac{\xi^{2}}{r^{3}}\right)\hat{\mathbf{r}}=\nabla\left(S\ln r-\frac{\xi^{2}}{2r^{2}}\right),$$

(20)

from which

$$p_{0}(r)=S\ln r-\frac{\xi^{2}}{2r^{2}}.$$

(21)

The pressure inside the passive drop is uniform and is dictated by the stress balance (17):

$$\overline{p}_{0}=S\ln\eta+\frac{\xi^{2}}{2\eta^{2}}+\frac{1}{\eta\mathrm{Ca}}.$$

(22)

Leveraging periodicity in the $\theta$ direction, we perturb the shape of the interface with an azimuthal Fourier mode

$$\zeta(\theta,t)=\epsilon\zeta_{1}(\theta,t)=\epsilon\widetilde{\zeta}_{1}(t)\cos m\theta,$$

(23)

where $|\epsilon|\ll 1$ and $m\in\mathbb{N}$. We assume $m$ to be given, and our goal is to analyze the growth of the perturbation amplitude $\widetilde{\zeta}_{1}(t)$. At linear order, the perturbed outward normal and curvature along the interface are, respectively,

	$\displaystyle\mathbf{n}$	$\displaystyle=\mathbf{n}_{0}+\epsilon\mathbf{n}_{1}+O(\epsilon^{2})=\hat{\mathbf{r}}+\epsilon\widetilde{\zeta}_{1}\frac{m}{\eta}\sin m\theta\,\hat{\bm{\theta}}+O(\epsilon^{2}),$		(24)
	$\displaystyle\kappa$	$\displaystyle=\kappa_{0}+\epsilon\kappa_{1}+O(\epsilon^{2})=\frac{1}{\eta}+\epsilon\widetilde{\zeta}_{1}\frac{m^{2}-1}{\eta^{2}}\cos m\theta+O(\epsilon^{2}).$		(25)

The instantaneous shape of the interface fully governs that of the nematic field lines. The director angle is expanded as $\beta(r,\theta)=\beta_{0}+\epsilon\beta_{1}+O(\epsilon^{2})$. The boundary conditions on the perturbation angle are $\beta_{1}=0$ at $r=1$, and

$$\beta_{1}=\widetilde{\zeta}_{1}\frac{m}{\eta}\sin m\theta\quad\mathrm{at}\quad r=\eta.$$

(26)

A solution of Laplace’s equation $\nabla^{2}\beta_{1}=0$ that satisfies both of these conditions is readily obtained by separation of variables,

$$\beta_{1}(r,\theta)=\widetilde{\zeta}_{1}\frac{m\eta^{m-1}}{\eta^{2m}-1}(r^{m}-r^{-m})\sin m\theta,$$

(27)

and the corresponding expansion for the nematic director follows as

$$\mathbf{q}=\hat{\bm{\theta}}-\epsilon\beta_{1}\hat{\mathbf{r}}+O(\epsilon^{2}).$$

(28)

Typical director fields corresponding to $m=1$ and $5$ are illustrated in Fig. 1(b,c).

The linearization of active and elastic stresses is obtained by inserting Eq. (28) into Eq. (15). Starting with the active stress,

$$\bm{\sigma}^{a}=\bm{\sigma}^{a}_{0}+\epsilon\bm{\sigma}^{a}_{1}+O(\epsilon^{2})=-S\hat{\bm{\theta}}\hat{\bm{\theta}}+\epsilon S\beta_{1}\big{(}\hat{\mathbf{r}}\hat{\bm{\theta}}+\hat{\bm{\theta}}\hat{\mathbf{r}}\big{)}+O(\epsilon^{2}).$$

(29)

To obtain the elastic stress, we first calculate the gradient of $\mathbf{q}$,

$$\nabla\mathbf{q}=\nabla\hat{\bm{\theta}}-\epsilon\left[(\nabla\beta_{1})\hat{\mathbf{r}}+\beta_{1}\nabla\hat{\mathbf{r}}\right]+O(\epsilon^{2})=-\frac{1}{r}\hat{\bm{\theta}}\hat{\mathbf{r}}-\epsilon\left[\partial_{r}\beta_{1}\hat{\mathbf{r}}\hat{\mathbf{r}}+\frac{\partial_{\theta}\beta_{1}}{r}\hat{\bm{\theta}}\hat{\mathbf{r}}+\frac{\beta_{1}}{r}\hat{\bm{\theta}}\hat{\bm{\theta}}\right]+O(\epsilon^{2}).$$

(30)

The linearized elastic stress then reads

$$\bm{\sigma}^{e}=\bm{\sigma}^{e}_{0}+\epsilon\bm{\sigma}^{e}_{1}+O(\epsilon^{2})=-\frac{\xi^{2}}{r^{2}}\hat{\mathbf{r}}\hat{\mathbf{r}}-\epsilon\frac{\xi^{2}}{r^{2}}\left[2\partial_{\theta}\beta_{1}\hat{\mathbf{r}}\hat{\mathbf{r}}+\beta_{1}\big{(}\hat{\mathbf{r}}\hat{\bm{\theta}}+\hat{\bm{\theta}}\hat{\mathbf{r}}\big{)}\right]+O(\epsilon^{2}).$$

(31)

The fluid velocity and pressure perturbations satisfy the Stokes equations (4) and (8) at linear order in $\epsilon$:

	$\displaystyle\nabla\cdot{\mathbf{u}}_{1}=0,\qquad\nabla\cdot(\bm{\sigma}^{h}_{1}+\bm{\sigma}^{a}_{1}+\bm{\sigma}^{e}_{1})=\mathbf{0},$	$\displaystyle\mathrm{for}\,\,\,r>\eta+\zeta,$		(32)
	$\displaystyle\nabla\cdot\overline{\mathbf{u}}_{1}=0,\qquad\nabla\cdot\overline{\bm{\sigma}}^{h}_{1}=\mathbf{0},$	$\displaystyle\mathrm{for}\,\,\,r<\eta+\zeta,$		(33)

where the active and elastic stresses in Eq. (32) are known from Eqs. (29) and (31). The flows in the two regions are coupled via the dynamic boundary condition (17), which applies at $r=\eta+\zeta$.

Our method of solution, described in Sec. III.5, will rely on the Lorentz reciprocal theorem for Stokes flow Masoud and Stone (2019). To this end, it is convenient to first recast the flow problem around a circular boundary ($r=\eta$). Note that Eqs. (32)–(33) remain unchanged around the circle to linear order. We also need to apply the dynamic boundary condition (17) on the circle, which is achieved by expanding the various fields involved in the boundary conditions as Taylor series in the neighborhood of $\zeta=0$. For example, the stress tensor can be expanded as:

$$\bm{\sigma}(r=\eta+\epsilon\zeta_{1})=\bm{\sigma}(r=\eta)+\epsilon\zeta_{1}\partial_{r}\bm{\sigma}(r=\eta)+O(\epsilon^{2}).$$

(34)

Inserting Eq. (34) into Eq. (17), linearizing, and simplifying using Eqs. (19)–(22), we obtain the following dynamic boundary condition at $O(\epsilon)$,

$$\hat{\mathbf{r}}\cdot\big{(}\bm{\sigma}^{h}_{1}+\bm{\sigma}^{a}_{1}+\bm{\sigma}^{e}_{1}-\overline{\bm{\sigma}}^{h}_{1}\big{)}=\frac{\kappa_{1}\hat{\mathbf{r}}}{\mathrm{Ca}}+\frac{\widetilde{\zeta_{1}}}{\eta}\left(S-\frac{\xi^{2}}{\eta^{2}}\right)\big{(}\cos m\theta\,\hat{\mathbf{r}}+m\sin m\theta\,\hat{\bm{\theta}}\big{)}\qquad\mathrm{at}\,\,\,r=\eta.$$

(35)

Finally, the linearized kinematic boundary condition (11) reads

$$\frac{\partial\zeta_{1}}{\partial t}=\hat{\mathbf{r}}\cdot\mathbf{u}_{1}\qquad\mathrm{at}\,\,\,r=\eta.$$

(36)

In this linearized problem, Eq. (36) is an amplitude evolution equation on the inner circle, driven by the solution of the two boundary value problems inside the inner circle and the outer annulus. This reduction allows the use of the Lorentz reciprocal theorem in the circular geometry as we describe next.

We apply the Lorentz reciprocal theorem Masoud and Stone (2019) to the linear problem of Eqs. (32)–(36) in the annular geometry with a circular boundary. To this end, we introduce an auxiliary flow problem in the same annular geometry, in which both fluids are passive and the flow is driven by a prescribed stress jump at the interface. In this problem, both fluids are Newtonian with viscosity ratio $\lambda=\overline{\mu}/\mu$ and the interface between them is circular and centered at the origin ($\zeta=0$). We use uppercase fonts ($\mathbf{U}$, $P$, $\bm{\Sigma}^{h}$,…) for all variables associated with the auxiliary problem. The flow in both fluid regions is governed by the homogeneous Stokes equations,

	$\displaystyle\nabla\cdot{\mathbf{U}}=0,\qquad\nabla\cdot\bm{\Sigma}^{h}=\mathbf{0},$	$\displaystyle\mathrm{for}\,\,\,r>\eta,$		(37)
	$\displaystyle\nabla\cdot\overline{\mathbf{U}}=0,\qquad\nabla\cdot\overline{\bm{\Sigma}}^{h}=\mathbf{0},$	$\displaystyle\mathrm{for}\,\,\,r<\eta.$		(38)

Here, $\bm{\Sigma}^{h}=-P\mathbf{I}+2\,\mathbf{E}$ and $\overline{\bm{\Sigma}}^{h}=-\overline{P}\mathbf{I}+2\lambda\,\overline{\mathbf{E}}$ are the Newtonian stress tensors, where $\mathbf{E}$ and $\overline{\mathbf{E}}$ are the respective rate-of-strain tensors. The velocity satisfies the no-slip condition at the outer boundary: $\mathbf{U}(r=1)=\mathbf{0}$, and is continuous at the interface: $\mathbf{U}(r=\eta)=\overline{\mathbf{U}}(r=\eta)$. The flow is driven by a prescribed jump in normal tractions at the interface:

$$\hat{\mathbf{r}}\cdot\big{(}\bm{\Sigma}^{h}-\overline{\bm{\Sigma}}^{h}\big{)}=\Sigma_{0}\cos m\theta\,\hat{\mathbf{r}},$$

(39)

where $\Sigma_{0}$ is an arbitrary constant. In two dimensions, this Newtonian problem is readily solved in terms of two streamfunctions $\psi(r,\theta)$ and $\overline{\psi}(r,\theta)$ such that

$$\psi(r,\theta)=\Sigma_{0}G_{m}(r)\sin m\theta,\qquad U_{r}=\frac{1}{r}\frac{\partial\psi}{\partial\theta}=\Sigma_{0}\frac{G_{m}(r)}{r}m\cos m\theta,\qquad U_{\theta}=-\frac{\partial\psi}{\partial r}=-\Sigma_{0}G_{m}^{\prime}(r)\sin m\theta,$$

(40)

with similar expressions for the streamfunction and its derivatives inside the drop. The two dimensionless functions $G_{m}(r)$ and $\overline{G}_{m}(r)$ can be obtained analytically by the solution of the Stokes equations; details of their solutions are presented in the Appendix.

Using Eqs. (32) and (37) in the outer nematic fluid region $\Omega_{0}^{+}=\{(r,\theta)\,|\,r>\eta\}$ and relying on the symmetry of the stress tensors, it is straightforward to derive the two relations

	$\displaystyle\nabla\cdot\left[(\bm{\sigma}^{h}_{1}+\bm{\sigma}^{a}_{1}+\bm{\sigma}^{e}_{1})\cdot\mathbf{U}\right]$	$\displaystyle=(\bm{\sigma}^{h}_{1}+\bm{\sigma}^{a}_{1}+\bm{\sigma}^{e}_{1}):\mathbf{E},$		(41)
	$\displaystyle\nabla\cdot\left[\bm{\Sigma}^{h}\cdot\mathbf{u}_{1}\right]$	$\displaystyle=\bm{\Sigma}^{h}:\mathbf{e}_{1}.$		(42)

Subtracting Eq. (42) from Eq. (41), integrating over the outer region $\Omega_{0}^{+}$, and applying the divergence theorem then yields

$$\int_{\mathcal{C}_{0}}\hat{\mathbf{r}}\cdot\left[\bm{\Sigma}^{h}\cdot\mathbf{u}_{1}-(\bm{\sigma}^{h}_{1}+\bm{\sigma}^{a}_{1}+\bm{\sigma}^{e}_{1})\cdot\mathbf{U}\right]\mathrm{d}\ell=\int_{\Omega_{0}^{+}}(\bm{\sigma}^{a}_{1}+\bm{\sigma}^{e}_{1}):\mathbf{E}\,\mathrm{d}S,$$

(43)

where $\mathcal{C}_{0}$ denotes the circular interface at $r=\eta$. Note that the boundary term at $r=1$ vanished due to the no-slip condition, and that we used that $\bm{\sigma}^{h}_{1}:\mathbf{E}=\bm{\Sigma}^{h}:\mathbf{e}_{1}$ to simplify the right-hand side. Following an identical process in the inner region $\Omega_{0}^{-}=\{(r,\theta)\,|\,r<\eta\}$, we obtain,

$$\int_{\mathcal{C}_{0}}\hat{\mathbf{r}}\cdot\left[\overline{\bm{\Sigma}}^{h}\cdot\mathbf{u}_{1}-\overline{\bm{\sigma}}^{h}_{1}\cdot\mathbf{U}\right]\mathrm{d}\ell=0,$$

(44)

which does not involve any bulk integral since active and elastic stresses are absent inside $\Omega_{0}^{-}$. Combining Eqs. (43) and (44), we obtain

$$\int_{\mathcal{C}_{0}}\hat{\mathbf{r}}\cdot\left[\big{(}\bm{\Sigma}^{h}-\overline{\bm{\Sigma}}^{h}\big{)}\cdot\mathbf{u}_{1}-\big{(}\bm{\sigma}^{h}_{1}+\bm{\sigma}^{e}_{1}+\bm{\sigma}^{a}_{1}-\overline{\bm{\sigma}}^{h}_{1}\big{)}\cdot\mathbf{U}\right]\mathrm{d}\ell=\int_{\Omega_{0}^{+}}\big{(}\bm{\sigma}^{e}_{1}+\bm{\sigma}^{a}_{1}\big{)}:\mathbf{E}\,\mathrm{d}S,$$

(45)

which is the central statement of the reciprocal theorem. The first contribution to the contour integral on the left-hand side can be simplified using the linearized kinematic boundary condition (36) along with the auxiliary stress balance (39), yielding

$$\int_{\mathcal{C}_{0}}\hat{\mathbf{r}}\cdot\big{(}\bm{\Sigma}^{h}-\overline{\bm{\Sigma}}^{h}\big{)}\cdot\mathbf{u}_{1}\,\mathrm{d}\ell=\int_{0}^{2\pi}\Sigma_{0}\cos m\theta\,\left(\hat{\mathbf{r}}\cdot\mathbf{u}_{1}\right)\,\eta\,\mathrm{d}\theta=\pi\eta\Sigma_{0}\frac{\mathrm{d}\widetilde{\zeta}_{1}}{\mathrm{d}t}.$$

(46)

Similarly, the second term in the contour integral can be simplified using the linearized dynamic boundary condition (35) along with Eq. (40) for the auxiliary velocity:

	$\displaystyle\int_{\mathcal{C}_{0}}\hat{\mathbf{r}}\cdot\big{(}\bm{\sigma}^{h}_{1}+\bm{\sigma}^{e}_{1}+\bm{\sigma}^{a}_{1}-\overline{\bm{\sigma}}^{h}_{1}\big{)}\cdot\mathbf{U}\,\mathrm{d}\ell$	$\displaystyle=\int_{0}^{2\pi}\bigg{[}\frac{\kappa_{1}U_{r}}{\mathrm{Ca}}+\frac{\widetilde{\zeta}_{1}}{\eta}\left(S-\frac{\xi^{2}}{\eta^{2}}\right)\big{(}U_{r}\cos m\theta+U_{\theta}m\sin m\theta\big{)}\bigg{]}\eta\,\mathrm{d}\theta$		(47)
		$\displaystyle=m\pi\Sigma_{0}\widetilde{\zeta}_{1}\left[\frac{m^{2}-1}{\mathrm{Ca}}\frac{G_{m}(\eta)}{\eta^{2}}+\left(S-\frac{\xi^{2}}{\eta^{2}}\right)\left(\frac{G_{m}(\eta)}{\eta}-G_{m}^{\prime}(\eta)\right)\right].$		(48)

Finally, the surface integral on the right-hand side of Eq. (45) can be evaluated using Eqs. (29) and (31) for the linearized active and elastic stresses:

	$\displaystyle\int_{\Omega_{0}^{+}}\big{(}\bm{\sigma}^{e}_{1}+\bm{\sigma}^{a}_{1}\big{)}:\mathbf{E}\,\mathrm{d}S$	$\displaystyle=\int_{0}^{2\pi}\int_{\eta}^{1}\bigg{[}2S\beta_{1}E_{r\theta}-\frac{2\xi^{2}}{r^{2}}\big{(}\partial_{\theta}\beta_{1}E_{rr}+\beta_{1}E_{r\theta}\big{)}\bigg{]}r\mathrm{d}r\mathrm{d}\theta$		(49)
	$\displaystyle\begin{split}&=-\pi\Sigma_{0}\widetilde{\zeta}_{1}\frac{m\eta^{m-1}}{\eta^{2m}-1}\left\{S\int_{\eta}^{1}\big{(}r^{m}-r^{-m}\big{)}\bigg{[}rG_{m}^{\prime\prime}(r)-G_{m}^{\prime}(r)+m^{2}\frac{G_{m}(r)}{r}\bigg{]}\mathrm{d}r\right.\\ &\,\,\,\,\,\,-\xi^{2}\left.\int_{\eta}^{1}\big{(}r^{m}-r^{-m}\big{)}\bigg{[}\frac{G_{m}^{\prime\prime}(r)}{r}-(2m^{2}+1)\frac{G_{m}^{\prime}(r)}{r^{2}}+3m^{2}\frac{G_{m}(r)}{r^{3}}\bigg{]}\mathrm{d}r\right\}.\end{split}$			(50)

Upon inserting Eqs. (46), (48) and (50) into the reciprocal relation (45), the arbitrary stress jump magnitude $\Sigma_{0}$ can be eliminated, and we arrive at a first-order ODE for the interfacial perturbation amplitude,

$$\frac{\mathrm{d}\widetilde{\zeta}_{1}}{\mathrm{d}t}=\Lambda_{m}\widetilde{\zeta}_{1}.$$

(51)

The constant growth rate $\Lambda_{m}$ for normal mode $m$ includes contributions from capillary, active and elastic stresses:

$$\Lambda_{m}=\frac{1}{\mathrm{Ca}}\Xi^{c}_{m}(\eta,\lambda)+S\,\Xi^{a}_{m}(\eta,\lambda)+\xi^{2}\,\Xi^{e}_{m}(\eta,\lambda),$$

(52)

respectively given by

	$\displaystyle\Xi^{c}_{m}(\eta,\lambda)$	$\displaystyle=m(m^{2}-1)\frac{G_{m}(\eta)}{\eta^{3}},$		(53)
	$\displaystyle\Xi^{a}_{m}(\eta,\lambda)$	$\displaystyle={\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}-m\frac{\mathrm{d}}{\mathrm{d}\eta}\bigg{[}\frac{G_{m}(\eta)}{\eta}\bigg{]}}-\frac{m\eta^{m-2}}{\eta^{2m}-1}\int_{\eta}^{1}\big{(}r^{m}-r^{-m}\big{)}\bigg{[}rG_{m}^{\prime\prime}(r)-G_{m}^{\prime}(r)+m^{2}\frac{G_{m}(r)}{r}\bigg{]}\mathrm{d}r,$		(54)
	$\displaystyle\Xi^{e}_{m}(\eta,\lambda)$	$\displaystyle={\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\frac{m}{\eta^{2}}\frac{\mathrm{d}}{\mathrm{d}\eta}\bigg{[}\frac{G_{m}(\eta)}{\eta}\bigg{]}}+\frac{m\eta^{m-2}}{\eta^{2m}-1}\int_{\eta}^{1}\big{(}r^{m}-r^{-m}\big{)}\bigg{[}\frac{G_{m}^{\prime\prime}(r)}{r}-(2m^{2}+1)\frac{G_{m}^{\prime}(r)}{r^{2}}+3m^{2}\frac{G_{m}(r)}{r^{3}}\bigg{]}\mathrm{d}r.$		(55)

These various growth rates depend on the function $G_{m}(r)$ appearing in the solution of the Newtonian auxiliary problem, whose form differs between the translation ($m=1$) and deformation ($m\geq 2$) modes of the drop. We analyze and discuss these two cases separately.

Governing equations (1)–(8) still apply, the only difference being that the active and elastic stresses of Eqs. (6) and (7) are now applied inside the drop. As a result, they are decorated by an overbar and enter the momentum equation (8) instead of (4). The kinematic boundary condition (11) is unchanged, while the dynamic boundary condition, in dimensionless form, becomes:

$$\mathbf{n}\cdot(\bm{\sigma}^{h}-\overline{\bm{\sigma}}^{a}-\overline{\bm{\sigma}}^{e}-\overline{\bm{\sigma}}^{h})=\frac{1}{\mathrm{Ca}}\kappa\mathbf{n}.$$

(66)