2.1 Governing equations and boundary conditions

A liquid droplet of a Newtonian fluid is considered here, with density $\unicode[STIX]{x1D70C}_{i}$ and dynamic viscosity $\unicode[STIX]{x1D702}_{i}$ , submerged in a second Newtonian fluid of density $\unicode[STIX]{x1D70C}_{o}$ and viscosity $\unicode[STIX]{x1D702}_{o}$ , containing a surfactant solute of concentration $C$ , as sketched in figure 1. Note that different experimental set-ups employ different inner and outer fluids: for instance, Izri et al. (Reference Izri, van der Linden, Michelin and Dauchot2014) observed water droplets in oil, while others used oil droplets in water (Moerman et al. Reference Moerman, Moyses, van der Wee, Grie, van Blaaderen, Kegel, Groenewold and Brujic2017). To remain general, we denote in the following by subscripts $i$ and $o$ the quantities relevant to the inner and outer fluid, respectively. Far from the droplet, the concentration of surfactant molecules is ${\mathcal{C}}_{\infty }$ . In the following, we consider axisymmetric deformations of the droplet under flow and Marangoni stresses, whose surface is thus described in spherical polar coordinates by $r=R(t,\unicode[STIX]{x1D707})$ with $\unicode[STIX]{x1D707}=\cos \unicode[STIX]{x1D703}$ . Denoting the radius of the drop at rest (i.e. when it is spherical) as $R_{0}$ , the radius of the deformed droplet can be written as

(2.1) $$\begin{eqnarray}\displaystyle R(t,\unicode[STIX]{x1D707})=R_{0}(1+\unicode[STIX]{x1D709}(t,\unicode[STIX]{x1D707})). & & \displaystyle\end{eqnarray}$$

Naturally, the presence of deformations would contribute not only to the position of the droplet interface, but also to its curvature. That is, vectors normal and tangential to the interface are now functions of $t$ and $\unicode[STIX]{x1D707}$ , thus affecting the boundary conditions formulated below.

Figure 1.

Cross-section of an axisymmetric deformable active droplet undergoing gradual micellar dissolution. As the droplet dissolves, it adsorbs surfactant that powers the chemical reaction at the droplet interface, thus sustaining the dissolution.

Multiple physico-chemical mechanisms have been identified in experiments leading to the self-propulsion of active droplets (Herminghaus et al. Reference Herminghaus, Maass, Krüger, Thutupalli, Goehring and Bahr2014). In the following, we explicitly refer to the molecular pathway identified in the experiments of Moerman et al. (Reference Moerman, Moyses, van der Wee, Grie, van Blaaderen, Kegel, Groenewold and Brujic2017). Yet the formalism presented here is completely general and could easily be applied to the micellar pathway relevant to the experiments of Izri et al. (Reference Izri, van der Linden, Michelin and Dauchot2014). In the molecular pathway leading to the solubilization of the oil phase into the aqueous solution, surfactant molecules are absorbed at the surface and swollen micelles are released, leading to slow decrease of the droplet size. In experiments, the typical time of droplet dissolution is substantially longer than the time scale associated with self-propulsion (Izri et al. Reference Izri, van der Linden, Michelin and Dauchot2014; Moerman et al. Reference Moerman, Moyses, van der Wee, Grie, van Blaaderen, Kegel, Groenewold and Brujic2017), so that this dissolution process can be neglected and the volume of the droplet is assumed constant, $V=\text{const.}$ , while the droplet consumes surfactant molecules at a fixed rate ${\mathcal{A}}>0$ ,

(2.2) $$\begin{eqnarray}\displaystyle {\mathcal{D}}\boldsymbol{n}\boldsymbol{\cdot }\unicode[STIX]{x1D735}C={\mathcal{A}}\quad \text{at }r=R, & & \displaystyle\end{eqnarray}$$

where $\boldsymbol{n}$ is the outward normal to the droplet interface. Surfactant molecules do not penetrate into the droplet. Thus, advection–diffusion of surfactant should only be taken into account outside of the drop,

(2.3) $$\begin{eqnarray}\displaystyle \unicode[STIX]{x2202}_{t}C+\boldsymbol{u}_{o}\boldsymbol{\cdot }\unicode[STIX]{x1D735}C={\mathcal{D}}\unicode[STIX]{x1D6FB}^{2}C, & & \displaystyle\end{eqnarray}$$

where $\unicode[STIX]{x2202}_{t}$ denotes the partial derivative with respect to time $t$ , $\boldsymbol{u}_{o}$ is the flow velocity outside of the drop, and ${\mathcal{D}}$ denotes the molecular diffusivity of the surfactant in the outer fluid. Naturally, far away from the droplet, surfactant concentration reaches a constant value,

(2.4) $$\begin{eqnarray}\displaystyle C\rightarrow {\mathcal{C}}_{\infty }\quad \text{for }r\rightarrow \infty . & & \displaystyle\end{eqnarray}$$

The presence of surfactant at the droplet’s surface modifies its interfacial tension $\unicode[STIX]{x1D6FE}$ . Assuming that adsorption/desorption of surfactant molecules at the fluid–fluid interface occurs instantaneously compared to its transport in the outer fluid, the concentration of the adsorbed surfactants is $\propto C|_{r=R}$ (Baret Reference Baret1969). We further linearize the relationship between $\unicode[STIX]{x1D6FE}$ and $C|_{r=R}$ ,

(2.5) $$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D6FE}=\unicode[STIX]{x1D6FE}_{0}-\unicode[STIX]{x1D6FE}_{C}(C|_{r=R}-{\mathcal{C}}_{\infty }+{\mathcal{A}}R_{0}/{\mathcal{D}}), & & \displaystyle\end{eqnarray}$$

and note that $\unicode[STIX]{x1D6FE}_{C}\equiv -(\text{d}\unicode[STIX]{x1D6FE}/\text{d}C)|_{C={\mathcal{C}}_{\infty }-{\mathcal{A}}R_{0}/{\mathcal{D}}}>0$ with ${\mathcal{C}}_{\infty }-{\mathcal{A}}R_{0}/{\mathcal{D}}$ corresponding to the surfactant concentration at $r=R$ in the absence of flow. The stress balance at the interface can then be written in vector form (i.e. accounting for both normal and tangential stresses) as

(2.6) $$\begin{eqnarray}\displaystyle (\unicode[STIX]{x1D748}_{o}-\unicode[STIX]{x1D748}_{i})\boldsymbol{\cdot }\boldsymbol{n}+\unicode[STIX]{x1D735}\boldsymbol{\cdot }[(\unicode[STIX]{x1D644}-\boldsymbol{n}\boldsymbol{n})\unicode[STIX]{x1D6FE}]=0\quad \text{at }r=R, & & \displaystyle\end{eqnarray}$$

where $\unicode[STIX]{x1D644}$ is the identity tensor, $\unicode[STIX]{x1D748}_{i,o}=-P_{i,o}\unicode[STIX]{x1D644}+\unicode[STIX]{x1D749}_{i,o}$ is the hydrodynamic stress tensor with $\unicode[STIX]{x1D749}_{i,o}=\unicode[STIX]{x1D702}_{i,o}(\unicode[STIX]{x1D735}\boldsymbol{u}_{i,o}+\unicode[STIX]{x1D735}\boldsymbol{u}_{i,o}^{\text{T}})$ its viscous part. The continuity of the fluid’s velocity and impermeability of the droplet’s surface is written as

(2.7) $$\begin{eqnarray}\displaystyle R_{0}\frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D709}}{\unicode[STIX]{x2202}t}-R_{0}\frac{\sqrt{1-\unicode[STIX]{x1D707}^{2}}}{r}\frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D709}}{\unicode[STIX]{x2202}\unicode[STIX]{x1D707}}(\boldsymbol{u}_{o}\boldsymbol{\cdot }\boldsymbol{e}_{\unicode[STIX]{x1D703}})=\boldsymbol{u}_{o}\boldsymbol{\cdot }\boldsymbol{e}_{r},\quad \boldsymbol{u}_{o}=\boldsymbol{u}_{i}\quad \text{at }r=R. & & \displaystyle\end{eqnarray}$$

In experiments, typical droplet sizes and velocities are respectively $R_{0}\sim 10~\unicode[STIX]{x03BC}\text{m}$ and $V\sim 10~\unicode[STIX]{x03BC}\text{m}~\text{s}^{-1}$ (Izri et al. Reference Izri, van der Linden, Michelin and Dauchot2014; Moerman et al. Reference Moerman, Moyses, van der Wee, Grie, van Blaaderen, Kegel, Groenewold and Brujic2017), so that the Reynolds number $Re=\unicode[STIX]{x1D70C}_{o}VR_{0}/\unicode[STIX]{x1D702}_{o}$ is small and inertia can essentially be neglected so that the velocity $\boldsymbol{u}_{i,o}$ and pressure $P_{i,o}$ , both inside and outside the droplet, satisfy the Stokes equations

(2.8a,b ) $$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D735}\boldsymbol{\cdot }\boldsymbol{u}_{i,o}=0,\quad \unicode[STIX]{x1D735}P_{i,o}=\unicode[STIX]{x1D702}_{i,o}\unicode[STIX]{x1D6FB}^{2}\boldsymbol{u}_{i,o}, & & \displaystyle\end{eqnarray}$$

with subscripts $i$ and $o$ referring to the inner and outer fluids, respectively. In the reference frame of the droplet’s centre considered here, the flow at infinity is opposite to the droplet’s translation,

(2.9) $$\begin{eqnarray}\displaystyle \boldsymbol{u}\rightarrow -{\mathcal{U}}_{\infty }\boldsymbol{e}_{z}\quad \text{for }r\rightarrow \infty , & & \displaystyle\end{eqnarray}$$

and the system of equations above is closed by enforcing the mechanical equilibrium of the droplet in the absence of inertia, i.e. the force-free condition

(2.10) $$\begin{eqnarray}\displaystyle \int _{r=R}\unicode[STIX]{x1D748}_{o}\boldsymbol{\cdot }\boldsymbol{n}\,\text{d}S=\mathbf{0}. & & \displaystyle\end{eqnarray}$$

2.2 Axisymmetric Stokes flow

Axisymmetric Stokes flows inside and outside the droplet can be recast in terms of a streamfunction $\unicode[STIX]{x1D713}_{i,o}(t,r,\unicode[STIX]{x1D707})$ , such that

(2.11) $$\begin{eqnarray}\displaystyle \boldsymbol{u}=-\frac{1}{r^{2}}\frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D713}}{\unicode[STIX]{x2202}\unicode[STIX]{x1D707}}\boldsymbol{e}_{r}-\frac{1}{r\sqrt{1-\unicode[STIX]{x1D707}^{2}}}\frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D713}}{\unicode[STIX]{x2202}r}\boldsymbol{e}_{\unicode[STIX]{x1D703}}. & & \displaystyle\end{eqnarray}$$

The general solution of the Stokes equations (2.8) for the inner and outer flows in axisymmetric spherical coordinates is given by a superposition of orthogonal modes, the so-called Lamb solution (Lamb Reference Lamb1945; Happel & Brenner Reference Happel and Brenner1983; Leal Reference Leal2007). The flow outside the droplet must converge to a finite unidirectional flow as $r\rightarrow \infty$ , and the flow inside the droplet must be regular at the origin, so the streamfunction and pressure can be written generally in the outer fluid as

(2.12) $$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D713}_{o}(t,r,\unicode[STIX]{x1D707})=\left(\frac{a_{o,1}(t)}{r}-b_{o,1}(t)r^{2}\right)(1-\unicode[STIX]{x1D707}^{2})+\mathop{\sum }_{n=2}^{\infty }\left(\frac{a_{o,n}(t)}{r^{n}}-\frac{b_{o,n}(t)}{r^{n-2}}\right)(1-\unicode[STIX]{x1D707}^{2})L_{n}^{\prime }(\unicode[STIX]{x1D707}), & \displaystyle \nonumber\\ \displaystyle & & \displaystyle\end{eqnarray}$$

(2.13) $$\begin{eqnarray}\displaystyle & \displaystyle P_{o}(t,r,\unicode[STIX]{x1D707})=\unicode[STIX]{x1D702}_{o}\mathop{\sum }_{n=2}^{\infty }2n(1-2n)b_{o,n}(t)\frac{L_{n}(\unicode[STIX]{x1D707})}{r^{n+1}}, & \displaystyle\end{eqnarray}$$

and within the droplet as

(2.14) $$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D713}_{i}(t,r,\unicode[STIX]{x1D707}) & = & \displaystyle \mathop{\sum }_{n=1}^{\infty }(a_{i,n}(t)r^{n+1}-b_{i,n}(t)r^{n+3})(1-\unicode[STIX]{x1D707}^{2})L_{n}^{\prime }(\unicode[STIX]{x1D707}),\end{eqnarray}$$

(2.15) $$\begin{eqnarray}\displaystyle P_{i}(t,r,\unicode[STIX]{x1D707}) & = & \displaystyle -2\unicode[STIX]{x1D702}_{i}\mathop{\sum }_{n=0}^{\infty }b_{i,n}(t)(n+1)(2n+3)r^{n}L_{n}(\unicode[STIX]{x1D707}),\end{eqnarray}$$

where $L_{n}$ denotes the $n$ th Legendre polynomial and the prime denotes the derivative. Note that the Stokeslet term is omitted in (2.12) since the droplet is force-free (Blake Reference Blake1971), which effectively enforces (2.10) automatically. The different modes in (2.12) correspond to singularities of increasing order in the hydrodynamic signature of the swimming droplet, and are associated with specific physical characteristics of the self-propulsion and associated fluid motion. For instance, the mode with $n=1$ carries information about the droplet self-propulsion velocity (since it is the only mode with non-zero velocity as $r\rightarrow \infty$ ), the mode with $n=2$ corresponds to a symmetric extensile flow akin to the flow excited by a force dipole (i.e. a stresslet), and so on. Note that the inner and outer hydrodynamic problems are fully determined by computing the intensity of the different modes $(a_{n,i},b_{n,i})$ and $(a_{n,o},b_{n,o})$ , respectively.

2.3 Non-dimensionalization

In the following, all quantities are non-dimensionalized using $R_{0}$ , ${\mathcal{A}}R_{0}/{\mathcal{D}}$ and $R_{0}/{\mathcal{V}}$ as reference scales for length, relative concentration of surfactant (i.e. $C-C_{\infty }$ ) and time, respectively. Here, ${\mathcal{V}}$ is the typical Marangoni velocity of a droplet in a surfactant gradient ${\mathcal{A}}$ (Anderson Reference Anderson1989):

(2.16) $$\begin{eqnarray}\displaystyle {\mathcal{V}}\equiv \frac{\unicode[STIX]{x1D6FE}_{C}{\mathcal{A}}R_{0}}{{\mathcal{D}}(2\unicode[STIX]{x1D702}_{o}+3\unicode[STIX]{x1D702}_{i})}. & & \displaystyle\end{eqnarray}$$

The pressure and viscous stress tensors are further rescaled as

(2.17a-c ) $$\begin{eqnarray}\displaystyle P_{i}\rightarrow P_{\infty }+{\displaystyle \frac{2\unicode[STIX]{x1D6FE}_{0}}{R_{0}}}+\frac{\unicode[STIX]{x1D702}_{o}{\mathcal{V}}}{R_{0}}\unicode[STIX]{x1D702}P_{i},\quad P_{o}\rightarrow P_{\infty }+\frac{\unicode[STIX]{x1D702}_{o}{\mathcal{V}}}{R_{0}}P_{o},\quad (\unicode[STIX]{x1D749}_{i},\unicode[STIX]{x1D749}_{o})\rightarrow \frac{\unicode[STIX]{x1D702}_{o}{\mathcal{V}}}{R_{0}}(\unicode[STIX]{x1D702}\unicode[STIX]{x1D749}_{i},\unicode[STIX]{x1D749}_{o}), & & \displaystyle \nonumber\\ \displaystyle & & \displaystyle\end{eqnarray}$$

where $P_{\infty }$ is the constant background pressure, $\unicode[STIX]{x1D6FE}_{0}$ the surface tension of the same spherical droplet in the absence of flow, and $\unicode[STIX]{x1D702}\equiv \unicode[STIX]{x1D702}_{i}/\unicode[STIX]{x1D702}_{o}$ is the viscosity ratio. We further denote by $U_{\infty }\equiv {\mathcal{U}}_{\infty }/{\mathcal{V}}$ the swimming velocity of the droplet. Besides $\unicode[STIX]{x1D702}$ , the physical problem is entirely characterized by two additional non-dimensional parameters, the Péclet and capillary numbers,

(2.18a,b ) $$\begin{eqnarray}\displaystyle Pe\equiv \frac{{\mathcal{V}}R_{0}}{{\mathcal{D}}},\quad Ca\equiv \frac{\unicode[STIX]{x1D6FE}_{C}{\mathcal{A}}R_{0}}{{\mathcal{D}}\unicode[STIX]{x1D6FE}_{0}}, & & \displaystyle\end{eqnarray}$$

that characterize the relative magnitude of surfactant advection and diffusion, and the relative magnitude of Marangoni and Laplace stresses, respectively.

2.4 Isotropic motionless base state

Equations (2.1)–(2.10) feature a motionless isotropic steady state given by

(2.19a-d ) $$\begin{eqnarray}\displaystyle \bar{\unicode[STIX]{x1D709}}=0,\quad \bar{\boldsymbol{u}}_{i}=\bar{\boldsymbol{u}}_{o}=\mathbf{0},\quad \bar{P_{i}}=\bar{P_{o}}=0,\quad \bar{C}=-1/r. & & \displaystyle\end{eqnarray}$$

Note that this solution (2.19) exists for any values of $Pe$ , $Ca$ and $\unicode[STIX]{x1D702}$ ; it features an isotropic surfactant distribution, no Marangoni stresses and, therefore, no droplet motion. In the following, we are interested in the existence and stability of additional non-isotropic steady states emerging in the vicinity of the base state (2.19). This amounts mathematically to finding the fundamental eigenmodes of the system. The next section focuses on finding these eigenmodes and their existence condition (i.e. the corresponding value of $Pe$ for given $Ca$ and $\unicode[STIX]{x1D702}$ ), while in § 4 we investigate the steady flows sustained by nonlinear saturation of the eigenmodes. Finally, § 5 analyses the stability of the trivial state (2.19) – the stability of the non-isotropic steady state is presented in appendix B.

3.1 Linearized equations

The linearized advection–diffusion equation reads

(3.1) $$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D6FB}^{2}C^{(1)}=-\frac{Pe}{r^{4}}\frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D713}_{o}^{(1)}}{\unicode[STIX]{x2202}\unicode[STIX]{x1D707}}, & & \displaystyle\end{eqnarray}$$

and linearized boundary conditions at the droplet interface can be written using domain perturbation, i.e. $f|_{r=R}\approx f|_{r=1}+\unicode[STIX]{x1D709}f^{\prime }|_{r=1}$ , as

(3.2a-c ) $$\begin{eqnarray}\displaystyle \frac{\unicode[STIX]{x2202}C^{(1)}}{\unicode[STIX]{x2202}r}-2\unicode[STIX]{x1D709}^{(1)}=0,\quad \unicode[STIX]{x1D713}_{i}^{(1)}=\unicode[STIX]{x1D713}_{o}^{(1)}=0,\quad \frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D713}_{i}^{(1)}}{\unicode[STIX]{x2202}r}=\frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D713}_{o}^{(1)}}{\unicode[STIX]{x2202}r}, & & \displaystyle\end{eqnarray}$$

(3.3) $$\begin{eqnarray}\displaystyle \left(\frac{\unicode[STIX]{x2202}^{2}}{\unicode[STIX]{x2202}r^{2}}-2\frac{\unicode[STIX]{x2202}}{\unicode[STIX]{x2202}r}-(1-\unicode[STIX]{x1D707}^{2})\frac{\unicode[STIX]{x2202}^{2}}{\unicode[STIX]{x2202}\unicode[STIX]{x1D707}^{2}}\right)(\unicode[STIX]{x1D713}_{o}^{(1)}-\unicode[STIX]{x1D702}\unicode[STIX]{x1D713}_{i}^{(1)})=(2+3\unicode[STIX]{x1D702})(1-\unicode[STIX]{x1D707}^{2})\frac{\unicode[STIX]{x2202}}{\unicode[STIX]{x2202}\unicode[STIX]{x1D707}}(C^{(1)}+\unicode[STIX]{x1D709}^{(1)}), & & \displaystyle \nonumber\\ \displaystyle & & \displaystyle\end{eqnarray}$$

(3.4) $$\begin{eqnarray}\displaystyle & & \displaystyle \frac{Ca}{2+3\unicode[STIX]{x1D702}}\left[-\unicode[STIX]{x1D702}P_{i}^{(1)}+P_{o}^{(1)}+2\left(\frac{\unicode[STIX]{x2202}}{\unicode[STIX]{x2202}r\unicode[STIX]{x2202}\unicode[STIX]{x1D707}}-2\frac{\unicode[STIX]{x2202}}{\unicode[STIX]{x2202}\unicode[STIX]{x1D707}}\right)(\unicode[STIX]{x1D713}_{o}^{(1)}-\unicode[STIX]{x1D702}\unicode[STIX]{x1D713}_{i}^{(1)})\right.\nonumber\\ \displaystyle & & \displaystyle \quad -\left.2(2+3\unicode[STIX]{x1D702})(C^{(1)}-2\unicode[STIX]{x1D709}^{(1)})\right]=2\unicode[STIX]{x1D709}^{(1)}-2\unicode[STIX]{x1D707}\frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D709}^{(1)}}{\unicode[STIX]{x2202}\unicode[STIX]{x1D707}}+(1-\unicode[STIX]{x1D707}^{2})\frac{\unicode[STIX]{x2202}^{2}\unicode[STIX]{x1D709}^{(1)}}{\unicode[STIX]{x2202}\unicode[STIX]{x1D707}^{2}},\end{eqnarray}$$

where superscript (1) denotes small perturbations of the base state. Because of the linearity of the Stokes equations, $\unicode[STIX]{x1D713}_{i,o}^{(1)}$ and $P_{i,o}^{(1)}$ assume the same form as (2.12)–(2.15).

The streamfunction and pressure fields can be decomposed into orthogonal modes (2.12)–(2.15). The form of the linearized equations (3.1)–(3.4) suggests that the linearized concentration $C^{(1)}(r,\unicode[STIX]{x1D707})$ and displacement $\unicode[STIX]{x1D709}^{(1)}(\unicode[STIX]{x1D707})$ also decompose in orthogonal modes of the form

(3.5a,b ) $$\begin{eqnarray}\displaystyle C^{(1)}(r,\unicode[STIX]{x1D707})=\mathop{\sum }_{n=0}^{\infty }C_{n}^{(1)}(r)L_{n}(\unicode[STIX]{x1D707}),\quad \unicode[STIX]{x1D709}^{(1)}(\unicode[STIX]{x1D707})=\unicode[STIX]{x1D709}_{0}^{(1)}+\mathop{\sum }_{n=2}^{\infty }\unicode[STIX]{x1D709}_{n}^{(1)}L_{n}(\unicode[STIX]{x1D707}), & & \displaystyle\end{eqnarray}$$

where the radial part of the basis functions of the concentration field $C_{n}^{(1)}(r)$ and constant amplitudes $\unicode[STIX]{x1D709}_{n}^{(1)}$ , together with the coefficients ( $a_{n,i}^{(1)},b_{n,o}^{(1)},b_{n,i}^{(1)},a_{n,o}^{(1)}$ ) are to be determined below and characterize each orthogonal eigenmode. In the expansion of $\unicode[STIX]{x1D709}^{(1)}(\unicode[STIX]{x1D707})$ presented in (3.5), we deliberately ignore the term $\propto L_{1}(\unicode[STIX]{x1D707})$ , since in the limit of small deformations this term corresponds to translation of the droplet, rather than deformation.

3.2 Asymptotic structure of the concentration field

Substitution of (2.12) and (3.5) into (3.1) yields an equation for $C_{n}^{(1)}(r)$ admitting the following solution:

(3.6a,b ) $$\begin{eqnarray}\displaystyle C_{0}^{(1)}(r)=\frac{c_{0}^{(1)}}{r}+d_{0}^{(1)},\quad C_{1}^{(1)}(r)=\frac{c_{1}^{(1)}}{r^{2}}+d_{1}^{(1)}r+Pe\frac{a_{o,1}^{(1)}+2b_{o,1}^{(1)}r^{3}}{2r^{3}}, & & \displaystyle\end{eqnarray}$$

(3.7) $$\begin{eqnarray}\displaystyle C_{n}^{(1)}(r)\big|_{n>1}=\frac{c_{n}^{(1)}}{r^{n+1}}+d_{n}^{(1)}r^{n}+Pe\frac{na_{o,n}^{(1)}+(n+1)b_{o,n}^{(1)}r^{2}}{2r^{n+2}}, & & \displaystyle\end{eqnarray}$$

where $c_{n}^{(1)}$ and $d_{n}^{(1)}$ are unknown constants. In the case of a self-propelling drop (i.e.  $b_{o,1}^{(1)}\neq 0$ ), this solution explicitly violates the far-field boundary condition (2.9). This is a well-known feature of advection–diffusion problems in the presence of a weak advective far-field flow: Acrivos & Taylor (Reference Acrivos and Taylor1962) used matching asymptotic expansions to demonstrate that perturbations of the solute concentration field due to particle motion are dissipated in a boundary layer located at $r\sim 1/\unicode[STIX]{x1D716}\gg 1$ , where $0<\unicode[STIX]{x1D716}\ll 1$ quantifies the velocity of the particle with respect to the surrounding fluid. This framework applies here, since asymptotically small perturbations of a motionless base state are considered, and we thus aim to construct a composite solution for the concentration field consisting of two parts: (i) a near-field part, corresponding to the immediate surroundings of the drop, $r\sim 1$ ,

(3.8) $$\begin{eqnarray}\displaystyle C(r,\unicode[STIX]{x1D707})=-1/r+\unicode[STIX]{x1D716}C^{(1)}+\unicode[STIX]{x1D716}^{2}C^{(2)}+\cdots \,, & & \displaystyle\end{eqnarray}$$

and (ii) a far-field part valid away from the drop, $\unicode[STIX]{x1D70C}\equiv r/\unicode[STIX]{x1D716}\sim 1$ ( $r\gg 1$ ),

(3.9) $$\begin{eqnarray}\displaystyle H(\unicode[STIX]{x1D70C},\unicode[STIX]{x1D707})=\unicode[STIX]{x1D716}H^{(1)}+\unicode[STIX]{x1D716}^{2}H^{(2)}+\cdots \,, & & \displaystyle\end{eqnarray}$$

as shown in figure 2. Here $H(\unicode[STIX]{x1D70C},\unicode[STIX]{x1D707})$ satisfies the rescaled advection–diffusion equation, namely,

(3.10) $$\begin{eqnarray}\displaystyle -\unicode[STIX]{x1D716}Pe\left(\frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D713}_{o}}{\unicode[STIX]{x2202}\unicode[STIX]{x1D707}}\frac{\unicode[STIX]{x2202}H}{\unicode[STIX]{x2202}\unicode[STIX]{x1D70C}}-\frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D713}_{o}}{\unicode[STIX]{x2202}\unicode[STIX]{x1D70C}}\frac{\unicode[STIX]{x2202}H}{\unicode[STIX]{x2202}\unicode[STIX]{x1D707}}\right)=\frac{\unicode[STIX]{x2202}}{\unicode[STIX]{x2202}\unicode[STIX]{x1D70C}}\left(\unicode[STIX]{x1D70C}^{2}\frac{\unicode[STIX]{x2202}H}{\unicode[STIX]{x2202}\unicode[STIX]{x1D70C}}\right)+\frac{\unicode[STIX]{x2202}}{\unicode[STIX]{x2202}\unicode[STIX]{x1D707}}\left((1-\unicode[STIX]{x1D707}^{2})\frac{\unicode[STIX]{x2202}H}{\unicode[STIX]{x2202}\unicode[STIX]{x1D707}}\right). & & \displaystyle\end{eqnarray}$$

Here $\unicode[STIX]{x1D716}$ is a small parameter that quantifies the distance to the isotropic steady state. Similarly to the work of Acrivos & Taylor (Reference Acrivos and Taylor1962), $C(r,\unicode[STIX]{x1D707})$ encapsulates the dynamics of the droplet interface, whereas $H(\unicode[STIX]{x1D70C},\unicode[STIX]{x1D707})$ is determined by advection of the concentration disturbances imparted by the droplet. Naturally, $C(r,\unicode[STIX]{x1D707})$ and $H(\unicode[STIX]{x1D70C},\unicode[STIX]{x1D707})$ must yield identical results in the matching region occurring at $1\ll r\ll 1/\unicode[STIX]{x1D716}$ (or, equivalently, $\unicode[STIX]{x1D716}\ll \unicode[STIX]{x1D70C}\ll 1$ ) (Holmes Reference Holmes1995). Note that although (3.8) and (3.9) include higher-order terms in $\unicode[STIX]{x1D716}$ , only $C^{(1)}$ and $H^{(1)}$ are relevant in the context of linear analysis.

Figure 2.

Asymptotic structure of the surfactant concentration field around the droplet. The near-field solution $C(r,\unicode[STIX]{x1D707})$ is matched with the far-field solution $H(r,\unicode[STIX]{x1D707})$ in the intermediate region where $1\ll r\ll \unicode[STIX]{x1D716}^{-1}$ .

3.3 Concentration field far from the translating drop and asymptotic matching

Linearization of the rescaled advection–diffusion equation (3.10) about the base state, namely, about $\bar{H}=0$ , yields

(3.11) $$\begin{eqnarray}\displaystyle F(H^{(1)}) & \equiv & \displaystyle 2Pe\,b_{o,1}^{(1)}\left(\unicode[STIX]{x1D707}\frac{\unicode[STIX]{x2202}H^{(1)}}{\unicode[STIX]{x2202}\unicode[STIX]{x1D70C}}+\frac{1-\unicode[STIX]{x1D707}^{2}}{\unicode[STIX]{x1D70C}}\frac{\unicode[STIX]{x2202}H^{(1)}}{\unicode[STIX]{x2202}\unicode[STIX]{x1D707}}\right)\nonumber\\ \displaystyle & & \displaystyle +\,\frac{1}{\unicode[STIX]{x1D70C}^{2}}\left[\frac{\unicode[STIX]{x2202}}{\unicode[STIX]{x2202}\unicode[STIX]{x1D70C}}\left(\unicode[STIX]{x1D70C}^{2}\frac{\unicode[STIX]{x2202}H}{\unicode[STIX]{x2202}\unicode[STIX]{x1D70C}}\right)+\frac{\unicode[STIX]{x2202}}{\unicode[STIX]{x2202}\unicode[STIX]{x1D707}}\left((1-\unicode[STIX]{x1D707}^{2})\frac{\unicode[STIX]{x2202}H}{\unicode[STIX]{x2202}\unicode[STIX]{x1D707}}\right)\right]=0.\end{eqnarray}$$

Its solution that decays as $\unicode[STIX]{x1D70C}\rightarrow \infty$ can be written as (Acrivos & Taylor Reference Acrivos and Taylor1962)

(3.12) $$\begin{eqnarray}\displaystyle H^{(1)}(\unicode[STIX]{x1D70C},\unicode[STIX]{x1D707})=\frac{\text{e}^{-\unicode[STIX]{x1D70C}_{s}\unicode[STIX]{x1D707}}}{\sqrt{|\unicode[STIX]{x1D70C}_{s}|}}\mathop{\sum }_{n=0}^{\infty }h_{n}^{(1)}\text{K}_{n+1/2}(|\unicode[STIX]{x1D70C}_{s}|)L_{n}(\unicode[STIX]{x1D707}), & & \displaystyle\end{eqnarray}$$

where $h_{n}^{(1)}$ are unknown constants to be determined in the matching process with the inner solution, $\unicode[STIX]{x1D70C}_{s}\equiv Pe\,b_{o,1}^{(1)}\unicode[STIX]{x1D70C}$ , and $\text{K}_{n}(x)$ denotes the modified Bessel function of the second kind of order $n$ . We note that the direction of the droplet motion along the symmetry axis is determined by the sign of the constant amplitude $b_{o,1}^{(1)}$ . Since there is no physical difference between the two directions of motion, we assume $b_{o,1}^{(1)}\geqslant 0$ in what follows.

We use Van Dyke’s matching rule (Holmes Reference Holmes1995) to match solutions (3.6), (3.7) and (3.12) in the region $\unicode[STIX]{x1D716}\ll \unicode[STIX]{x1D70C}\ll 1$ . More specifically, $C(r,\unicode[STIX]{x1D707})$ is expressed in terms of $\unicode[STIX]{x1D70C}$ and both $C$ and $H$ are expanded in powers of $\unicode[STIX]{x1D70C}\ll 1$ . Note that only the terms linear in $\unicode[STIX]{x1D716}$ , $\unicode[STIX]{x1D70C}$ or $\unicode[STIX]{x1D716}/\unicode[STIX]{x1D70C}$ can be matched at the leading order of expansion. As a consequence,

(3.13a-c ) $$\begin{eqnarray}\displaystyle d_{0}^{(1)}=Pe\,b_{o,1}^{(1)},\quad h_{0}^{(1)}=-\sqrt{\frac{2}{\unicode[STIX]{x03C0}}}\,Pe\,b_{o,1}^{(1)},\quad d_{n}^{(1)}=h_{n}^{(1)}=0\quad \text{for }n>0. & & \displaystyle\end{eqnarray}$$

3.4 Solvability condition

Substitution of $\unicode[STIX]{x1D713}_{i}^{(1)}$ , $\unicode[STIX]{x1D713}_{o}^{(1)}$ and $C^{(1)}$ given by (2.12), (2.14), (3.6) and (3.7) into the boundary conditions (3.2)–(3.4) and subsequent projection of the result onto the $n$ th Legendre polynomial yields a sequence of sets of homogeneous linear algebraic equations for the amplitudes $a_{i,n}^{(1)}$ , $b_{i,n}^{(1)}$ , $a_{o,n}^{(1)}$ , $b_{o,n}^{(1)}$ , $c_{n}^{(1)}$ and $\unicode[STIX]{x1D709}_{n}^{(1)}$ . For given $n$ , the solvability condition, i.e. the existence condition for a non-trivial solution to the linearized problem, reads

(3.14) $$\begin{eqnarray}\displaystyle Pe=Pe_{n}\equiv \left\{\begin{array}{@{}ll@{}}4,\quad & n=1,\\ {\displaystyle \frac{4}{2+3\unicode[STIX]{x1D702}}}\left[(n+1)(2n+1)(1+\unicode[STIX]{x1D702})-Ca\left(1+\unicode[STIX]{x1D702}{\displaystyle \frac{n-1}{n+2}}\right)\right],\quad & n>1.\end{array}\right. & & \displaystyle\end{eqnarray}$$

Equations (3.14) determine when the linearized problem yields a non-trivial time-independent (i.e. neutrally stable) solution. In § 5, we will show that $Pe_{n}$ is the minimum Péclet number below (respectively above) which the isotropic state is linearly stable (respectively unstable) to perturbations along the $n$ th mode presented above. In other words, $Pe_{n}$ represents the threshold of the $n$ th mode of monotonic instability. This motivates referring to $Pe_{n}$ as the instability threshold for mode $n$ in the following. Depending on the value of the capillary number $Ca$ , one of the first two modes in (3.14) features the lowest instability threshold. As a consequence, the next stages of the analysis are focused exclusively on these first two instability modes, namely, the cases of $n=1$ and $n=2$ , shown in figure 3.

The first eigenmode ( $n=1$ , $Pe_{1}=4$ ) reads

(3.15a,b ) $$\begin{eqnarray}\displaystyle b_{i,0}^{(1)}=\frac{Pe_{1}A_{1}}{3\unicode[STIX]{x1D702}}(2+3\unicode[STIX]{x1D702}),\quad \frac{3}{2}a_{i,1}^{(1)}=\frac{3}{2}b_{i,1}^{(1)}=a_{o,1}^{(1)}=b_{o,1}^{(1)}=A_{1}, & & \displaystyle\end{eqnarray}$$

(3.16a-d ) $$\begin{eqnarray}\displaystyle c_{0}^{(1)}=0,\quad d_{0}^{(1)}=Pe_{1}A_{1},\quad c_{1}^{(1)}=-\frac{3Pe_{1}A_{1}}{4},\quad d_{1}^{(1)}=0, & & \displaystyle\end{eqnarray}$$

and physically corresponds to a polar concentration field at the surface of the droplet $C^{(1)}\propto \unicode[STIX]{x1D707}$ , which maintains its steady translation. The corresponding flow field is that of a translating droplet, i.e. the superposition of a steady flow with a source dipole singularity to enforce the impermeability condition. Note that the instability threshold $Pe_{1}=4$ does not include the viscosity ratio $\unicode[STIX]{x1D702}$ , since we define dimensionless velocity based on the terminal velocity of a droplet in an imposed gradient of surfactant concentration (2.16) (Anderson Reference Anderson1989).

The second eigenmode ( $n=2$ ) can be written as

(3.17a,b ) $$\begin{eqnarray}\displaystyle a_{i,2}^{(1)}=b_{i,2}^{(1)}=a_{o,2}^{(1)}=b_{o,2}^{(1)}=A_{2},\quad \unicode[STIX]{x1D709}_{2}^{(1)}=Ca\,A_{2}\frac{4+\unicode[STIX]{x1D702}}{2(2+3\unicode[STIX]{x1D702})}, & & \displaystyle\end{eqnarray}$$

(3.18a,b ) $$\begin{eqnarray}\displaystyle c_{2}^{(1)}=-\frac{A_{2}}{3}\left(7Pe_{2}+Ca\frac{4+\unicode[STIX]{x1D702}}{2+3\unicode[STIX]{x1D702}}\right),\quad d_{2}^{(1)}=0. & & \displaystyle\end{eqnarray}$$

In that case, the concentration field $C^{(1)}\propto L_{2}(\unicode[STIX]{x1D707})$ is front–back symmetric and cannot drive any net droplet motion. Instead, an extensile flow is forced by the Marangoni stress. Outside the droplet, it takes the same form as the second mode of the classical squirmer model (Blake Reference Blake1971) and consists of a stresslet singularity (i.e. symmetric force dipole) and source quadrupole.

The amplitudes of the eigenmodes of the system, $A_{1}\geqslant 0$ in (3.15), (3.16) or $A_{2}$ in (3.17), (3.18), remain naturally undetermined within this linear framework. The weakly nonlinear analysis, which is the focus of the next section, will provide these saturation amplitudes in the vicinity of the critical conditions $Pe=Pe_{n}$ .

Figure 3.

Two instability modes competing for the minimal value of the instability threshold $Pe_{n}$ . (a) Evolution of the critical Péclet number with the capillary number for $\unicode[STIX]{x1D702}=1$ . (b) Flow field corresponding to the first and the second mode of instability, respectively.

Finally, we note that (3.14) echoes the result of the stability analysis in the case of chemically active isotropic particles developed by Michelin et al. (Reference Michelin, Lauga and Bartolo2013). In particular, Michelin et al. have demonstrated that the onset of spontaneous self-propulsion of active spherical particles also corresponds to $Pe=4$ . We argue that the persisting value of the instability threshold is related to the choice of dimensionless velocity. Both this paper and the work of Michelin et al. (Reference Michelin, Lauga and Bartolo2013) define dimensionless velocity based on the velocity of an active drop/particle in an external concentration gradient ${\mathcal{A}}/{\mathcal{D}}$ . In both cases, self-propulsion and the resulting advection of the concentration field generates a front–back concentration contrast. The onset of the drop or particle motion then corresponds in both cases to a fixed ratio of advective (i.e. destabilizing) and diffusive (i.e. stabilizing) terms, resulting in a fixed value of the Péclet number, $Pe=4$ .

4.1 Steady self-propulsion ( $n=1$ )

In the case of $n=1$ , $Pe_{1}\equiv 4$ and the leading-order flow is given by the first squirming mode with coefficients presented in (3.15) and (3.16). This squirming mode corresponds to a self-propelling droplet and, to determine the self-propulsion velocity, we now consider the problem at $\unicode[STIX]{x1D716}^{2}$ featuring quadratic interactions of the flow and concentration fields obtained at  $\unicode[STIX]{x1D716}$ .

4.1.1 Concentration field around the droplet

The quadratic approximation of the advection–diffusion equation (i.e. retaining only $O(\unicode[STIX]{x1D716}^{2})$ terms) reads

(4.3) $$\begin{eqnarray}\displaystyle -Pe_{n}\frac{1}{r^{4}}\frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D713}_{o}^{(2)}}{\unicode[STIX]{x2202}\unicode[STIX]{x1D707}}-\unicode[STIX]{x1D6FB}^{2}C^{(2)}=\frac{Pe_{n}}{r^{2}}\left(\frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D713}_{o}^{(1)}}{\unicode[STIX]{x2202}\unicode[STIX]{x1D707}}\frac{\unicode[STIX]{x2202}C^{(1)}}{\unicode[STIX]{x2202}r}-\frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D713}_{o}^{(1)}}{\unicode[STIX]{x2202}r}\frac{\unicode[STIX]{x2202}C^{(1)}}{\unicode[STIX]{x2202}\unicode[STIX]{x1D707}}\right)+\frac{\unicode[STIX]{x1D6FF}}{r^{4}}\frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D713}_{o}^{(1)}}{\unicode[STIX]{x2202}\unicode[STIX]{x1D707}}. & & \displaystyle\end{eqnarray}$$

Using (3.15), (3.16) and the resulting forms of $C^{(1)}$ and $\unicode[STIX]{x1D713}_{i,o}^{(1)}$ , the inhomogeneous right-hand side of (4.3) includes non-zero projections onto the zeroth, first and second Legendre harmonics. Accordingly, the angular component of $C^{(2)}(r,\unicode[STIX]{x1D707})$ is given by the first three Legendre polynomials, namely, $C^{(2)}(r,\unicode[STIX]{x1D707})=\sum _{n=0}^{2}C_{n}^{(2)}(r)L_{n}(\unicode[STIX]{x1D707})$ , with

(4.4) $$\begin{eqnarray}\displaystyle & \displaystyle C_{0}^{(2)}(r)=\frac{c_{0}^{(2)}}{r}+d_{0}^{(2)}-Pe_{1}^{2}A_{1}^{2}\frac{8-15r+20r^{3}+80r^{6}}{120r^{5}}, & \displaystyle\end{eqnarray}$$

(4.5) $$\begin{eqnarray}\displaystyle & \displaystyle C_{1}^{(2)}(r)=\frac{c_{1}^{(2)}}{r^{2}}+d_{1}^{(2)}r+\unicode[STIX]{x1D6FF}A_{1}\frac{1+2r^{3}}{2r^{3}}+Pe_{1}\frac{a_{o,1}^{(2)}+2b_{o,1}^{(2)}r^{3}}{2r^{3}}, & \displaystyle\end{eqnarray}$$

(4.6) $$\begin{eqnarray}\displaystyle & \displaystyle C_{2}^{(2)}(r)=\frac{c_{2}^{(2)}}{r^{3}}+d_{2}^{(2)}r^{2}+Pe_{1}\frac{2a_{o,2}^{(2)}+3b_{o,2}^{(2)}r^{2}}{2r^{4}}-Pe_{1}^{2}A_{1}^{2}\frac{10-21r+70r^{3}-42r^{4}+28r^{6}}{84r^{5}}, & \displaystyle \nonumber\\ \displaystyle & & \displaystyle\end{eqnarray}$$

where $c_{n}^{(2)}$ and $d_{n}^{(2)}$ are unknown constant amplitudes to be determined in the matching process with the far-field boundary layer.

Equation (3.10) reads at $O(\unicode[STIX]{x1D716}^{2})$ as

(4.7) $$\begin{eqnarray}\displaystyle F(H^{(2)})=-2(Pe_{1}b_{o,1}^{(2)}+\unicode[STIX]{x1D6FF}A_{1})\left(\unicode[STIX]{x1D707}\frac{\unicode[STIX]{x2202}H^{(1)}}{\unicode[STIX]{x2202}\unicode[STIX]{x1D70C}}+\frac{1-\unicode[STIX]{x1D707}^{2}}{\unicode[STIX]{x1D70C}}\frac{\unicode[STIX]{x2202}H^{(1)}}{\unicode[STIX]{x2202}\unicode[STIX]{x1D707}}\right), & & \displaystyle\end{eqnarray}$$

where the linear operator $F$ is defined in (3.11). Since $A_{1}\geqslant 0$ , the solution of (4.7) which decays as $\unicode[STIX]{x1D70C}\rightarrow \infty$ reads

(4.8) $$\begin{eqnarray}\displaystyle H^{(2)}(\unicode[STIX]{x1D70C},\unicode[STIX]{x1D707}) & = & \displaystyle \frac{Pe_{1}b_{o,1}^{(2)}+\unicode[STIX]{x1D6FF}A_{1}}{2}\left(2(1+\unicode[STIX]{x1D707})+3\unicode[STIX]{x1D707}\frac{1+\unicode[STIX]{x1D70C}_{s}}{\unicode[STIX]{x1D70C}_{s}^{2}}\right)\text{e}^{-\unicode[STIX]{x1D70C}_{s}(1+\unicode[STIX]{x1D707})}\nonumber\\ \displaystyle & & \displaystyle +\,\frac{\text{e}^{-\unicode[STIX]{x1D70C}_{s}\unicode[STIX]{x1D707}}}{\sqrt{\unicode[STIX]{x1D70C}_{s}}}\mathop{\sum }_{n=0}^{\infty }h_{n}^{(2)}\text{K}_{n+1/2}(\unicode[STIX]{x1D70C}_{s})L_{n}(\unicode[STIX]{x1D707}),\end{eqnarray}$$

where $h_{n}^{(2)}$ are unknown constant amplitudes to be determined in the matching process, and $\unicode[STIX]{x1D70C}_{s}\equiv Pe_{1}A_{1}\unicode[STIX]{x1D70C}\geqslant 0$ .

Asymptotic matching of $C(r,\unicode[STIX]{x1D707})$ and $H(r,\unicode[STIX]{x1D707})$ at $O(\unicode[STIX]{x1D716}^{2})$ in the region $\unicode[STIX]{x1D716}\ll \unicode[STIX]{x1D70C}\ll 1$ is achieved by expressing $C(r,\unicode[STIX]{x1D707})=-1/r+\unicode[STIX]{x1D716}C^{(1)}+\unicode[STIX]{x1D716}^{2}C^{(2)}$ in terms of $\unicode[STIX]{x1D70C}$ , and expanding both $C(\unicode[STIX]{x1D70C},\unicode[STIX]{x1D707})$ and $H(\unicode[STIX]{x1D70C},\unicode[STIX]{x1D707})$ in powers of $\unicode[STIX]{x1D70C}$ . Linear and quadratic terms in $\unicode[STIX]{x1D716}$ , $\unicode[STIX]{x1D70C}$ or $\unicode[STIX]{x1D716}/\unicode[STIX]{x1D70C}$ must now be matched, leading to

(4.9a-c ) $$\begin{eqnarray}\displaystyle d_{0}^{(2)}=Pe_{1}b_{o,1}^{(2)}+\unicode[STIX]{x1D6FF}A_{1},\quad d_{1}^{(2)}=-Pe_{1}^{2}A_{1}^{2},\quad d_{n}^{(2)}=0\quad \text{for }n>1, & & \displaystyle\end{eqnarray}$$

(4.10a,b ) $$\begin{eqnarray}\displaystyle h_{1}^{(2)}=-3\frac{Pe_{1}b_{o,1}^{(2)}+\unicode[STIX]{x1D6FF}A_{1}}{\sqrt{2\unicode[STIX]{x03C0}}},\quad h_{n}^{(2)}=0\quad \text{for }n\neq 1. & & \displaystyle\end{eqnarray}$$

4.1.2 Solvability condition

Combination of (4.9) and (4.4)–(4.6) yields $C^{(2)}(r,\unicode[STIX]{x1D707})$ . Expanding the boundary conditions (2.2), (2.6) and (2.7) at $O(\unicode[STIX]{x1D716}^{2})$ and projecting the result onto the first three Legendre polynomials provides a set of inhomogeneous linear algebraic equations for the amplitudes $a_{i,n}^{(2)}$ , $b_{i,n}^{(2)}$ , $a_{o,n}^{(2)}$ , $b_{o,n}^{(2)}$ , $c_{n}^{(2)}$ and $\unicode[STIX]{x1D709}_{n}^{(2)}$ . The solvability condition of this set of equations reads

(4.11) $$\begin{eqnarray}\displaystyle A_{1}(\unicode[STIX]{x1D6FF}-32A_{1})=0. & & \displaystyle\end{eqnarray}$$

Equation (4.11) implies that two branches of steady solutions of the nonlinear problem (2.3)–(2.10) exist near $Pe_{1}$ : the first branch, given by $A_{1}=0$ , corresponds to a motionless droplet and is in fact simply the isotropic steady state (2.19) already discussed; the second branch describes a self-propelling drop with finite velocity ( $A_{1}=\unicode[STIX]{x1D6FF}/32>0$ ). Note that in the case of self-propulsion, the velocity of the droplet grows linearly with $Pe-Pe_{1}$ (i.e. $A_{1}\propto \unicode[STIX]{x1D6FF}$ ), suggesting that the onset of the droplet motion is a transcritical bifurcation. Also recall that the self-propelling mode is associated with no droplet deformation. As a result, the solvability condition (4.11) does not include $Ca$ , i.e. in the leading order, deformability does not affect droplet self-propulsion velocity.

Following the steps of the analysis above, it is easy to demonstrate that, in the case of $A_{1}\leqslant 0$ , the counterpart of the solvability condition (4.11) reads

(4.12) $$\begin{eqnarray}\displaystyle A_{1}(\unicode[STIX]{x1D6FF}+32A_{1})=0, & & \displaystyle\end{eqnarray}$$

implying that, similarly to the case of $A_{1}\geqslant 0$ , the steady problem is solvable either for $A_{1}$ = 0 or for $\unicode[STIX]{x1D6FF}>0$ , i.e. above the threshold of the Marangoni instability.

4.1.3 Effect of deformability on the droplet’s self-propulsion

Equation (4.11) provides the leading-order evolution of the droplet velocity for the self-propelled steady state near the onset of propulsion, $Pe=Pe_{1}$ , and was obtained by considering the solvability condition of the problem including corrections up to $O(\unicode[STIX]{x1D716}^{2})$ . We now use the same approach to extend the expansion of the different equations up to $O(\unicode[STIX]{x1D716}^{3})$ and obtain the quadratic correction for the self-propulsion velocity.

When the solvability condition (4.11) is satisfied, the non-trivial solution of the problem at $O(\unicode[STIX]{x1D716}^{2})$ can be written as

(4.13a,b ) $$\begin{eqnarray}\displaystyle b_{i,0}^{(2)}={\displaystyle \frac{(2+3\unicode[STIX]{x1D702})(29\unicode[STIX]{x1D6FF}^{2}+10240B_{1})}{7680\unicode[STIX]{x1D702}}},\quad \frac{3}{2}a_{i,1}^{(2)}=\frac{3}{2}b_{i,1}^{(2)}=a_{o,1}^{(2)}=b_{o,1}^{(2)}=B_{1}, & & \displaystyle\end{eqnarray}$$

(4.14a,b ) $$\begin{eqnarray}\displaystyle a_{i,2}^{(2)}=b_{i,2}^{(2)}=a_{o,2}^{(2)}=b_{o,2}^{(2)}=-\frac{33\unicode[STIX]{x1D6FF}^{2}}{896(Pe_{2}-4)},\quad \unicode[STIX]{x1D709}_{2}^{(2)}=-\frac{33Ca\unicode[STIX]{x1D6FF}^{2}(4+\unicode[STIX]{x1D702})}{1792(2+3\unicode[STIX]{x1D702})(Pe_{2}-4)}, & & \displaystyle \nonumber\\ \displaystyle & & \displaystyle\end{eqnarray}$$

(4.15a-d ) $$\begin{eqnarray}\displaystyle c_{0}^{(2)}=-\frac{\unicode[STIX]{x1D6FF}^{2}}{128},\quad d_{0}^{(2)}=4B_{1}+\frac{\unicode[STIX]{x1D6FF}^{2}}{32},\quad c_{1}^{(2)}=-\left(\frac{\unicode[STIX]{x1D6FF}^{2}}{32}+3B_{1}\right),\quad d_{1}^{(2)}=-\frac{\unicode[STIX]{x1D6FF}^{2}}{64}, & & \displaystyle \nonumber\\ \displaystyle & & \displaystyle\end{eqnarray}$$

(4.16a,b ) $$\begin{eqnarray}\displaystyle c_{2}^{(2)}=3\unicode[STIX]{x1D6FF}^{2}\frac{20(12+17\unicode[STIX]{x1D702})+3Ca(4+\unicode[STIX]{x1D702})}{896(2+3\unicode[STIX]{x1D702})(Pe_{2}-4)},\quad d_{2}^{(2)}=0, & & \displaystyle\end{eqnarray}$$

where the remaining unknown constant $B_{1}$ is determined from the solvability condition at $O(\unicode[STIX]{x1D716}^{3})$ . It should be noted in (4.13)–(4.16) that the leading-order deformation of the self-propelling droplet is $O(\unicode[STIX]{x1D716}^{2})$ (recall that $\unicode[STIX]{x1D709}^{(1)}=0$ ) and always corresponds to an oblate shape ( $\unicode[STIX]{x1D709}_{2}^{(2)}<0$ ). The leading-order volume change is $O(\unicode[STIX]{x1D716}^{4})$ , so that volume conservation is automatically enforced up to quartic order in $\unicode[STIX]{x1D716}$ .

To obtain a correction to the droplet self-propulsion velocity, the weakly nonlinear analysis must be carried out up to the cubic order. The solution procedure of the problem at $\unicode[STIX]{x1D716}^{3}$ remains exactly the same as for the lower-order problems, and several intermediate results of the derivation are provided in appendix A. As for the $\unicode[STIX]{x1D716}^{2}$ problem, the solvability condition at $\unicode[STIX]{x1D716}^{3}$ provides information about the $O(\unicode[STIX]{x1D716}^{2})$ droplet velocity, namely,

(4.17) $$\begin{eqnarray}\displaystyle U_{\infty } & = & \displaystyle 2(\unicode[STIX]{x1D716}A_{1}+\unicode[STIX]{x1D716}^{2}B_{1})\nonumber\\ \displaystyle & = & \displaystyle \frac{Pe-4}{16}\nonumber\\ \displaystyle & & \displaystyle -\,\frac{(Pe-4)^{2}[4(2+3\unicode[STIX]{x1D702})(343\,073+325\,872\unicode[STIX]{x1D702})-35Ca(4+\unicode[STIX]{x1D702})(932+2883\unicode[STIX]{x1D702})]}{1254\,400(2+3\unicode[STIX]{x1D702})^{2}(Pe_{2}-4)},\nonumber\\ \displaystyle & & \displaystyle\end{eqnarray}$$

with $Pe_{2}(Ca,\unicode[STIX]{x1D702})=[60(1+\unicode[STIX]{x1D702})-Ca(4+\unicode[STIX]{x1D702})]/(2+3\unicode[STIX]{x1D702})$ (see (3.14)). One can immediately observe that at this order $\unicode[STIX]{x2202}_{Ca}U_{\infty }>0$ (respectively $\unicode[STIX]{x2202}_{Ca}U_{\infty }<0$ ) when $\unicode[STIX]{x1D702}>1.0134$ (respectively $\unicode[STIX]{x1D702}<1.0134$ ). The deformability of the droplet’s interface therefore affects differently the self-propulsion of droplets that are more or less viscous than the surrounding fluid: roughly speaking, (4.17) states that deformability enhances self-propulsion of viscous droplets, but hinders propulsion of drops that are less viscous than the surrounding fluid (figure 4).

Figure 4.

Effect of deformability ( $Ca$ ) and relative viscosity ( $\unicode[STIX]{x1D702}$ ) on the self-propulsion velocity of an active droplet near the onset of the Marangoni instability as given in (4.17). The solid line represents the case of a non-deformable droplet ( $Ca=0$ ) with $\unicode[STIX]{x1D702}=0.5$ or $\unicode[STIX]{x1D702}=2$ (the two cases are essentially indistinguishable due to weak dependence of $U_{\infty }$ on $\unicode[STIX]{x1D702}$ in the limit of $Ca=0$ ). The dotted and dashed lines represent deformable droplets ( $Ca=10$ ), with $\unicode[STIX]{x1D702}=0.5$ and $\unicode[STIX]{x1D702}=2$ , respectively.

This result can be interpreted as follows. Steady self-propulsion of the droplet occurs when the Stokes drag balances the thrust generated by Marangoni stresses, which result from concentration gradients at the surface. In our analysis, this balance is represented by a saturated self-propelling eigenmode, where the saturation comes from nonlinear terms implementing weakly nonlinear interaction of the different components of the solution. At $\unicode[STIX]{x1D716}^{3}$ , the cumulative Marangoni stress on the drop is $\propto C_{1}^{(3)}(r=1)$ , see (A 2), which includes the term $\propto 1/(Pe_{2}-4)$ resulting from the interaction of the front–back symmetric component ( $n=2$ ) of the $O(\unicode[STIX]{x1D716}^{2})$ concentration mode with the leading-order (i.e. $O(\unicode[STIX]{x1D716})$ ) flow associated with self-propulsion. The second critical Péclet number $Pe_{2}$ is a decreasing function of $Ca$ ; thus increasing deformability enhances the term $\propto 1/(Pe_{2}-4)$ ; moreover, when $r=1$ , this term is positive (respectively negative) for $\unicode[STIX]{x1D702}<39/29$ (respectively $\unicode[STIX]{x1D702}>39/29$ ). That is, an increase in capillary number tends to increase the front–back concentration gradient of surfactant (and Marangoni forcing) for viscous droplets, while deformability tends to reduce them for less viscous droplets.

In addition, (4.13)–(4.16) establish that self-propulsion at $Ca>0$ is always accompanied by droplet deformations of order $\unicode[STIX]{x1D716}^{2}$ with an oblate shape, $\unicode[STIX]{x1D709}^{(2)}<0$ . Oblate deformations are known to increase the Stokes drag on a steadily moving droplet (Matunobu Reference Matunobu1966), namely,

(4.18) $$\begin{eqnarray}\displaystyle F_{D}\propto 1-\unicode[STIX]{x1D716}^{2}\unicode[STIX]{x1D709}_{2}^{(2)}\frac{3\unicode[STIX]{x1D702}^{2}-\unicode[STIX]{x1D702}+8}{5(1+\unicode[STIX]{x1D702})(2+3\unicode[STIX]{x1D702})}. & & \displaystyle\end{eqnarray}$$

In our case, the magnitude of droplet deformations given by $|\unicode[STIX]{x1D709}^{(2)}|$ decreases with increasing $\unicode[STIX]{x1D702}$ , i.e. less viscous droplets deform more and experience higher drag and a lower velocity for a given Marangoni forcing. Both effects of $Ca$ (modified Marangoni forcing and modified viscous drag) therefore reinforce each other: drops with smaller $\unicode[STIX]{x1D702}$ deform more than their viscous counterparts and, thus, experience higher drag, and they also experience a reduced Marangoni forcing, leading to the non-monotonic effect of deformability illustrated in figure 4. Interestingly, oblate deformations of active droplets were also predicted in the limit of $Pe\gg 1$ by Golovin et al. (Reference Golovin, Gupalo and Ryazantsev1989) and in the presence of a chemical reaction in the bulk fluid by Yoshinaga (Reference Yoshinaga2014). An alternative approach based on reaction–diffusion equations also yields oblate deformations of the reacting domain (Shitara et al. Reference Shitara, Hiraiwa and Ohta2011).

4.2 Symmetric extensile flow ( $n=2$ )

The previous section focused on the analysis of the self-propelling mode of instability ( $n=1$ ). The present section focuses now on the second mode characterized by a symmetric extensile flow (see figure 3). In the case of $n=2$ , the solvability condition (3.14) yields $Pe=Pe_{2}$ , and the leading-order flow is given by the second squirming mode with an unknown amplitude $A_{2}$ introduced in (3.17) and (3.18). In order to determine $A_{2}$ , the solvability condition must be established for the $\unicode[STIX]{x1D716}^{2}$ problem featuring quadratic interactions of the flow and concentration fields obtained in § 3.

4.2.1 Concentration distribution around the droplet

For $n=2$ , $C^{(1)}(r,\unicode[STIX]{x1D707})\propto L_{2}(\unicode[STIX]{x1D707})$ , and quadratic nonlinearities in (4.3) produce non-zero projections of the $O(\unicode[STIX]{x1D716}^{2})$ concentration field onto the zeroth, second and fourth Legendre harmonics, that is,

(4.19) $$\begin{eqnarray}\displaystyle & C^{(2)}(r,\unicode[STIX]{x1D707})=\mathop{\sum }_{n=0}^{2}C_{2n}^{(2)}(r)L_{2n}(\unicode[STIX]{x1D707}), & \displaystyle\end{eqnarray}$$

with

(4.20) $$\begin{eqnarray}\displaystyle C_{0}^{(2)}(r) & = & \displaystyle \frac{c_{0}^{(2)}}{r}-Pe_{2}^{2}A_{2}^{2}\frac{180-490r+126r^{2}+735r^{3}-630r^{4}}{1050r^{7}}\nonumber\\ \displaystyle & & \displaystyle +\,Pe_{2}Ca\,A_{2}^{2}\frac{(4+\unicode[STIX]{x1D702})(2-3r^{2})}{30r^{6}(2+3\unicode[STIX]{x1D702})},\end{eqnarray}$$

(4.21) $$\begin{eqnarray}\displaystyle C_{2}^{(2)}(r) & = & \displaystyle \frac{c_{2}^{(2)}}{r^{3}}+\unicode[STIX]{x1D6FF}A_{2}\frac{2+3r^{2}}{2r^{4}}+Pe_{2}\frac{2a_{o,2}^{(2)}+3b_{o,2}^{(2)}r^{2}}{2r^{4}}+2Pe_{2}Ca\,A_{2}^{2}\frac{(4+\unicode[STIX]{x1D702})(1-3r^{2})}{21r^{6}(2+3\unicode[STIX]{x1D702})}\nonumber\\ \displaystyle & & \displaystyle -\,Pe_{2}^{2}A_{2}^{2}\frac{875-2450r+225r^{2}+7350r^{3}+756r^{4}+3780r^{4}\log r}{3675r^{7}},\end{eqnarray}$$

(4.22) $$\begin{eqnarray}\displaystyle C_{4}^{(2)}(r) & = & \displaystyle \frac{c_{4}^{(2)}}{r^{5}}+Pe_{2}\frac{4a_{o,4}^{(2)}+5b_{o,4}^{(2)}r^{2}}{2r^{6}}+3Pe_{2}Ca\,A_{2}^{2}\frac{(4+\unicode[STIX]{x1D702})(4+9r^{2})}{70r^{6}(2+3\unicode[STIX]{x1D702})}\nonumber\\ \displaystyle & & \displaystyle -\,Pe_{2}^{2}A_{2}^{2}\frac{2016-6468r+616r^{2}-14\,553r^{3}+3564r^{4}+5544r^{2}\log r}{5390r^{7}},\qquad \quad\end{eqnarray}$$

where $c_{n}^{(2)}$ are constant amplitudes to be determined. In the absence of any translation of the droplet, the rescaled advection–diffusion equation (3.10) at $\unicode[STIX]{x1D716}^{2}$ reduces to Laplace’s equation, implying that for $n=2$ advection of the perturbations in the far field is negligible (i.e. no boundary layer is needed here). Equivalently, it is now possible to satisfy the far-field and near-field boundary conditions for $C^{(2)}(r,\unicode[STIX]{x1D707})$ , which must decay as $r\rightarrow \infty$ .

4.2.2 Boundary and solvability conditions

Expanding the boundary conditions (2.2) and (2.6), (2.7) up to $O(\unicode[STIX]{x1D716}^{2})$ , and substituting for $C^{(1)}$ , $C^{(2)}$ and the corresponding flow field, projection of the result onto the first five Legendre polynomials provides a set of inhomogeneous linear algebraic equations for the amplitudes $a_{i,n}^{(2)}$ , $b_{i,n}^{(2)}$ , $a_{o,n}^{(2)}$ , $b_{o,n}^{(2)}$ , $c_{n}^{(2)}$ and $\unicode[STIX]{x1D709}_{n}^{(2)}$ . In particular, projection of the kinematic boundary condition onto $L_{0}(\unicode[STIX]{x1D707})$ yields

(4.23) $$\begin{eqnarray}\displaystyle Ca\,A_{2}^{2}\frac{6(4+\unicode[STIX]{x1D702})}{5(2+3\unicode[STIX]{x1D702})}=0. & & \displaystyle\end{eqnarray}$$

That is, a symmetric steady flow field around an active droplet is not possible for $Ca=O(1)$ , but may exist in the limit of a weakly deformable droplet, i.e. $Ca\sim O(\unicode[STIX]{x1D716})$ . Indeed, a weakly deformable droplet is spherical in the leading order, $\unicode[STIX]{x1D709}^{(1)}=0$ , and features instability thresholds independent of $Ca$ . In the particular case of $Ca\sim \unicode[STIX]{x1D716}$ , the $O(Ca)$ terms are pushed to the next order of asymptotic expansion. For instance, the $O(Ca)$ terms in (3.2)–(3.4) appear in the corresponding boundary conditions of the problem at $\unicode[STIX]{x1D716}^{2}$ . Consequently, the solvability condition of the $\unicode[STIX]{x1D716}^{2}$ problem is

(4.24) $$\begin{eqnarray}\displaystyle A_{2}(49\unicode[STIX]{x1D6FF}(2+3\unicode[STIX]{x1D702})^{2}+49Ca(4+\unicode[STIX]{x1D702})(2+3\unicode[STIX]{x1D702})+16\,320A_{2}(1+\unicode[STIX]{x1D702})^{2})=0. & & \displaystyle\end{eqnarray}$$

Equation (4.24) establishes that in the limit of a weakly deformable active droplet, two branches of steady solutions exist in the $\unicode[STIX]{x1D716}$ -neighbourhood of $Pe_{2}$ : the first branch, given by $A_{2}=0$ , corresponds to a motionless state; whereas the second, featuring $A_{2}<0$ , describes a symmetric extensional flow field akin to that of a force dipole. Recall that (3.18) connects the coefficient $A_{2}$ with a particular type of droplet deformation: $A_{2}>0$ (respectively $A_{2}<0$ ) corresponds to prolate (respectively oblate) deformations. In general, prolate and oblate deformations are not symmetric, so it is natural that our analysis yields $A_{2}$ with a particular sign.

4.3 Simultaneous onset of the two dominant instability modes ( $Pe_{1}\sim Pe_{2}$ )

The definition of the instability thresholds (3.14) implies that for

(4.25) $$\begin{eqnarray}\displaystyle Ca=4(13+12\unicode[STIX]{x1D702})/(4+\unicode[STIX]{x1D702})+O(\unicode[STIX]{x1D716}) & & \displaystyle\end{eqnarray}$$

the thresholds of the first two eigenmodes coincide, namely, $|Pe_{1}-Pe_{2}|=O(\unicode[STIX]{x1D716})$ . The purpose of this section is to investigate how the potential interaction of these two modes may impact the self-propulsion of the droplet. In this case, the leading-order flow field is given by the first two squirming modes with coefficients presented in (3.15)–(3.18), respectively. We now consider the problem at $\unicode[STIX]{x1D716}^{2}$ featuring quadratic interactions of the two instability modes shown in figure 3.

The weakly nonlinear dynamics of the system is investigated in the case of

(4.26) $$\begin{eqnarray}\displaystyle Ca=4(13+12\unicode[STIX]{x1D702})/(4+\unicode[STIX]{x1D702})+\unicode[STIX]{x1D716}Ca_{1}, & & \displaystyle\end{eqnarray}$$

where the system admits two linearly independent eigenmodes and $C^{(1)}(r,\unicode[STIX]{x1D707})$ is represented by a combination of $L_{0}(\unicode[STIX]{x1D707})$ , $L_{1}(\unicode[STIX]{x1D707})$ and $L_{2}(\unicode[STIX]{x1D707})$ . Thus, at $\unicode[STIX]{x1D716}^{2}$ quadratic nonlinearities in (4.3) produce non-zero projections onto the first five Legendre harmonics. Following the steps of the analysis presented in §§ 4.1 and 4.2, the rescaled advection–diffusion equation (3.10) reduces at $\unicode[STIX]{x1D716}^{2}$ to (4.7) with a solution given by (4.8), and the near- and far-field solutions match in the region $\unicode[STIX]{x1D716}\ll \unicode[STIX]{x1D70C}\ll 1$ , when conditions (4.9) are met.

Expanding the boundary conditions at the droplet interface to $O(\unicode[STIX]{x1D716}^{2})$ , the projection of the kinematic boundary condition onto $L_{0}(\unicode[STIX]{x1D707})$ yields

(4.27) $$\begin{eqnarray}\displaystyle \frac{24(13+12\unicode[STIX]{x1D702})}{5(2+3\unicode[STIX]{x1D702})}A_{2}^{2}=0, & & \displaystyle\end{eqnarray}$$

and can only be satisfied when $A_{2}=0$ , as for the case of mode $n=2$ alone, see (4.23). Consequently, the solvability condition for the problem at $\unicode[STIX]{x1D716}^{2}$ in the case of $Pe_{1}\sim Pe_{2}$ reads

(4.28) $$\begin{eqnarray}\displaystyle A_{1}=0. & & \displaystyle\end{eqnarray}$$

No non-trivial solution can be found within the $O(\unicode[STIX]{x1D716})$ neighbourhood of $Pe_{1}\sim Pe_{2}$ , and in particular the regime of steady self-propulsion discussed in § 4.1 ceases to exist due to the competition between the first and the second modes of the Marangoni instability. In other words, we have arrived at a conclusion that deformability may cause a qualitative change in droplet dynamics: droplets with $Ca<4(13+12\unicode[STIX]{x1D702})/(4+\unicode[STIX]{x1D702})$ exhibit a regime of steady self-propulsion with $U_{\infty }\propto Pe-Pe_{1}$ , whereas highly deformable drops with $Ca\geqslant 4(13+12\unicode[STIX]{x1D702})/(4+\unicode[STIX]{x1D702})$ have no steady regime to reach in the vicinity of the base state (2.19). As demonstrated in the following section, this qualitative change in droplet behaviour can be linked to the asymptotic disparity of the time scales associated with the first two modes of the Marangoni instability.