2.1 Physical problem

We consider a spherical active droplet of radius $R$ composed of a Newtonian fluid of viscosity $\unicode[STIX]{x1D702}_{i}$ and density $\unicode[STIX]{x1D70C}_{i}$ swimming normally to a rigid wall ($z=0$) with velocity $\boldsymbol{V}=V(t)\boldsymbol{e}_{z}$ in a second Newtonian fluid of viscosity $\unicode[STIX]{x1D702}_{o}$ and density $\unicode[STIX]{x1D70C}_{o}$. The droplet is chemically active, i.e. it exchanges a chemical solute with the external fluid, whose molecular diffusivity as well as local and background concentrations are noted respectively $D$, $C(\boldsymbol{r},t)$ and $C_{\infty }$.

In experiments, two main routes have been identified for this chemical activity (Herminghaus et al. Reference Herminghaus, Maass, Krüger, Thutupalli, Goehring and Bahr2014; Maass et al. Reference Maass, Krger, Herminghaus and Bahr2016): (i) a micellar pathway, in which empty micelles in the continuous (outer) phase capture some of the inner fluid when approaching the droplet’s boundary to form swollen micelles (Izri et al. Reference Izri, Van Der Linden, Michelin and Dauchot2014) and (ii) a molecular pathway, in which suspended surfactant molecules in the continuous phase incorporate some of the droplet’s fluid to form swollen micelles (Moerman et al. Reference Moerman, Moyses, van der Wee, Grie, van Blaaderen, Kegel, Groenewold and Brujic2017). In the following, inspired by the former mechanism, this chemical activity is generically represented as the net release of a chemical solute (e.g. swollen micelles) from the droplet’s surface into the outer fluid with a fixed and constant chemical flux per unit time and area, ${\mathcal{A}}>0$ (note that a similar model is obtained from the molecular pathway by considering a fixed capture rate of a surfactant species, see Morozov & Michelin (Reference Morozov and Michelin2019c)).

Additionally, the presence of the solute modifies the local tension of the droplet’s interface and for small enough concentration differences, $\unicode[STIX]{x1D6FE}=\unicode[STIX]{x1D6FE}_{0}+\unicode[STIX]{x1D6FE}_{1}C$ where $\unicode[STIX]{x1D6FE}_{0}$ and $\unicode[STIX]{x1D6FE}_{1}$ are two positive constants (Izri et al. Reference Izri, Van Der Linden, Michelin and Dauchot2014; Morozov & Michelin Reference Morozov and Michelin2019c), surface gradients of the solute concentration therefore lead to Marangoni stresses at the droplet’s interface.

The solute transport, which involves both molecular diffusion and advection by the self-generated Marangoni flows, plays therefore a key role in understanding the self-propulsion of a single droplet and also its interaction with a fixed solid wall (figure 1a) or a second droplet (figure 1b). In the following, and unless stated otherwise, we focus primarily on the former problem (droplet–wall interaction), and briefly analyse the latter in § 4.4. We denote by $d$ the distance at a given time between the droplet’s surface and the wall (figure 1).

2.2 Non-dimensional equations and boundary conditions

As self-propulsion arises from the Marangoni flows resulting from solute concentration gradients, a natural velocity scale for the problem is given by the drift velocity of a chemically passive droplet in an externally imposed gradient of concentration ${\mathcal{A}}/D$, i.e. $V^{\ast }={\mathcal{A}}R\unicode[STIX]{x1D6FE}_{1}/[D(2\unicode[STIX]{x1D702}_{o}+3\unicode[STIX]{x1D702}_{i})]$ (Anderson Reference Anderson1989). In the following, the problem and all quantities of interest are made dimensionless by choosing $R$, $V^{\ast }$, $R/V^{\ast }$ and ${\mathcal{A}}R/D$ as characteristic length, velocity, time and concentration scales, respectively.

Noting by $\boldsymbol{u}^{o,i}(\boldsymbol{r},t)$ the Eulerian velocity field in the outer and inner fluids measured in the fixed laboratory frame, the evolution of the dimensionless relative concentration $c=(C-C_{\infty })D/({\mathcal{A}}R)$ in the outer phase follows an advection–diffusion equation,

(2.1)$$\begin{eqnarray}Pe\left(\frac{\unicode[STIX]{x2202}c}{\unicode[STIX]{x2202}t}+\boldsymbol{u}^{o}\boldsymbol{\cdot }\unicode[STIX]{x1D735}c\right)=\unicode[STIX]{x1D6FB}^{2}c\quad \text{with }Pe=\frac{V^{\ast }R}{D}=\frac{{\mathcal{A}}R^{2}\unicode[STIX]{x1D6FE}_{1}}{(2\unicode[STIX]{x1D702}_{o}+3\unicode[STIX]{x1D702}_{i})D^{2}}.\end{eqnarray}$$

A critical parameter in our study is the Péclet number, $Pe$, which physically quantifies the relative importance of Marangoni advection and molecular diffusion in the solute transport dynamics, and can range between $O(1)$ and $O(10^{3})$ depending on the particular physical system considered in experiments (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; Suga et al. Reference Suga, Suda, Ichikawa and Kimura2018).

The droplet’s chemical activity is modelled here as a solute emission, which takes the form of a fixed chemical flux and we further assume that the rigid wall is unable to exchange any solute with the fluid. Boundary conditions for the relative concentration $c$ are therefore obtained as

(2.2a-c)$$\begin{eqnarray}\boldsymbol{n}\boldsymbol{\cdot }\left.\unicode[STIX]{x1D735}c\right|_{S}=-1,\quad \left.\boldsymbol{n}\boldsymbol{\cdot }\unicode[STIX]{x1D735}c\right|_{W}=0,\quad \left.c\right|_{r\rightarrow \infty }=0,\end{eqnarray}$$

where $S$ (respectively $W$) denotes the surface of the droplet (respectively the wall) and $\boldsymbol{n}$ the unit normal vector pointing into the outer fluid domain.

Experimentally, the typical radii of active droplets are of the order of $10$–$100~\unicode[STIX]{x03BC}\text{m}$, and their characteristic self-propulsion velocities are ${\sim}10~\unicode[STIX]{x03BC}\text{m}~\text{s}^{-1}$. At such scales, inertia is negligible (i.e. characteristic Reynolds number is $Re=\unicode[STIX]{x1D70C}_{o}V^{\ast }R/\unicode[STIX]{x1D702}_{o}\sim 10^{-3}$) and the flow velocity and pressure satisfy a steady Stokes equation in both phases,

(2.3a,b)$$\begin{eqnarray}\unicode[STIX]{x1D735}^{2}\boldsymbol{u}^{o,i}=\unicode[STIX]{x1D735}p^{o,i},\quad \unicode[STIX]{x1D735}\boldsymbol{\cdot }\boldsymbol{u}^{o,i}=0.\end{eqnarray}$$

In the laboratory frame, $\boldsymbol{u}^{o}$ vanishes far away from the droplet. At the droplet’s interface, the velocity field is continuous and Marangoni stresses result in a jump in tangential hydrodynamic stresses. Additionally, a no-slip condition is imposed at the wall surface. Noting by $\unicode[STIX]{x1D748}^{o,i}$ the Newtonian stress tensor in each fluid and $\tilde{\unicode[STIX]{x1D702}}=\unicode[STIX]{x1D702}_{i}/\unicode[STIX]{x1D702}_{o}$ the viscosity ratio, the boundary conditions for the hydrodynamic problem are therefore obtained as

(2.4a,b)$$\begin{eqnarray}\left.\boldsymbol{u}^{o}\right|_{S}=\left.\boldsymbol{u}^{i}\right|_{S},\quad \left.(\boldsymbol{I}-\boldsymbol{n}\boldsymbol{n})\boldsymbol{\cdot }(\unicode[STIX]{x1D748}^{o}-\tilde{\unicode[STIX]{x1D702}}\unicode[STIX]{x1D748}^{i})\right|_{S}\boldsymbol{\cdot }\boldsymbol{n}=-(\boldsymbol{I}-\boldsymbol{n}\boldsymbol{n})\boldsymbol{\cdot }\left(2+3\tilde{\unicode[STIX]{x1D702}}\right)\left.\unicode[STIX]{x1D735}c\right|_{S},\end{eqnarray}$$

(2.5a,b)$$\begin{eqnarray}\left.\boldsymbol{u}^{o}\right|_{W}=\mathbf{0},\quad \left.\boldsymbol{u}^{o}\right|_{r\rightarrow \infty }=\mathbf{0},\end{eqnarray}$$

and recalling that $\boldsymbol{V}$ is the droplet’s surface velocity, the impermeability condition further imposes

(2.6)$$\begin{eqnarray}\left.(\boldsymbol{u}^{o}\boldsymbol{\cdot }\boldsymbol{n})\right|_{S}=\boldsymbol{V}\boldsymbol{\cdot }\boldsymbol{n}.\end{eqnarray}$$

It should be noted here that we assume that the droplet remains spherical, since the mean surface tension $\unicode[STIX]{x1D6FE}_{0}$ is such that the capillary pressure $\unicode[STIX]{x1D6FE}_{0}/R$ is typically much larger than characteristic hydrodynamic stresses, $\unicode[STIX]{x1D702}_{o}V^{\ast }/R$ , i.e. the typical capillary number, $Ca=\unicode[STIX]{x1D702}_{o}V^{\ast }/\unicode[STIX]{x1D6FE}_{0}\approx 10^{-4}$, is much lower than unity (Izri et al. Reference Izri, Van Der Linden, Michelin and Dauchot2014). Neglecting inertia, the droplet must remain force free, which provides an additional implicit relation to determine $\boldsymbol{V}$,

(2.7)$$\begin{eqnarray}\boldsymbol{F}=\int _{S}\unicode[STIX]{x1D748}^{o}\boldsymbol{\cdot }\boldsymbol{n}\,\text{d}S=\mathbf{0}.\end{eqnarray}$$

Note that, due to the axisymmetric nature of the problem, the droplet is necessarily also torque free.

2.3 Droplet velocity and polarity

Once $c$ is known, the hydrodynamic problem is linear and is essentially similar to any surface-driven swimming problem. The total hydrodynamic force $\boldsymbol{F}$ on the droplet can therefore be decomposed into two distinct parts: (i) a Marangoni force $\boldsymbol{F}_{m}$ experienced by a fixed droplet (i.e. $\boldsymbol{V}=0$) in response to the concentration distribution at its surface, and (ii) a hydrodynamic drag $\boldsymbol{F}_{d}$ resulting from its translation and obtained from the classical Hadamard–Rybczynski problem of a translating droplet with uniform surface tension (Leal Reference Leal2007). Due to the linearity of the latter problem with respect to the droplet velocity, $\boldsymbol{F}_{d}=-\unicode[STIX]{x1D64D}\boldsymbol{\cdot }\boldsymbol{V}$, where $\unicode[STIX]{x1D64D}$, the resistance matrix, is a function only of the geometry of the problem (i.e. the droplet radius and its distance to the wall). Using $\boldsymbol{F}=\boldsymbol{F}_{d}+\boldsymbol{F}_{m}$, equation (2.7) becomes

(2.8)$$\begin{eqnarray}\displaystyle \boldsymbol{V}=\unicode[STIX]{x1D64D}^{-1}\boldsymbol{\cdot }\boldsymbol{F}_{m}. & & \displaystyle\end{eqnarray}$$

While the hydrodynamic problem for fixed $c$ is linear, the overall dynamics of the solute concentration and droplet motion is not. Indeed, advection by the self-induced Marangoni flows introduces a nonlinearity which is key to understanding the self-propulsion of an isotropic droplet or chemically active particle (Michelin et al. Reference Michelin, Lauga and Bartolo2013): when diffusion is slow enough, advection maintains a front–back polarity of the concentration gradient that drives the necessary Marangoni flows. In other words, beyond a critical Péclet number ($Pe_{c}=4$) the isotropic state of a single active droplet becomes unstable and self-propulsion can develop (Izri et al. Reference Izri, Van Der Linden, Michelin and Dauchot2014; Morozov & Michelin Reference Morozov and Michelin2019c).

It is therefore clear that the front–back asymmetry of the concentration field plays a critical role in the self-propulsion dynamics and in the following, this asymmetry is characterised by the polarity $\unicode[STIX]{x1D72B}$ of the surface concentration,

(2.9)$$\begin{eqnarray}\unicode[STIX]{x1D72B}=-\frac{1}{2\unicode[STIX]{x03C0}}\int _{S}c\boldsymbol{n}\,\text{d}S.\end{eqnarray}$$

In the case of a single droplet in an unbounded fluid, Lorentz’ reciprocal theorem for Stokes flow can be used to obtain the droplet velocity directly in terms of the surface tension (or solute concentration) gradient (Pak, Feng & Stone Reference Pak, Feng and Stone2014; Masoud & Stone Reference Masoud and Stone2019),

(2.10)$$\begin{eqnarray}\boldsymbol{V}=-\frac{1}{4\unicode[STIX]{x03C0}}\int _{S}(\boldsymbol{I}-\boldsymbol{n}\boldsymbol{n})\boldsymbol{\cdot }\unicode[STIX]{x1D735}c\,\text{d}S,\end{eqnarray}$$

which demonstrates, after integration by parts, that the droplet’s velocity and chemical polarity match exactly for a single droplet

(2.11)$$\begin{eqnarray}\boldsymbol{V}=\unicode[STIX]{x1D72B}.\end{eqnarray}$$

It should be noted that this result stems purely from hydrodynamics, and therefore does not depend on the solute dynamics (or the Péclet number $Pe$). In the axisymmetric configuration considered here, both velocity and polarity are along $\boldsymbol{e}_{z}$ and we thus focus in the following on their respective axial projections $V=\boldsymbol{V}\boldsymbol{\cdot }\boldsymbol{e}_{z}$ and $\unicode[STIX]{x1D6F1}=\unicode[STIX]{x1D72B}\boldsymbol{\cdot }\boldsymbol{e}_{z}$.

3.1 Bi-spherical coordinate system

Figure 2.

Bi-spherical coordinate system. Contours of fixed $\unicode[STIX]{x1D707}$ (solid red) and fixed $\unicode[STIX]{x1D709}$ (dashed blue) are shown at a given time $t$. The surfaces of the wall and droplet are given by $\unicode[STIX]{x1D709}=0$ and $\unicode[STIX]{x1D709}=1$, respectively.

Cartesian coordinates are well adapted to describe fluid motion or solute transport above a flat wall, yet spherical coordinates are typically more convenient to describe the flow and solute dynamics near the droplet surface. A body-fitted mesh is thus defined to describe simply both boundaries using bi-spherical coordinates, an approach that is convenient for the application of boundary conditions (figure 2). In contrast with many studies using such coordinates (Stimson & Jeffery Reference Stimson and Jeffery1926; Popescu, Tasinkevych & Dietrich Reference Popescu, Tasinkevych and Dietrich2011; Michelin & Lauga Reference Michelin and Lauga2015; Reigh & Kapral Reference Reigh and Kapral2015), the droplet is not fixed with respect to the wall so that the bi-spherical system needs to be modified at each time to match the evolving boundaries.

At a given time $t$, a point located at $(z,\unicode[STIX]{x1D70C},\unicode[STIX]{x1D719})$ in the fixed cylindrical coordinate system (with $\boldsymbol{e}_{z}$ the axis of symmetry of the problem and the origin located on the wall) has bi-spherical coordinates $(\unicode[STIX]{x1D709},\unicode[STIX]{x1D707},\unicode[STIX]{x1D719})$ defined by

(3.1a,b)$$\begin{eqnarray}\unicode[STIX]{x1D70C}=\frac{a(t)\sqrt{1-\unicode[STIX]{x1D707}^{2}}}{\unicode[STIX]{x1D6E4}},\quad z=\frac{a(t)\sinh (\unicode[STIX]{x1D706}(t)\unicode[STIX]{x1D709})}{\unicode[STIX]{x1D6E4}}\quad \text{with }\unicode[STIX]{x1D6E4}(\unicode[STIX]{x1D709},\unicode[STIX]{x1D707},t)=\cosh (\unicode[STIX]{x1D706}(t)\unicode[STIX]{x1D709})-\unicode[STIX]{x1D707},\end{eqnarray}$$

where $a(t)$ and $\unicode[STIX]{x1D706}(t)$ are functions of time to account for the time-dependent stretching of the grid. Surfaces of constant $\unicode[STIX]{x1D709}$ represent a set of non-intersecting spheres (figure 2). At any time $t$, the wall and droplet surfaces correspond to $\unicode[STIX]{x1D709}=0$ and $\unicode[STIX]{x1D709}=1$, respectively. The functions $a(t)$ and $\unicode[STIX]{x1D706}(t)$ are determined uniquely from the droplet’s radius and distance to the wall

(3.2a,b)$$\begin{eqnarray}\unicode[STIX]{x1D706}(t)=\cosh ^{-1}\left(\frac{d(t)}{R}+1\right),\quad a(t)=R\sqrt{d(t)(d(t)+2R)}.\end{eqnarray}$$

In addition, the unit vectors of the bi-spherical basis are defined as $(\boldsymbol{e}_{\unicode[STIX]{x1D709}},\boldsymbol{e}_{\unicode[STIX]{x1D707}},\boldsymbol{e}_{\unicode[STIX]{x1D719}})$ with

(3.3)$$\begin{eqnarray}\displaystyle & \displaystyle \boldsymbol{e}_{\unicode[STIX]{x1D709}}=\frac{1-\unicode[STIX]{x1D707}\cosh (\unicode[STIX]{x1D706}\unicode[STIX]{x1D709})}{\unicode[STIX]{x1D6E4}}\boldsymbol{e}_{z}-\frac{\sqrt{1-\unicode[STIX]{x1D707}^{2}}\sinh (\unicode[STIX]{x1D706}\unicode[STIX]{x1D709})}{\unicode[STIX]{x1D6E4}}\boldsymbol{e}_{\unicode[STIX]{x1D70C}}, & \displaystyle\end{eqnarray}$$

(3.4)$$\begin{eqnarray}\displaystyle & \displaystyle \boldsymbol{e}_{\unicode[STIX]{x1D707}}=\frac{\sqrt{1-\unicode[STIX]{x1D707}^{2}}\sinh (\unicode[STIX]{x1D706}\unicode[STIX]{x1D709})}{\unicode[STIX]{x1D6E4}}\boldsymbol{e}_{z}+\frac{1-\unicode[STIX]{x1D707}\cosh (\unicode[STIX]{x1D706}\unicode[STIX]{x1D709})}{\unicode[STIX]{x1D6E4}}\boldsymbol{e}_{\unicode[STIX]{x1D70C}}, & \displaystyle\end{eqnarray}$$

and the corresponding metric coefficients are

(3.5a-c)$$\begin{eqnarray}h_{\unicode[STIX]{x1D709}}=\frac{a\unicode[STIX]{x1D706}}{\unicode[STIX]{x1D6E4}},\quad h_{\unicode[STIX]{x1D707}}=\frac{a}{\unicode[STIX]{x1D6E4}\sqrt{1-\unicode[STIX]{x1D707}^{2}}},\quad h_{\unicode[STIX]{x1D719}}=\frac{a\sqrt{1-\unicode[STIX]{x1D707}^{2}}}{\unicode[STIX]{x1D6E4}}.\end{eqnarray}$$

3.2 Grid adaptation for unsteady problems

Because of the motion of the droplet (and resulting grid adaption), a point of fixed $(\unicode[STIX]{x1D709},\unicode[STIX]{x1D707},\unicode[STIX]{x1D719})$ is not fixed in the laboratory frame, i.e. it has time-dependent $(\unicode[STIX]{x1D70C},z)$-coordinates. This has consequences when solving time-dependent equations such as (2.1). Indeed, considering the local change in time of the concentration field $c$ at a fixed point ($\unicode[STIX]{x1D70C},z$) now introduces a material derivative when considering $c$ as a function of $(\unicode[STIX]{x1D709},\unicode[STIX]{x1D707},t)$, and we must thus replace

(3.6)$$\begin{eqnarray}\displaystyle \left.\frac{\unicode[STIX]{x2202}c}{\unicode[STIX]{x2202}t}\right|_{\unicode[STIX]{x1D70C},z}=\left.\frac{\unicode[STIX]{x2202}c}{\unicode[STIX]{x2202}t}\right|_{\unicode[STIX]{x1D709},\unicode[STIX]{x1D707}}-\unicode[STIX]{x1D74C}\boldsymbol{\cdot }\unicode[STIX]{x1D735}c, & & \displaystyle\end{eqnarray}$$

where $\unicode[STIX]{x1D74C}$ is the velocity in the physical space of a point with fixed $(\unicode[STIX]{x1D709},\unicode[STIX]{x1D707})$, and is obtained from (3.1). The advection–diffusion equation for $c(\unicode[STIX]{x1D709},\unicode[STIX]{x1D707},t)$ is therefore obtained as

(3.7)$$\begin{eqnarray}\left.\frac{\unicode[STIX]{x2202}c}{\unicode[STIX]{x2202}t}\right|_{\unicode[STIX]{x1D709},\unicode[STIX]{x1D707}}+\left(\boldsymbol{u}^{o}-\unicode[STIX]{x1D74C}\right)\boldsymbol{\cdot }\unicode[STIX]{x1D735}c=\frac{1}{Pe}\unicode[STIX]{x1D6FB}^{2}c,\end{eqnarray}$$

where, noting time derivatives of single-variable functions with a dot symbol,

(3.8)$$\begin{eqnarray}\unicode[STIX]{x1D74C}\boldsymbol{\cdot }\unicode[STIX]{x1D735}c=\left(\frac{\dot{\unicode[STIX]{x1D706}}\unicode[STIX]{x1D709}}{\unicode[STIX]{x1D706}}-\frac{{\dot{a}}\unicode[STIX]{x1D707}\sinh (\unicode[STIX]{x1D706}\unicode[STIX]{x1D709})}{\unicode[STIX]{x1D706}a}\right)\frac{\unicode[STIX]{x2202}c}{\unicode[STIX]{x2202}\unicode[STIX]{x1D709}}+\frac{{\dot{a}}}{a}(1-\unicode[STIX]{x1D707}^{2})\cosh (\unicode[STIX]{x1D706}\unicode[STIX]{x1D709})\frac{\unicode[STIX]{x2202}c}{\unicode[STIX]{x2202}\unicode[STIX]{x1D707}}.\end{eqnarray}$$

3.3 Hydrodynamical problem

Solving the transport equation (3.7) requires knowing the velocity field $\boldsymbol{u}$. As noted before, the hydrodynamic (Stokes) problem is instantaneous and linear and the classical method to obtain Stokes flow solutions in bi-spherical geometries can be used. The (inner and outer) flow fields are obtained in terms of streamfunctions $\unicode[STIX]{x1D713}^{i,o}$ (Michelin & Lauga Reference Michelin and Lauga2015)

(3.9a,b)$$\begin{eqnarray}u_{\unicode[STIX]{x1D709}}^{i,o}=\frac{\unicode[STIX]{x1D6E4}^{2}}{a^{2}}\frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D713}^{i,o}}{\unicode[STIX]{x2202}\unicode[STIX]{x1D707}},\quad u_{\unicode[STIX]{x1D707}}^{i,o}=-\frac{\unicode[STIX]{x1D6E4}^{2}}{a^{2}\unicode[STIX]{x1D706}\sqrt{1-\unicode[STIX]{x1D707}^{2}}}\frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D713}^{i,o}}{\unicode[STIX]{x2202}\unicode[STIX]{x1D709}},\end{eqnarray}$$

which can be written for an axisymmetric problem as (Stimson & Jeffery Reference Stimson and Jeffery1926)

(3.10)$$\begin{eqnarray}\unicode[STIX]{x1D713}^{i,o}(\unicode[STIX]{x1D709},\unicode[STIX]{x1D707},t)=\unicode[STIX]{x1D6E4}^{-3/2}\mathop{\sum }_{n=1}^{\infty }(1-\unicode[STIX]{x1D707}^{2})L_{n}^{\prime }(\unicode[STIX]{x1D707})U_{n}^{i,o}(\unicode[STIX]{x1D709},t),\end{eqnarray}$$

where $L_{n}^{\prime }$ is the first derivative of $L_{n}$ the Legendre polynomial of degree $n$, and the functions $U_{n}^{i,o}$ are given by

(3.11)$$\begin{eqnarray}\displaystyle U_{n}^{o}(\unicode[STIX]{x1D709},t) & = & \displaystyle \unicode[STIX]{x1D6FC}_{n}\cosh \left[\left(n+{\textstyle \frac{3}{2}}\right)\unicode[STIX]{x1D706}\unicode[STIX]{x1D709}\right]+\unicode[STIX]{x1D6FD}_{n}\sinh \left[\left(n+{\textstyle \frac{3}{2}}\right)\unicode[STIX]{x1D706}\unicode[STIX]{x1D709}\right]+\unicode[STIX]{x1D6FE}_{n}\cosh \left[\left(n-{\textstyle \frac{1}{2}}\right)\unicode[STIX]{x1D706}\unicode[STIX]{x1D709}\right]\nonumber\\ \displaystyle & & \displaystyle +\,\unicode[STIX]{x1D6FF}_{n}\sinh \left[\left(n-{\textstyle \frac{1}{2}}\right)\unicode[STIX]{x1D706}\unicode[STIX]{x1D709}\right],\end{eqnarray}$$

(3.12)$$\begin{eqnarray}\displaystyle U_{n}^{i}(\unicode[STIX]{x1D709},t) & = & \displaystyle \tilde{\unicode[STIX]{x1D6FC}}_{n}\text{e}^{-(n+3/2)\unicode[STIX]{x1D706}|\unicode[STIX]{x1D709}|}+\tilde{\unicode[STIX]{x1D6FD}}_{n}\text{e}^{-(n-1/2)\unicode[STIX]{x1D706}|\unicode[STIX]{x1D709}|},\end{eqnarray}$$

where $\unicode[STIX]{x1D6FC}_{n}$, $\unicode[STIX]{x1D6FD}_{n}$, $\unicode[STIX]{x1D6FE}_{n}$, $\unicode[STIX]{x1D6FF}_{n}$, $\tilde{\unicode[STIX]{x1D6FC}}_{n}$ and $\tilde{\unicode[STIX]{x1D6FD}}_{n}$ are determined at each instant $t$ independently from the kinematic and dynamic boundary conditions on the droplet and the wall surfaces. The continuity and impermeability conditions at the droplet’s boundary, equations (2.4) and (2.6), become, using both (3.3) and (3.9),

(3.13)$$\begin{eqnarray}\left.\frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D713}^{i}}{\unicode[STIX]{x2202}\unicode[STIX]{x1D707}}\right|_{\unicode[STIX]{x1D709}=1}=\left.\frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D713}^{o}}{\unicode[STIX]{x2202}\unicode[STIX]{x1D707}}\right|_{\unicode[STIX]{x1D709}=1}=\frac{a^{2}(1-\unicode[STIX]{x1D707}\cosh \unicode[STIX]{x1D706})}{(\cosh \unicode[STIX]{x1D706}-\unicode[STIX]{x1D707})^{3}}V.\end{eqnarray}$$

Integrating with respect to $\unicode[STIX]{x1D707}$ along the droplet’s boundary (and imposing $\unicode[STIX]{x1D713}^{i}=\unicode[STIX]{x1D713}^{o}=0$ on the axis of symmetry which is a streamline of the problem) one obtains

(3.14)$$\begin{eqnarray}\displaystyle \left.\unicode[STIX]{x1D713}^{i}\right|_{\unicode[STIX]{x1D709}=1}=\left.\unicode[STIX]{x1D713}^{o}\right|_{\unicode[STIX]{x1D709}=1}=\frac{(1-\unicode[STIX]{x1D707}^{2})a^{2}}{2(\cosh \unicode[STIX]{x1D706}-\unicode[STIX]{x1D707})^{2}}V. & & \displaystyle\end{eqnarray}$$

The continuity of the velocity field at the droplet’s boundary further imposes

(3.15)$$\begin{eqnarray}\left.\frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D713}^{o}}{\unicode[STIX]{x2202}\unicode[STIX]{x1D709}}\right|_{\unicode[STIX]{x1D709}=1}=\left.\frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D713}^{i}}{\unicode[STIX]{x2202}\unicode[STIX]{x1D709}}\right|_{\unicode[STIX]{x1D709}=1},\end{eqnarray}$$

and the Marangoni condition in (2.4) at the surface of the droplet becomes

(3.16)$$\begin{eqnarray}\displaystyle \left.(\unicode[STIX]{x1D70E}_{\unicode[STIX]{x1D709}\unicode[STIX]{x1D707}}^{o}-\tilde{\unicode[STIX]{x1D702}}\unicode[STIX]{x1D70E}_{\unicode[STIX]{x1D709}\unicode[STIX]{x1D707}}^{i})\right|_{\unicode[STIX]{x1D709}=1}=-\left.\frac{(2+3\tilde{\unicode[STIX]{x1D702}})(\cosh \unicode[STIX]{x1D706}-\unicode[STIX]{x1D707})\sqrt{1-\unicode[STIX]{x1D707}^{2}}}{a}\frac{\unicode[STIX]{x2202}c}{\unicode[STIX]{x2202}\unicode[STIX]{x1D707}}\right|_{\unicode[STIX]{x1D709}=1}. & & \displaystyle\end{eqnarray}$$

Finally, the no-slip boundary condition at the wall, equation (2.5), becomes

(3.17)$$\begin{eqnarray}\left.\unicode[STIX]{x1D713}^{o}\right|_{\unicode[STIX]{x1D709}=0}=0,\quad \left.\frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D713}^{o}}{\unicode[STIX]{x2202}\unicode[STIX]{x1D709}}\right|_{\unicode[STIX]{x1D709}=0}=0.\end{eqnarray}$$

Equations (3.14)–(3.17) projected in the polar direction along the $n$th Legendre polynomial (§ A.2) provide $n$ sets of 6 linear equations, that determine $\unicode[STIX]{x1D6FC}_{n}$, $\unicode[STIX]{x1D6FD}_{n}$, $\unicode[STIX]{x1D6FE}_{n}$, $\unicode[STIX]{x1D6FF}_{n}$, $\tilde{\unicode[STIX]{x1D6FC}}_{n}$ and $\tilde{\unicode[STIX]{x1D6FD}}_{n}$ (and thus the streamfunction) uniquely in terms of the surface concentration distribution.

3.4 Transport problem

Solving the transport equation, equation (3.7), also exploits a spectral decomposition of $c$ along the Legendre polynomials in the polar direction. Inspired by the separated form of the solution for Laplace’s equation in bi-spherical coordinates, the relative concentration field $c$ (which vanishes at infinity here) is thus decomposed as

(3.18)$$\begin{eqnarray}\displaystyle c(\unicode[STIX]{x1D709},\unicode[STIX]{x1D707},t)=\unicode[STIX]{x1D6E4}^{1/2}\mathop{\sum }_{n=0}^{\infty }c_{n}(\unicode[STIX]{x1D709},t)L_{n}(\unicode[STIX]{x1D707}), & & \displaystyle\end{eqnarray}$$

where the $c_{n}(\unicode[STIX]{x1D709},t)$ functions are yet to be determined. Using (3.5) and (3.9), the advection–diffusion equation, equation (3.7), becomes

(3.19)$$\begin{eqnarray}\displaystyle & & \displaystyle \left.\frac{\unicode[STIX]{x2202}c}{\unicode[STIX]{x2202}t}\right|_{\unicode[STIX]{x1D709},\unicode[STIX]{x1D707}}+\!\left[\frac{{\dot{a}}\unicode[STIX]{x1D707}\sinh (\unicode[STIX]{x1D706}\unicode[STIX]{x1D709})}{\unicode[STIX]{x1D706}a}-\frac{\dot{\unicode[STIX]{x1D706}}\unicode[STIX]{x1D709}}{\unicode[STIX]{x1D706}}+\frac{\unicode[STIX]{x1D6E4}^{3}}{\unicode[STIX]{x1D706}a^{3}}\frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D713}^{o}}{\unicode[STIX]{x2202}\unicode[STIX]{x1D707}}\right]\frac{\unicode[STIX]{x2202}c}{\unicode[STIX]{x2202}\unicode[STIX]{x1D709}}-\!\left[\frac{{\dot{a}}}{a}(1-\unicode[STIX]{x1D707}^{2})\cosh (\unicode[STIX]{x1D706}\unicode[STIX]{x1D709})+\frac{\unicode[STIX]{x1D6E4}^{3}}{\unicode[STIX]{x1D706}a^{3}}\frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D713}^{o}}{\unicode[STIX]{x2202}\unicode[STIX]{x1D709}}\right]\!\frac{\unicode[STIX]{x2202}c}{\unicode[STIX]{x2202}\unicode[STIX]{x1D707}}\nonumber\\ \displaystyle & & \displaystyle \qquad =\frac{1}{Pe}\frac{\unicode[STIX]{x1D6E4}^{3}}{\unicode[STIX]{x1D706}a^{2}}\left[\frac{1}{\unicode[STIX]{x1D706}}\frac{\unicode[STIX]{x2202}}{\unicode[STIX]{x2202}\unicode[STIX]{x1D709}}\left(\frac{1}{\unicode[STIX]{x1D6E4}}\frac{\unicode[STIX]{x2202}c}{\unicode[STIX]{x2202}\unicode[STIX]{x1D709}}\right)+\frac{\unicode[STIX]{x2202}}{\unicode[STIX]{x2202}\unicode[STIX]{x1D707}}\left(\frac{(1-\unicode[STIX]{x1D707}^{2})}{\unicode[STIX]{x1D6E4}}\frac{\unicode[STIX]{x2202}c}{\unicode[STIX]{x2202}\unicode[STIX]{x1D707}}\right)\right],\end{eqnarray}$$

and substituting $\unicode[STIX]{x1D713}^{o}$ and $c$ from (3.10) and (3.18) yields

(3.20)$$\begin{eqnarray}\displaystyle & & \displaystyle \mathop{\sum }_{n=0}^{\infty }\left\{\frac{L_{n}}{\unicode[STIX]{x1D6E4}^{1/2}}\frac{\unicode[STIX]{x2202}c_{n}}{\unicode[STIX]{x2202}t}+\frac{{\dot{a}}}{a}\left[\frac{(1+\unicode[STIX]{x1D707}\cosh (\unicode[STIX]{x1D706}\unicode[STIX]{x1D709}))L_{n}-2\cosh (\unicode[STIX]{x1D706}\unicode[STIX]{x1D709})(1-\unicode[STIX]{x1D707}^{2})L_{n}^{\prime }}{2\unicode[STIX]{x1D6E4}^{1/2}}\right]c_{n}\right.\nonumber\\ \displaystyle & & \displaystyle \qquad +\,\left(\frac{{\dot{a}}\unicode[STIX]{x1D707}\sinh (\unicode[STIX]{x1D706}\unicode[STIX]{x1D709})-\dot{\unicode[STIX]{x1D706}}\unicode[STIX]{x1D709}a}{\unicode[STIX]{x1D706}a}\right)\frac{L_{n}}{\unicode[STIX]{x1D6E4}^{1/2}}\frac{\unicode[STIX]{x2202}c_{n}}{\unicode[STIX]{x2202}\unicode[STIX]{x1D709}}\nonumber\\ \displaystyle & & \displaystyle \qquad +\,\frac{1}{\unicode[STIX]{x1D706}a^{3}}\mathop{\sum }_{k=1}^{\infty }\left[\left(\frac{3}{2}(1-\unicode[STIX]{x1D707}^{2})L_{k}^{\prime }L_{n}-k(k+1)\unicode[STIX]{x1D6E4}L_{k}L_{n}\right)U_{k}\frac{\unicode[STIX]{x2202}c_{n}}{\unicode[STIX]{x2202}\unicode[STIX]{x1D709}}\right.\nonumber\\ \displaystyle & & \displaystyle \qquad +\,\frac{\unicode[STIX]{x1D706}\sinh (\unicode[STIX]{x1D706}\unicode[STIX]{x1D709})}{2}\left[3(1-\unicode[STIX]{x1D707}^{2})L_{k}^{\prime }L_{n}^{\prime }-k(k+1)L_{n}L_{k}\right]U_{k}c_{n}\nonumber\\ \displaystyle & & \displaystyle \qquad \left.\left.+\,(1-\unicode[STIX]{x1D707}^{2})L_{k}^{\prime }\left(\frac{L_{n}}{2}-\unicode[STIX]{x1D6E4}L_{n}^{\prime }\right)\frac{\unicode[STIX]{x2202}U_{k}}{\unicode[STIX]{x2202}\unicode[STIX]{x1D709}}c_{n}\right]\right\}\nonumber\\ \displaystyle & & \displaystyle \quad =\frac{\unicode[STIX]{x1D6E4}^{3/2}}{a^{2}Pe}\mathop{\sum }_{n=0}^{\infty }\left(\frac{1}{\unicode[STIX]{x1D706}^{2}}\frac{\unicode[STIX]{x2202}^{2}c_{n}}{\unicode[STIX]{x2202}\unicode[STIX]{x1D709}^{2}}-\left(n+\frac{1}{2}\right)^{2}c_{n}\right)L_{n}.\end{eqnarray}$$

Projecting the advection–diffusion equation, equation (3.20), onto $L_{p}(\unicode[STIX]{x1D707})$ provides a set of coupled partial differential equations for $\boldsymbol{C}(\unicode[STIX]{x1D709},t)=[c_{0}(\unicode[STIX]{x1D709},t),c_{1}(\unicode[STIX]{x1D709},t),\ldots ,c_{N}(\unicode[STIX]{x1D709},t)]$ (see appendix A) which can be formally written as

(3.21)$$\begin{eqnarray}\displaystyle & & \displaystyle \boldsymbol{H}\boldsymbol{\cdot }\frac{\unicode[STIX]{x2202}\boldsymbol{C}}{\unicode[STIX]{x2202}t}+\left(\boldsymbol{B}^{1}\boldsymbol{\cdot }\boldsymbol{U}+\boldsymbol{B}^{2}\boldsymbol{\cdot }\frac{\unicode[STIX]{x2202}\boldsymbol{U}}{\unicode[STIX]{x2202}\unicode[STIX]{x1D709}}+\boldsymbol{G}^{1}\right)\boldsymbol{\cdot }\boldsymbol{C}\nonumber\\ \displaystyle & & \displaystyle \quad +\,\left(\boldsymbol{B}^{3}\boldsymbol{\cdot }\boldsymbol{U}+\boldsymbol{G}^{2}\right)\boldsymbol{\cdot }\frac{\unicode[STIX]{x2202}\boldsymbol{C}}{\unicode[STIX]{x2202}\unicode[STIX]{x1D709}}=\frac{1}{\text{Pe}}\left(\boldsymbol{A}^{1}\boldsymbol{\cdot }\boldsymbol{C}+\boldsymbol{A}^{2}\boldsymbol{\cdot }\frac{\unicode[STIX]{x2202}^{2}\boldsymbol{C}}{\unicode[STIX]{x2202}\unicode[STIX]{x1D709}^{2}}\right),\end{eqnarray}$$

where the second-order tensors $\boldsymbol{H}$, $\boldsymbol{A}^{i}$, $\boldsymbol{G}^{i}$ ($i=1,2$) and the third-order tensor $\boldsymbol{B}^{j}$ ($j=1,2,3$) have coefficients that are obtained in terms of integrals of appropriate combinations of Legendre polynomials and depend on $\unicode[STIX]{x1D709}$ (see appendix A). Physically, terms in $\boldsymbol{H}$, $\boldsymbol{G}^{i}$, $\boldsymbol{B}^{j}$ and $\boldsymbol{A}^{i}$ are related to the local time derivative (for fixed $\unicode[STIX]{x1D709}$ and $\unicode[STIX]{x1D707}$), the grid adaptation and convection by the Marangoni flow and diffusion, respectively. In (3.21), $\boldsymbol{U}(\unicode[STIX]{x1D709},t)=[U_{1}^{o}(\unicode[STIX]{x1D709},t),U_{2}^{o}(\unicode[STIX]{x1D709},t),\ldots ,U_{N}^{o}(\unicode[STIX]{x1D709},t)]$ is a linear and instantaneous function of $\boldsymbol{C}$ (see § 3.3).

Finally, this set of evolution equations for $c_{n}(\unicode[STIX]{x1D709},t)$, equation (3.21), which are second-order in space, must be complemented by appropriate boundary conditions. After substitution of (3.18), and projection onto the $n$th Legendre polynomial the flux conditions, equations (2.2), at the droplet ($\unicode[STIX]{x1D709}=1$) and wall ($\unicode[STIX]{x1D709}=0$) surfaces become

(3.22)$$\begin{eqnarray}\displaystyle & \hspace{-24.0pt}\displaystyle \left.\left(\frac{\unicode[STIX]{x1D706}\sinh \unicode[STIX]{x1D706}}{2}\,c_{n}+\cosh \unicode[STIX]{x1D706}\frac{\unicode[STIX]{x2202}c_{n}}{\unicode[STIX]{x2202}\unicode[STIX]{x1D709}}-\frac{n+1}{2n+3}\frac{\unicode[STIX]{x2202}c_{n+1}}{\unicode[STIX]{x2202}\unicode[STIX]{x1D709}}-\frac{n}{2n-1}\frac{\unicode[STIX]{x2202}c_{n-1}}{\unicode[STIX]{x2202}\unicode[STIX]{x1D709}}\right)\right|_{\unicode[STIX]{x1D709}=1} & \displaystyle =\sqrt{2}a\unicode[STIX]{x1D706}\text{e}^{-(n+1/2)|\unicode[STIX]{x1D706}|},\nonumber\\ \displaystyle & & \displaystyle\end{eqnarray}$$

(3.23)$$\begin{eqnarray}\displaystyle & \hspace{-24.0pt}\displaystyle \left.\left(\frac{\unicode[STIX]{x2202}c_{n}}{\unicode[STIX]{x2202}\unicode[STIX]{x1D709}}-\frac{n+1}{2n+3}\frac{\unicode[STIX]{x2202}c_{n+1}}{\unicode[STIX]{x2202}\unicode[STIX]{x1D709}}-\frac{n}{2n-1}\frac{\unicode[STIX]{x2202}c_{n-1}}{\unicode[STIX]{x2202}\unicode[STIX]{x1D709}}\right)\right|_{\unicode[STIX]{x1D709}=0}=0. & \displaystyle\end{eqnarray}$$

At each time step, the functions $(U_{n}^{o})_{n}$ are obtained from $(c_{n})_{n}$ following the results of § 3.3. Note that only the outer flow needs to be known explicitly. The set of nonlinear partial differential equations, equation (3.21), is then solved using finite differences and a uniform grid of $N_{\unicode[STIX]{x1D709}}$ points with $0\leqslant \unicode[STIX]{x1D709}\leqslant 1$. Advective terms, i.e. those involving $U_{n}^{o}(\unicode[STIX]{x1D709})$, are treated explicitly while a Crank–Nicholson scheme is used to account for the diffusive terms.

It should be noted here that the present framework can easily be adapted to consider the collision of two droplets rather than a droplet with a wall: for two symmetric droplets, the $\unicode[STIX]{x1D709}$-grid is simply extended to $-1\leqslant \unicode[STIX]{x1D709}\leqslant 1$ with the boundary conditions at the wall ($\unicode[STIX]{x1D709}=0$) being replaced by appropriate conditions on the second droplet ($\unicode[STIX]{x1D709}=-1$). Alternatively, for a purely symmetric situation (see § 4.4), we only need to solve for the right half-plane and the right-most droplet by imposing the symmetry of the concentration and velocity fields on $\unicode[STIX]{x1D709}=0$.

3.5 Validation of the numerical model

The validity and accuracy of the present approach are tested by comparison with the case of a single self-propelled droplet considered by Izri et al. (Reference Izri, Van Der Linden, Michelin and Dauchot2014). To this end, we consider the case of two droplets initially separated by a distance $d=48$, i.e. far enough that one expects their interaction to be only weak and to recover the single-droplet results. The concentration field is initialised using the purely diffusive ($Pe=0$) solution for which an analytical solution is available for $c$ (Michelin & Lauga Reference Michelin and Lauga2015). At $t>0$, the previous simulation framework is used for a fixed non-zero value of $Pe$; both droplets are initially forced to move at a fixed positive velocity $V=0.1$ until $t=2$ and are allowed to evolve force free for $t>2$; after a transient regime, their velocity relaxes toward a fixed and common value identified as their self-propulsion velocity $V_{0}$ when isolated. The mean long-time velocity of the droplets is measured at $t=3000$; this is reported in figure 3 and compared to the results of Izri et al. (Reference Izri, Van Der Linden, Michelin and Dauchot2014) for a single droplet. The results are in excellent agreement and validate the present framework: the maximum relative errors obtained are approximately $2\%$ for $6\leqslant Pe\leqslant 20$, when using $N=60$ polar modes and $N_{\unicode[STIX]{x1D709}}=100$ regularly spaced grid points in $\unicode[STIX]{x1D709}$, a resolution precise enough to guarantee an accurate description of the physical processes, yet light enough to analyse the time-dependent rebound dynamics of the droplet as well as the influence of $Pe$ on the detailed chemo-hydrodynamic interaction between the droplet and the wall.

Figure 3.

(a) Mean self-propulsion velocity $V_{0}$ of two active droplets with $\tilde{\unicode[STIX]{x1D702}}=1/36$ initially located at a distance $d=48$ from each other (red crosses). The results of Izri et al. (Reference Izri, Van Der Linden, Michelin and Dauchot2014) for a single active droplet with the same viscosity ratio are also reported (black solid line). (b) The concentration distribution around the self-propelling droplets is shown for $Pe=6$ and $Pe=20$.

4.1 Collision dynamics for moderate $Pe$

We first analyse the collision dynamics for moderate advection, i.e. for $Pe$ slightly above the instability threshold ($Pe=6$ is chosen here). Figure 4 displays the evolution of the axial velocity $V$ and polarity $\unicode[STIX]{x1D6F1}$ of the concentration field with the distance, $d$, between the droplet and the wall. A corresponding movie (Movie 1) is also provided as Supplementary Material available at https://doi.org/10.1017/jfm.2019.1055. Initially, the droplet swims toward the wall with a constant self-propulsion velocity $-V_{0}$, where $V_{0}(Pe)$ is the magnitude of self-propulsion of an isolated droplet. As it approaches the wall, it decelerates and reverses direction ($V=0$) at a finite distance $d_{rb}(Pe)$ from the wall ($d_{rb}=1.4$ for $Pe=6$). In a second phase, the droplet accelerates away from the wall and eventually reaches again its self-propulsion velocity $V_{0}$. A main observation of figure 4 is that the velocity and polarity are almost equal throughout the collision and the consequence of this is discussed in more depth below. In the following, we analyse each sequence of the interaction to identify the roles of the different mechanisms.

Figure 4.

Collision of an active droplet with a rigid wall at $Pe=6$. The evolution of the droplet velocity $V$ (blue solid) and polarity $\unicode[STIX]{x1D6F1}$ (red dashed) during the collision are reported in terms of the distance $d$ between the droplet’s surface and the wall (b). Snapshots of the concentration (a) and streamfunction (c) are also shown for four representative positions indicated as (A–D) on panel (b).

4.1.2 Near field and re-acceleration toward self-propulsion

When the droplet velocity vanishes (instant B in figure 4), the polarity of the concentration field also comes close to zero. A closer look at the distribution of solute on the surface at that instant in fact reveals that the concentration distribution is almost homogeneous (its variance is reduced by an order of magnitude, when compared with the initial self-propelling state): this results in the droplet’s arrest as there is no longer a Marangoni effect acting on the surface of the droplet and the fluid is at rest (see figure 4).

Figure 5.

Evolution in time of the droplet’s velocity for $Pe=6$ (a) and $Pe=20$ (b) during its re-acceleration after its collision with the wall (solid blue). This evolution is compared to the acceleration of an isolated droplet initially forced with a positive velocity $V=0.1$ before being released force free at $t=0$ (dashed red).

However, this equilibrium is only ephemeral as the presence of the wall promptly breaks this uniform distribution: the chemical flux at the droplet’s surface being spatially uniform, the confinement on the side of the wall leads to an increased solute concentration there and a repelling Marangoni force. As a result, the droplet drifts away from the wall (figure 4, instant C). Because $Pe$ is greater than the instability threshold, this perturbation of the concentration field simultaneously leads to the development of the same instability phenomenon that conferred the droplet its initial velocity, until it reaches $V_{0}$ as the droplet moves far away from the wall.

In order to study the influence of the wall in the droplet ‘forced’ re-acceleration, figure 5(a) compares this second phase of the motion with the situation of a single isolated droplet initially pushed at a finite velocity in the positive $z$ direction before being left force free. The acceleration is initially greater in the second part of the droplet collision with the wall than in the reference isolated droplet case. Indeed, the presence of the wall reinforces the droplet acceleration (by accumulating more solute at its back) than in the case where the solute is able to diffuse freely.

4.2 Collision at higher $Pe$

The results obtained for moderate $Pe$ presented a rather simple picture of the collision dynamics: dominated by the chemical interactions with the wall, it amounts to a slowing down and repulsion of the droplet under the effect of the accumulating chemical solute in front of it due to the confining presence of the wall, the wall’s hydrodynamic effect being mostly subdominant. The picture becomes however much more complex as the importance of convection of solute vs. diffusion is increased, and the focus of this section is to analyse how wall interactions and collision dynamics are modified as $Pe$ is increased.

Figure 6.

Collision of an active droplet with a rigid wall at $Pe$$=20$. The evolution of the droplet velocity $V$ (blue solid) and polarity $\unicode[STIX]{x1D6F1}$ (red dashed) during the collision are reported in terms of the distance $d$ between the droplet’s surface and the wall (b). Snapshots of the concentration (a) and streamfunction (c) are also shown for four representative positions indicated as (A–D) on panel (b).

Figure 6 presents the evolution of velocity and polarity throughout the collision for $Pe=20$. A corresponding movie (Movie 2) is also provided as Supplementary Material. As for $Pe=6$, the droplet initially propels at $-V_{0}$ toward the wall (A) and decelerates up to a stopping point (B). This minimum rebound distance, $d_{rb}$, is however much lower ($d_{rb}=0.3$ for $Pe=20$). Also, unlike for moderate $Pe$, the re-acceleration of the droplet is not a smooth process. In particular, it displays a clear velocity plateau right after the rebound during which the droplet velocity remains almost constant (C). Eventually, and as expected, the droplet reaches once again its self-propulsion velocity as it moves away from the wall, whose influence becomes negligible (D). In contrast with the moderate-$Pe$ collision (figure 4), we note that the polarity and velocity do not match one another anymore during most of the near-field interactions with the wall, suggesting a stronger hydrodynamic influence of the confinement.

Figure 7.

Evolution as a function of the droplet–wall distance $d$ of the velocity $V$ (solid blue), the polarity $\unicode[STIX]{x1D6F1}$ (red dashed) and front–back surface concentration contrast $\unicode[STIX]{x0394}c=c_{front}-c_{back}$ (black dashed-dotted) at $Pe=20$. The inset shows the evolution of $c_{back}$ and $c_{front}$ with $d$ individually and the definition of these two quantities is given at the top of the figure.

4.3 Rebound distance for varying $Pe$

Figure 8.

Evolution of the rebound distance $d_{rb}$ (measured when $V=0$) for the droplet–wall collision (red crosses) and the two-droplet collision (blue stars, see § 4.4). The inset reproduces the same data in log–log scales to identify the asymptotic scaling for near-critical and large $Pe$. The concentration fields around an active droplet approaching the wall (i.e. propelling to the left) are also shown for three representative Péclet numbers: $Pe=6$, 13 and 20.

In the previous sections, the droplet–wall interaction was analysed in detail for two different values of $Pe$. One key feature was that, due to the structure of the concentration field ahead of the moving droplet, the rebound distance (i.e. the minimum distance of approach of the droplet to the wall) is reduced when advection plays a more important role in the solute transport. The goal of this section is to provide a more complete characterisation of this phenomenon and we thus now focus on the evolution with $Pe$ of the distance $d_{rb}$ between the wall and the front of the droplet at the time it reverses direction (figure 8). A monotonic decrease of $d_{rb}$ is observed, as expected from the structure of the concentration field ahead of the moving droplet (see figure 8a).

We also note the existence of two distinct regimes in this decrease. For moderate $Pe$, and as $Pe$ approaches the minimum value for self-propulsion ($Pe_{c}=4$), $d_{rb}$ diverges as $(Pe-Pe_{c})^{-1}$: close to the self-propulsion threshold, the droplet is more sensitive to the wall’s influence and is repelled at much greater distances. The asymptotic analysis of the collision near $Pe_{c}$, presented in § 5, confirms this scaling and provides more insight into the interaction and rebound dynamics. For larger $Pe$ (typically $Pe\gtrsim 10$), a slower decrease is observed as $(Pe-Pe_{c})^{-1/2}$.

4.4 Comparison to the two-droplet collision

Figure 9.

Collision between two active droplets at $Pe=6$ (a) and $Pe=20$ (b). In each case, the evolution with $d$, the half-distance between the two droplets, of the droplet velocity $V$ (solid blue) and polarity $\unicode[STIX]{x1D6F1}$ (dashed red) is shown. For reference the corresponding results for the collision of a single droplet with a no-slip wall are also shown for the velocity (solid black) and polarity (dashed black).

The interaction of a droplet with a chemically inert wall shares many similarities to the symmetric interaction of two identical droplets. Mathematically, the chemical problem is in fact identical, the no-flux condition at the wall in (2.2) being strictly equivalent to a perfect symmetry of the concentration field as for the case of two symmetric droplets. The only difference therefore lies in the hydrodynamic problem and resulting flow field: a no-slip condition is applied for the case of rebound on a wall, equation (2.5), while symmetry conditions on the velocity field would hold on $z=0$ for the case of two droplets (effectively amounting to the presence of a free surface rather than a rigid wall). The modification of the hydrodynamic field may nevertheless have significant consequences, in particular at larger $Pe$ due to the importance of advection by this flow field in the solute dynamics.

Before closing this section, we therefore briefly analyse the difference between the two situations in more detail. Defining $2d$ as the minimum distance between the surfaces of the two droplets (figure 1b), and using the same numerical approach as for the droplet–wall interaction, we report in figure 9 the evolution with $d$ of the velocity and polarity of the right-hand droplet for moderate and higher $Pe$, and compare these results to those obtained for the wall collision. Both configurations lead to the same dynamics except in the immediate vicinity of the stopping point where some small variations can be identified. These differences are more pronounced for larger $Pe$, which is likely due to the closer proximity to the wall. Indeed, the no-slip condition applied at the wall surface is relaxed for the two-droplet collision and therefore the hydrodynamical interactions close to the wall are not identical (see Kim & Karrila Reference Kim and Karrila2013, ch. 12).

Figure 8 further compares the evolution of $d_{rb}$ with $Pe$ for the droplet–wall and droplet–droplet collisions. Once again the evolutions are similar, although the rebound distance is systematically smaller in the case of two droplets (or of a free surface) than for a droplet and a no-slip wall. This may be due to the different structure of the flow field between the droplet and wall/symmetry plane: a no-slip condition on the wall prevents the fluid flowing along the $z=0$ axis and it thus cannot convect away as much solute as in the case of two droplets.

5.1 Problems at $\unicode[STIX]{x1D716}^{0}$ and $\unicode[STIX]{x1D716}^{1}$: isolated drops

To compute the terms $N_{i}^{(0)}$, $N_{i}^{(1)}$ and $F_{i}^{(1)}$ of expansions (5.3)–(5.4), we consider each of the droplets separately. That is, we substitute (5.3)–(5.4) into the problem formulated by (2.1)–(2.7) in the case of a solitary spherical active drop. We then collect $O(1)$ terms and recover the isotropic solution,

(5.5)$$\begin{eqnarray}N_{i}^{(0)}(\boldsymbol{r}_{i})=1/r_{i},\end{eqnarray}$$

that corresponds to a motionless droplet and quiescent fluid.

We proceed with the solution of the problem formulated by the $O(\unicode[STIX]{x1D716})$ terms. Since the flow field does not vanish at $\unicode[STIX]{x1D716}^{1}$, we need to write a solution to Stokes equations (2.3) within and around a spherical drop. In contrast to the advection–diffusion equation that requires a composite solution, equations (5.3)–(5.4), Stokes equations are linear and admit a solution that is uniformly valid. At each order in $\unicode[STIX]{x1D716}$, the axisymmetric solution of the Stokes problem, equation (2.3), within and outside of the spherical droplet is given by a superposition of orthogonal modes (Lamb Reference Lamb1945; Happel & Brenner Reference Happel and Brenner1983; Leal Reference Leal2007),

(5.6)$$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D713}_{n}^{i}(\boldsymbol{r})=r^{n+1}\left(1-r^{2}\right)\left(1-\unicode[STIX]{x1D707}^{2}\right)L_{n}^{\prime }(\unicode[STIX]{x1D707}), & \displaystyle\end{eqnarray}$$

(5.7)$$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D713}_{n}^{o}(\boldsymbol{r})=\left\{\begin{array}{@{}ll@{}}\left(1-r^{3}\right)\left(1-\unicode[STIX]{x1D707}^{2}\right)/r\quad \quad & n=1,\\ \left(1-r^{2}\right)\left(1-\unicode[STIX]{x1D707}^{2}\right)L_{n}^{\prime }(\unicode[STIX]{x1D707})/r^{n}\quad \quad & n>1,\end{array}\right. & \displaystyle\end{eqnarray}$$

where $\unicode[STIX]{x1D713}_{n}^{i,o}$ denote the streamfunctions corresponding to the $n$th mode of the flow decomposition within and outside of the drop, respectively.

Equation (5.7) implies that in the reference frame co-moving with the drop, the far-field flow at $O(\unicode[STIX]{x1D716})$ is unidirectional and the far-field advection–diffusion equation for each droplet may thus be rewritten as

(5.8)$$\begin{eqnarray}-Pe\boldsymbol{V}_{i}^{(1)}\boldsymbol{\cdot }\unicode[STIX]{x1D735}_{\unicode[STIX]{x1D746}}F_{i}^{(1)}(\unicode[STIX]{x1D746}_{i})=\unicode[STIX]{x1D6FB}_{\unicode[STIX]{x1D746}}^{2}F_{i}^{(1)}(\unicode[STIX]{x1D746}_{i}),\end{eqnarray}$$

where $\unicode[STIX]{x1D735}_{\unicode[STIX]{x1D746}}$ denotes the gradient in stretched coordinates. Since we disregard the droplet interactions at this order, we may adopt the corresponding solutions for $N_{i}^{(1)}$, and $F_{i}^{(1)}$ obtained by Morozov & Michelin (Reference Morozov and Michelin2019c,Reference Morozov and Michelina),

(5.9a,b)$$\begin{eqnarray}N_{i}^{(1)}(\boldsymbol{r}_{i})=-2V_{i}^{(1)}\left(1+\unicode[STIX]{x1D707}_{i}+\unicode[STIX]{x1D707}_{i}\frac{2-3r_{i}}{4r_{i}^{3}}\right),\quad F_{i}^{(1)}(\unicode[STIX]{x1D746}_{i})=\frac{\text{e}^{-2V_{i}^{(1)}\unicode[STIX]{x1D70C}_{i}(1+\unicode[STIX]{x1D707}_{i})}}{\unicode[STIX]{x1D70C}_{i}},\end{eqnarray}$$

where $V_{i}^{(1)}=\big|\boldsymbol{V}_{i}^{(1)}\big|$ and $\unicode[STIX]{x1D707}_{i}\equiv \cos \unicode[STIX]{x1D703}_{i}$. Note this effectively implies that $V_{i}^{(1)}>0$ and that we define $\unicode[STIX]{x1D703}_{i}$ (and therefore $\unicode[STIX]{x1D707}_{i}$) from an axis of reference that is oriented along the direction of propulsion of each droplet, i.e. toward (respectively away from) the second droplet for the approaching (respectively departing) case (figure 10).

Also note that (5.3), (5.7) and (5.5) indicate that in the laboratory frame the leading-order flow at a distance $1/\unicode[STIX]{x1D716}$ from an active droplet is $O(\unicode[STIX]{x1D716}^{4})$ (recall that at $\unicode[STIX]{x1D716}^{0}$ the fluid is motionless), while the droplet’s chemical footprint is $O(\unicode[STIX]{x1D716})$. That is, at a distance $1/\unicode[STIX]{x1D716}$ from both droplets, the advection—diffusion equation, equation (2.1), reduces to an unsteady diffusion equation in the laboratory frame,

(5.10)$$\begin{eqnarray}Pe\frac{\unicode[STIX]{x2202}F}{\unicode[STIX]{x2202}t}=\unicode[STIX]{x1D6FB}^{2}F+O(\unicode[STIX]{x1D716}^{4}),\end{eqnarray}$$

where $F$ denotes the combined far-field concentration footprint of a pair of active drops, the time derivative accounts for the displacement of the droplets in the laboratory frame and $O(\unicode[STIX]{x1D716}^{4})$ corresponds to the contribution of the flow field. Equation (5.10) is linear, therefore $F$ can be found as a superposition of the contributions from the individual droplets,

(5.11)$$\begin{eqnarray}F=F_{1}+F_{2}=\unicode[STIX]{x1D716}(F_{1}^{(1)}+F_{2}^{(1)})+\unicode[STIX]{x1D716}^{2}(F_{1}^{(2)}+F_{2}^{(2)})+\cdots\end{eqnarray}$$

Finally, we demonstrate that interaction between the droplets appears only in the problem at $O(\unicode[STIX]{x1D716}^{2})$. To this end, we write the concentration field of droplet 2 in the coordinate system of droplet 1 and expand the result in powers of $\unicode[STIX]{x1D716}$ in the cases of approaching and departing drops, respectively,

(5.12)$$\begin{eqnarray}\displaystyle & \displaystyle F_{2,approach}=\unicode[STIX]{x1D716}\frac{\text{e}^{-8DV^{(1)}}}{2D}-\unicode[STIX]{x1D716}^{2}\frac{\text{e}^{-8DV^{(1)}}(1+8DV^{(1)})}{4D^{2}}r_{1}\unicode[STIX]{x1D707}_{1}+O(\unicode[STIX]{x1D716}^{3}), & \displaystyle\end{eqnarray}$$

(5.13)$$\begin{eqnarray}\displaystyle & \displaystyle F_{2,departure}=\unicode[STIX]{x1D716}\frac{1}{2D}+\unicode[STIX]{x1D716}^{2}\frac{r_{1}\unicode[STIX]{x1D707}_{1}}{4D^{2}}+O(\unicode[STIX]{x1D716}^{3}), & \displaystyle\end{eqnarray}$$

where $V^{(1)}\equiv V_{1}^{(1)}=V_{2}^{(1)}$. Note that the first term on the right-hand side of (5.12)–(5.13) is constant and, thus, cannot produce any Marangoni stresses. In turn, the second term on the right-hand side of (5.12) constitutes a unidirectional concentration gradient that implements the effect of droplet 2 on droplet 1. This term is quadratic in $\unicode[STIX]{x1D716}$, which justifies that the droplets may be considered separately at $O(\unicode[STIX]{x1D716}^{0})$ and $O(\unicode[STIX]{x1D716}^{1})$.

5.2 Problem at $\unicode[STIX]{x1D716}^{2}$: droplet interaction

Section 5.1 established that the interaction of a pair of distant droplets is implemented by a linear concentration gradient in (5.12). On the other hand, in the reference frame co-moving with the drop, the near-field advection–diffusion, equation (2.1), reduces to an inhomogeneous steady diffusion equation (Morozov & Michelin Reference Morozov and Michelin2019a,Reference Morozov and Michelinc) that admits (5.12) as a solution. Therefore, the influence of the concentration gradient imposed by droplet 2 in the vicinity of droplet 1 only appears in the boundary conditions for droplet 1 and vice versa. As a result, general solutions for $N_{i}^{(2)}$ and $F_{i}^{(2)}$ may be adapted directly from Morozov & Michelin (Reference Morozov and Michelin2019a,Reference Morozov and Michelinc),

(5.14)$$\begin{eqnarray}\displaystyle & & \displaystyle \hspace{-36.0pt}N_{i}^{(2)}(\boldsymbol{r}_{i})=\frac{C_{0,i}}{r_{i}}+\frac{(V^{(1)})^{2}}{30r_{i}^{5}}(8-15r_{i}+20r_{i}^{3}+80r_{i}^{6})-\frac{1}{2}(\unicode[STIX]{x1D6FF}V^{(1)}+8A_{1,i})\nonumber\\ \displaystyle & & \displaystyle \hspace{-36.0pt}\qquad +\,L_{1}(\unicode[STIX]{x1D707}_{i})\left(\frac{C_{1,i}}{r_{i}^{2}}+4r_{i}(V^{(1)})^{2}-\frac{1+2r_{i}^{3}}{4r_{i}^{3}}(\unicode[STIX]{x1D6FF}V^{(1)}+8A_{1,i})\right)\nonumber\\ \displaystyle & & \displaystyle \hspace{-36.0pt}\qquad +\,L_{2}(\unicode[STIX]{x1D707}_{i})\left(\frac{C_{2,i}}{r_{i}^{3}}+\frac{(V^{(1)})^{2}}{21r_{i}^{5}}(10-21r_{i}+70r_{i}^{3}-42r_{i}^{4}+28r_{i}^{6})-\frac{4+6r_{i}^{2}}{r_{i}^{4}}A_{2,i}\right)\!,\end{eqnarray}$$

(5.15)$$\begin{eqnarray}F_{i}^{(2)}(\unicode[STIX]{x1D746}_{i})=-\frac{1+\unicode[STIX]{x1D707}_{i}}{2}\text{e}^{-2V^{(1)}\unicode[STIX]{x1D70C}_{i}(1+\unicode[STIX]{x1D707}_{i})}(\unicode[STIX]{x1D6FF}V^{(1)}+8A_{1,i}),\end{eqnarray}$$

where $A_{1,i}$, $A_{2,i}$, $C_{0,i}$, $C_{1,i}$ and $C_{2,i}$ are unknown coefficients to be determined from the boundary conditions. Following the procedure described by Morozov & Michelin (Reference Morozov and Michelin2019c), we substitute $N_{i}^{(2)}(\boldsymbol{r}_{i})$ and $F_{i}^{(2)}(\unicode[STIX]{x1D746}_{i})$ into the boundary conditions at $\unicode[STIX]{x1D716}^{2}$ and obtain a set of algebraic equations in terms of the coefficients $A_{1,i}$, $A_{2,i}$, $C_{0,i}$, $C_{1,i}$ and $C_{2,i}$. Solvability condition of this set of equations yields the leading-order self-propulsion velocity of the droplet exposed to a weak concentration field,

(5.16)$$\begin{eqnarray}V^{(1)}=\frac{\unicode[STIX]{x1D6FF}\pm \sqrt{\unicode[STIX]{x1D6FF}^{2}+256G}}{32},\end{eqnarray}$$

where

(5.17a,b)$$\begin{eqnarray}G=G_{approach}\equiv -\frac{\text{e}^{-8DV^{(1)}}\left(1+8DV^{(1)}\right)}{4D^{2}}\quad \text{or}\quad G=G_{departure}\equiv \frac{1}{4D^{2}}\end{eqnarray}$$

for approaching and departing drops, respectively. We emphasise that $V^{(1)}\geqslant 0$ by definition. As a result, the departing configuration for which $G>0$ is used in (5.16) always admits a single solution, while the approaching configuration for which $G<0$ may have zero or two solutions depending on the magnitude of the concentration gradient.

Reconstructing the absolute relative velocity of the droplets is achieved as $V=-\unicode[STIX]{x1D716}V_{approach}^{(1)}$ (respectively $V=\unicode[STIX]{x1D716}V_{departure}^{(1)}$) for the approaching (respectively departing) case. In our numerical analysis, droplet approach corresponds to negative velocity (as shown in figures 4 and 6, 7 and 9), and the bifurcation diagram (5.16) is plotted in figure 11 to match this convention. Equation (5.16) corresponds to an imperfect transcritical bifurcation implying that not all of the branches of the bifurcation diagram (5.16) are stable, as shown in figure 11. Also note that the bifurcation diagram corresponding to (5.16) is quasistatic, i.e. it depends on time through the separation distance between the droplets, $D$, only.

Figure 11.

Bifurcation diagram (5.16): leading-order self-propulsion velocities of a pair of identical active droplets located at a distance $2d=2D/\unicode[STIX]{x1D716}$ with $D=10$. Top and bottom branches correspond to departing and approaching drops, respectively. A solid (respectively dashed) line denotes a stable (respectively unstable) propulsion regime.

In essence, equation (5.16) establishes that the regime of steadily approaching droplets does not exist when $G_{approach}<G_{c}\equiv -\unicode[STIX]{x1D6FF}^{2}/256<0$. As the droplets approach each other, $G$ is negative and increases in magnitude as $D$ decreases, up to a point where the quasi-steady approach branch ceases to exist, which is identified to the droplet’s velocity reversal and rebound dynamics. This event is associated with a critical centre-to-centre distance $\tilde{D}/\unicode[STIX]{x1D716}=\tilde{d}_{c}$ which provides an estimate for the rebound distance $d_{rb}$ and satisfies

(5.18)$$\begin{eqnarray}\left(\tilde{D}\unicode[STIX]{x1D6FF}/4\right)^{2}=4\left(1+\tilde{D}\unicode[STIX]{x1D6FF}/4\right)\text{e}^{-\tilde{D}\unicode[STIX]{x1D6FF}/4}\end{eqnarray}$$

and

(5.19)$$\begin{eqnarray}\tilde{d}_{c}\approx \frac{5.98}{Pe-4}.\end{eqnarray}$$

Estimation of the droplet rebound distance (5.19) reproduces the scaling obtained from our numerical simulations (§ 4.3). The prefactor of this asymptotic estimate $\tilde{d}_{c}$ of the rebound distance differs however slightly from the value estimated form our numerical results (figure 8). Such discrepancy is to be expected, since in numerical simulation rearrangement of the concentration field in the course of the rebound is not immediate, and in fact even in the quasi-steady framework presented here, $|V|=V_{0}/2\neq 0$ when the approaching branch ceases to exist (figure 11). Instead, droplets take some time to slow down. This regime of transitory approach is beyond the reach of our asymptotic analysis and, thus, the real rebound distance is shorter, compared to the theoretical prediction, equation (5.19).

The analysis presented here is focused on the case of two symmetric droplets. Yet, as emphasised throughout the analysis, the coupling between the droplets is purely chemical. Indeed, the flow field contribution is limited to the near-field dynamics of the flow field and does not influence the far-field signature (only the displacement of the droplet does). This further validates that, near the self-propulsion threshold, hydrodynamic interactions play a subdominant role. While chemical interactions are mediated through chemical gradients which decay as $1/d^{2}$, direct hydrodynamic interactions (i.e. the drift of a droplet in the flow field of the second one) would be dominated by the stresslet flow created by each droplet, which also decays as $1/d^{2}$. However, the intensity of the stresslet itself is weak for $Pe\approx Pe_{c}$, namely scaling as $\unicode[STIX]{x1D716}^{2}$ (Morozov & Michelin Reference Morozov and Michelin2019c), so that hydrodynamic interactions are $O(\unicode[STIX]{x1D716}^{4})$ in contrast to $O(\unicode[STIX]{x1D716}^{2})$ chemical interactions. As a consequence, the present asymptotic analysis applies exactly to the case of a droplet collision with a no-slip wall.

Finally, we note that the present approach and bifurcation diagram in (5.16) applies to any active droplet exposed to a concentration field, $C_{e}$, that allows for an expansion in powers of $\unicode[STIX]{x1D716}$, as shown in (5.12)–(5.13). The physical meaning of this mathematical requirement is twofold: (i) the evolution of $C_{e}$ must be slow, that is, in the laboratory frame $C_{e}$ should satisfy the steady diffusion equation up to $O(\unicode[STIX]{x1D716}^{3})$,

(5.20)$$\begin{eqnarray}\unicode[STIX]{x1D6FB}^{2}C_{e}=O(\unicode[STIX]{x1D716}^{3}),\end{eqnarray}$$

and (ii) the gradient of $C_{e}$ must be weak, namely, $C_{e}$ may only contribute to $N_{i}^{(2)}$, that is,

(5.21)$$\begin{eqnarray}\unicode[STIX]{x1D735}C_{e}=O(\unicode[STIX]{x1D716}^{2}).\end{eqnarray}$$

Any $C_{e}$ that satisfies these requirements, equations (5.20)–(5.21), may be seamlessly included in the superposition in (5.11) and subsequently expanded in powers of $\unicode[STIX]{x1D716}$ to obtain the corresponding value of $G$ for the bifurcation diagram (5.16).