Buoyancy-driven attraction of active droplets

Abstract

For dissolving active oil droplets in an ambient liquid, it is generally assumed that the Marangoni effect results in repulsive interactions, while the buoyancy effects caused by the density difference between the droplets, diffusing product and the ambient fluid are usually neglected. However, it has been observed in recent experiments that active droplets can form clusters due to buoyancy-driven convection (Krüger et al., Eur. Phys. J. E, vol. 39, 2016, pp. 1–9). In this study we numerically analyse the buoyancy effect, in addition to the propulsion caused by Marangoni flow (with its strength characterized by the Péclet number $Pe$). The buoyancy effects have their origin in (i) the density difference between the droplet and the ambient liquid, which is characterized by the Galileo number $Ga$; and (ii) the density difference between the diffusing product (i.e. filled micelles) and the ambient liquid, which can be quantified by a solutal Rayleigh number $Ra$. We analyse how the attracting and repulsing behaviour of neighbouring droplets depends on the control parameters $Pe$, $Ga$ and $Ra$. We find that while the Marangoni effect leads to the well-known repulsion between the interacting droplets, the buoyancy effect of the reaction product leads to buoyancy-driven attraction. At sufficiently large $Ra$, even collisions between the droplets can take place. Our study on the effect of $Ga$ further shows that with increasing $Ga$, the collision becomes delayed. Moreover, we derive that the attracting velocity of the droplets, which is characterized by a Reynolds number $Re_d$, is proportional to $Ra^{1/4}/( \ell /R)$, where $\ell /R$ is the distance between the neighbouring droplets normalized by the droplet radius. Finally, we numerically obtain the repulsive velocity of the droplets, characterized by a Reynolds number $Re_{rep}$, which is proportional to $PeRa^{-0.38}$. The balance of attractive and repulsive effect leads to $Pe\sim Ra^{0.63}$, which agrees well with the transition curve between the regimes with and without collision.

1. Introduction

The fundamental principles of microorganisms propulsion have gained attention across disciplines over the last few decades (Brennen & Winet 1977; Stone & Samuel 1996; Lauga & Powers 2009; Marchetti et al. 2013; Li & Ardekani 2016; Blackiston et al. 2021). Given the abundance of such microorganisms like bacteria and plankton within ecosystems (Hays, Richardson & Robinson 2005), studying both their individual and collective motion is critical for understanding ecosystem dynamics (Guasto, Rusconi & Stocker 2012). Interactions between microorganisms can stem from either physical factors, like hydrodynamics (Ramia, Tullock & Phan-Thien 1993; Ishikawa, Simmonds & Pedley 2006) or biological factors, such as visual (Trushin 2004) or chemical signals (Adler 1975). Disentangling these effects complicates the analysis of interactions in real microorganism colonies. To reduce complexity, researchers have turned to artificial microswimmers as simplified models, aiming to elucidate living microorganism interactions (Maass et al. 2016; Pedley 2016; Datt & Elfring 2019; Hokmabad et al. 2019; Gompper et al. 2020; Chen et al. 2021; Li 2022). Artificial microswimmers are designed to propel themselves by converting free energy from the environment into kinetic energy (Ogrin, Petrov & Winlove 2008), showing similar interactions observed among living microorganisms, such as chemotaxis, collective entrainment and cluster formation (Lozano et al. 2016; Maass et al. 2016; Jin, Krüger & Maass 2017; Lohse & Zhang 2020; Jin et al. 2021).

One extensively investigated type of artificial microswimmers is a dissolving active oil droplet suspended in water (Maass et al. 2016; Birrer, Cheon & Zarzar 2022; Dwivedi, Pillai & Mangal 2022). Propulsion of such active droplets is the Marangoni effect. The basic feature is that an imbalance in surfactant concentration at the droplet surface induces a tangential Marangoni flow, propelling the droplet forward (Herminghaus et al. 2014; Morozov & Michelin 2019a,b; Michelin 2023). This phenomenon extends to other coupled systems such as particles with catalytic surfaces, producing flow termed diffusiophoretic flow (Anderson 1989). With the Marangoni effect or diffusiophoresis being present, active droplets exhibit repulsion, because the concentration of surfactant molecules is lower between two close-by droplets or particles than that in the periphery, and the Marangoni flow or the diffusiophoretic flow will propel the active droplet or particles towards a higher surfactant concentration direction (where the surface tension is lower).

Repulsive interactions induced by Marangoni effects between active droplets have been well studied in numerous experimental and theoretical works. Moerman et al. (2017) performed a comprehensive experimental study, quantifying repulsive velocity and analysing the force-distance relationship between active droplets. Lippera, Benzaquen & Michelin (2020a) subsequently provided a theoretical analysis of the repulsive interaction, categorizing distinct motion modes. Further exploration by Lippera, Benzaquen & Michelin (2021) focused on obliquely colliding droplets, and identified whether the droplets interact directly or through their chemical wake. Besides the direct interactions, Jin et al. (2017) elucidated trail avoidance behaviour in active droplets, where the emitted filled micelles acted as chemical repellents, causing trajectory avoidance. This led to the identification of transit self-caging behaviour, reported by Daftari & Newhall (2022), where the active droplets are trapped due to trail avoidance.

While Marangoni effects from surface tension drive repulsion, the buoyancy effects arising from density inhomogeneity can introduce additional complexities via natural convection. The ‘jumping droplet’ phenomenon, examined by Li et al. (2019, 2021) and Li, Meijer & Lohse (2022), offers an example of intricate dynamics resulting from interplay of Marangoni and buoyancy forces. The droplet alternates between jumping due to Marangoni effects and sinking due to buoyancy. Another example is the evaporation of binary micro-droplets, where buoyancy and Marangoni forces compete, influencing convection and impacting the evaporation process (Edwards et al. 2018; Li et al. 2018, 2020; Diddens, Li & Lohse 2021). Finally, for a pair of fixed droplets, Lopez de la Cruz et al. (2022) reported an oscillatory flow near the droplet, which was attributed to the competition between Marangoni and buoyancy effects.

Coming back to dissolving an active droplet, Krueger et al. (2016) uncovered contrasting collective behaviours influenced by buoyancy in the presence of significant density differences. Active droplets were observed to attract and form hovering clusters, when solvent-droplet density difference exceeded a threshold. In a further study, Hokmabad et al. (2022) investigated the spontaneous rotation of the cluster formed by the attracted droplets, and attributed the attraction to the pusher-type squirmer. Very recently, Théry, Maaß & Lauga (2023) provided a theoretical hydrodynamic analysis of multi-squirmer interactions. Their work indicated that, although a pair of pusher swimmers exhibited attraction thanks to the far-field dynamics, the pusher clusters are hydrodynamically unstable. In any case, a thorough elucidation of the cluster mechanism is still missing, especially in terms of how the buoyancy effects drive attraction and override the Marangoni-induced repulsion.

Motivated by the above-mentioned studies, we focus on the collective behaviour of active droplets with buoyancy effects. The buoyancy effects stem from either the density difference between the droplet and ambient fluid or the solutal density difference between the dissolving product (i.e. filled micelles) and the ambient fluid. These effects are referred to as the ‘droplet buoyancy effect’ and the ‘product buoyancy effect’, respectively. In this study we quantitatively analyse the interplay between the Marangoni effect, the droplet buoyancy effect and the product buoyancy effect. We first simulate the interaction between a pair of active droplets, which allows us to elucidate the induced flow field and resultant droplet interactions. Subsequently, we present a model designed to predict attracting velocity using the method of reflections and Faxen's law. The model is validated with systems involving pairs or trios of interacting droplets. Then we analyse the repulsive velocity (based on the Marangoni flow) from simulations of a pair of fixed droplets. Finally, by comparing the repulsive and attractive velocity model, we obtain a good prediction for the regime transition of droplet collision.

The paper is organized as follows. Firstly, we describe the problem set-up in § 2. The numerical method and validations of the numerical scheme are provided in § 3. We first qualitatively analyse the role of the diffusiophoretic effect (characterized by $Pe$) and the product buoyancy effect in § 4 and then we analyse the role of the droplet buoyancy effect in § 5. Next, we develop a model to explain the attraction and calculate the attracting velocity and how it scales in §§ 6.1 and 6.2. The model is then further tested with the system of a pair of droplets and three droplets in § 6.3. Then the repulsive effect is analysed with cases of a pair of fixed droplets in § 7. Next, we discuss the influence of the product buoyancy effect to the horizontal self-propulsion of a single active droplet in § 8. Finally, concluding remarks are given in § 9.

2. Set-up and control parameters

We start with a pair of active droplets in the surfactant solution sketched in figure 1. The gradual solubilization of the oil into the surfactant micelles triggers a repulsive interaction due to the Marangoni effect (Jin et al. 2017). Simultaneously, the oil-filled micelles form in proximity to the droplet's surface. Additionally, we incorporate both the droplet buoyancy and the product buoyancy effects into our simulations.

Figure 1. The set-up of the system. (a) The two droplets are initially located in the middle of the domain. The droplet buoyancy effect, the product (indicated by yellow tails under the droplets) buoyancy effect and the diffusiophoretic effect are taken into consideration. The radius of the droplets is taken as the characteristic length. The domain size expressed in this length is then $16\times 16 \times 24$. The numerical grid resolution is $161 \times 161 \times 241$. The top and bottom boundary conditions are set as a solid wall (marked by grey plane) and the boundary conditions at $x$ and $y$ directions are periodic. (b) Because of the periodic boundary condition in the $y$ direction, the two droplets in the domain align with a series of droplets. In the $x$ direction the periodic boundary condition results in a balanced force. The distance between the two neighbouring droplets inside the domain is $\ell _1$ and the distance of neighbouring particles between inside and outside of the domain is $\ell _2$.

Considering similarities between diffusiophoresis and the Marangoni effect (Desai & Michelin 2021), for simplicity, we focus on the phoretic effect induced by the concentration gradient of the filled micelles and will use the corresponding terminology. The physical variables describing the system are the solutal concentration $\hat {c}$ and the velocity $\hat {\boldsymbol {u}}$. Note that all dimensional physical quantities are denoted with a hat (e.g. $\hat {c}$, $\hat {\boldsymbol {u}}$).

The droplets emit solute (filled micelles) at a rate $\alpha >0$. The concentration boundary condition at the droplet surface is given by

(2.1)\begin{equation} D\frac{\partial {\hat{c}}}{\partial {\hat{n}}}={-}\alpha, \end{equation}

where $D$ is the diffusion coefficient of the dissolution product, $\alpha$ the dissolution rate at the surface and $\partial \hat {c}/\partial \hat {n}$ the concentration gradient perpendicular to the surface. The same concentration boundary condition has widely been applied in previous studies of active droplets (Hokmabad et al. 2021; Michelin 2023).

The tangential concentration gradient at the surface induces a slip velocity, known as diffusiophoretic flow. The slip velocity magnitude $u_s$ is directly proportional to the local tangential concentration gradient, given by

(2.2)\begin{equation} \hat{u}_s=M\boldsymbol{\nabla}_s{\hat{c}}, \end{equation}

where $M$ is the mobility and $\boldsymbol {\nabla }_s$ represents the tangential gradient. In our analysis, we focus on the case where $M>0$, given that the filled micelles act as a chemical repellent.

We define $\hat {\rho }_0$ as the density of the ambient fluid without any dissolved product and $\hat {\rho }_d$ as the density of the droplet itself. Note that in the experiments by Krueger et al. (2016), the density difference among $\hat {\rho }_0$, $\hat {\rho }_d$ and the density of the dissolving product (filled micelles) $\hat {\rho }$ is lower than $3\,\%$. Accordingly, we adopt the Boussinesq approximation for the density difference, wherein the fluid's density $\hat {\rho }$ is assumed to be linearly proportional to the filled micelle concentration, expressed as

(2.3)\begin{equation} \hat{\rho}({\hat{c}})=\hat{\rho}_0(1+\beta \hat{c}), \end{equation}

where $\beta$ is the proportionality constant between the density and the product concentration.

The velocity field outside the droplets is governed by the Navier–Stokes equations and the product concentration field by the advection-diffusion equation. The equations are non-dimensionalized in the same way as in previous studies (Michelin, Lauga & Bartolo 2013; Chen et al. 2021), with $R$ for lengths, $c_0=\alpha R/D$ for concentrations and $\alpha M/D$ for velocities. Note that all dimensionless quantities are represented without a hat (e.g. $c$, $\boldsymbol {u}$). The non-dimensional governing equations are presented as

(2.4)$$\begin{gather} \frac{\partial c}{\partial t}+\boldsymbol{u}\boldsymbol{\cdot}\boldsymbol{\nabla} c= \frac{1}{Pe}\nabla^2 c, \end{gather}$$

(2.5a,b)$$\begin{gather}\frac{\partial \boldsymbol{u}}{\partial t}+(\boldsymbol{u}\boldsymbol{\cdot}\boldsymbol{\nabla})\boldsymbol{u}={-}\boldsymbol{\nabla} p+\frac{Sc}{Pe}\nabla^2 \boldsymbol{u}-\frac{RaSc}{Pe^2}c\boldsymbol{e}_z,\quad \boldsymbol{\nabla} \boldsymbol{\cdot} \boldsymbol{u}=0, \end{gather}$$

and the velocity of the droplet $\boldsymbol {U}_d$ satisfies

(2.6)\begin{equation} \frac{{\rm d}\boldsymbol{U}_d}{{\rm d}t} = \frac{3Sc}{4{\rm \pi} Pe}\int (\boldsymbol{\tau} \boldsymbol{\cdot} \boldsymbol{n})\,{\rm d}S-\frac{Ga^2Sc^2}{Pe^2}\boldsymbol{e}_z, \end{equation}

where $\int (\boldsymbol{\tau} \boldsymbol{\cdot} \boldsymbol{n}) \,{\rm d}S$ is the force integrated over the surface of the droplet and $\boldsymbol {e}_z$ is the unit vector of the $z$ axis.

The dimensionless control parameters of these equations are the Rayleigh number $Ra$, which represents the strength of the product buoyancy effect,

(2.7)\begin{equation} Ra=\frac{c_0 \beta R^3 g}{\nu D} , \end{equation}

the Péclet number $Pe$, which indicates the strength of the diffusiophoretic effect,

(2.8)\begin{equation} Pe=\frac{\alpha M R}{D^2}, \end{equation}

the Schmidt number $Sc$

(2.9)\begin{equation} Sc=\frac{\nu}{D}, \end{equation}

which is kept as a constant in our study, and the Galileo number $Ga$, which represents the strength of the droplet buoyancy effect,

(2.10)\begin{equation} Ga=\frac{\sqrt{|\hat{\rho}_d/\hat{\rho}_0-1|gR^3}}{\nu}. \end{equation}

The concentration boundary condition (2.1) at the droplet surface reads in non-dimensional form

(2.11)\begin{equation} \frac{\partial {c}}{\partial {n}}={-}1. \end{equation}

The non-dimensional version of the velocity boundary condition at the droplet surface (2.2) is

(2.12)\begin{equation} u_s=\boldsymbol{\nabla}_s c. \end{equation}

We apply periodic boundary conditions along the horizontal directions ($x$ and $y$) of the domain (figure 1b), and a solid wall boundary condition at the top and bottom of the domain. The concentration and velocity boundary conditions at the top and bottom are, respectively,

(2.13)\begin{equation} \frac{\partial c}{\partial z}=0 \end{equation}

and

(2.14)\begin{equation} \boldsymbol{u}=0. \end{equation}

3. Numerical methods and validation

The Navier–Stokes equations and advection–diffusion equation are solved using direct numerical simulation in Cartesian coordinates. The spatial discretization employs a central second-order finite difference scheme. The time integration is achieved using a fractional-step method. Nonlinear terms are explicitly computed using a low-storage third-order Runge–Kutta scheme and the viscous and diffusion terms are addressed with a Crank–Nicolson scheme (Verzicco & Orlandi 1996; Ostilla-Mónico et al. 2015; van der Poel et al. 2015; Spandan et al. 2018). The model for the droplet–droplet and droplet–wall collision is based on the spring-dashpot model (Costa et al. 2015).

Because the vicinity of the surface is adopted to satisfy the constant normal fluxes and slip velocity boundary condition, we cannot allow the gap between the droplets to reach zero in collisions. Therefore, droplets are considered to be in contact when the gap width is below 2 grid spacings.

The numerical set-up is shown in figure 1. The radius $R$ of the droplet is the characteristic length of the system and the domain size is $L_x\times L_y\times L_z =16 \times 16 \times 24$. Two droplets of unit radius are initially aligned along the $y$ axis at the centre of the domain with an initial distance $L_0=4$. Uniform grids $N_x\times N_y\times N_z=161\times 161\times 241$ are employed. Since there is the periodic boundary condition along the $y$ axis, two droplets in simulations are actually a part of a series of droplets aligning along the $y$ axis. In the $x$ direction the periodic boundary conditions ensure a balanced force. We use $\ell _1(t)$ and $\ell _2(t)$ to denote the distance between the droplet and its two neighbouring droplets along the $y$ axis, and we further define $\ell (t)=\text {min}(\ell _1(t), \ell _2(t))$ as the droplet's distance to its nearest neighbour, as shown in figure 1(b).

We explore a range of parameters based on the experimental data of Krueger et al. (2016): the Péclet number $0.5 \leq Pe \le 10$, Rayleigh number over the range $0.1\leq Ra\leq 245$ and the Galileo number $0 \leq Ga \le 0.19$. The Schmidt number in the experiments is of the order of $10^4$. However, the simulations at such high $Sc$ are challenging due to the very small diffusivity compared with viscosity. Therefore, in our study we set the Schmidt number to $Sc=100$.

Our code has been used to simulate diffusiophoretic particles. For the corresponding code validation, we refer the reader to our previous work (Chen et al. 2021). As a further validation, we test our code by simulating particle-laden flow incorporating both droplet buoyancy and product buoyancy effects, and compare with the existing results from the literature (Majlesara et al. 2020). The reference study considers the sedimentation of cold/hot (fixed temperature) spherical particles within a long vertical fluid channel and evaluate their terminal velocity. In that work, the product buoyancy effect is induced by the temperature variation, characterized by the Grashof number $Gr=Ra/Pr$. The particle buoyancy effect is quantified by the Reynolds number $Re$, which is linearly related to the Galileo number, $Re=\frac {2\sqrt {3}}{3}Ga$. Our code is applied to simulate the same cases as those presented by Majlesara et al. (2020), facilitating a direct comparison of the normalized terminal velocity $U_t/U_0$ across various $Re$ and $Gr$ values, where $U_0$ is the characteristic buoyancy velocity. The numerical results are plotted in figure 2; they agree very well with those by Majlesara et al. (2020).

Figure 2. Code validation for a settling particle at fixed temperature in a long vertical channel. The terminal velocity $U_t$, normalized by the reference velocity $U_0$, versus the Reynolds number $Re$, which is linearly correlated to the Galileo number, $Re=\frac {2\sqrt {3}}{3}Ga$. We show results for three Grashof numbers $Gr=Ra/Pr$. The results obtained by Majlesara et al. (2020) are indicated by filled symbols with the dashed lines. Our simulations are represented by the opened symbols, showing excellent agreement.

4. Effect of Péclet number and of Rayleigh number

In this section we first investigate the role of the Péclet number and of the Rayleigh number by simulating the interaction between a pair of droplets with $Sc=100$, $0.5 \leq Pe \leq 10$ and $0.1 \leq Ra \leq 245$. The droplet buoyancy effect is absent in this subsection ($Ga=0$) and will be analysed in § 5.

To demonstrate different interaction modes, we first focus on the cases $Pe=5$ and $Ra= 0.1$, $2$ and $245$ in figure 3, all for $Sc=100$ as throughout in this paper. For $Ra=0.1$, the diffusiophoretic effect is dominant. The mutual repulsion drives the droplets to the horizontally balanced positions ($1/4L_y$ and $3/4L_y$). The repulsion of the neighbouring droplets acts as a restoring force to the balanced position, while the droplets also experience a damping force due to the viscous drag. Therefore, the droplets perform a damped oscillation in the horizontal direction around the balanced positions. In the vertical direction the droplets rebound from the walls because of the concentration accumulation in between.

Figure 3. Concentration contours (a–c) and distance between droplets (d–f) as a function of time for a pair of droplets in the domain with parameters $Ga=0$, $Sc=100$, $Pe=5$ and various $Ra=$ $0.1$ (a,d), $Ra=2$ (b,e), $Ra=245$ (c,f). At the right we plot the distances $\ell _1$ and $\ell _2$ defined in figure 1 as a function of time. The droplet distances corresponding to the concentration contours are indicated as red filled circles in the plots.

As $Ra$ increases to $2$, the droplets approach their neighouring droplets. The reason is that the equi-distance balanced position becomes an unstable equilibrium due to the attraction between the droplets. In the end, the two droplets reach a new balanced point of finite distance $\ell < L_y/2$. Along the vertical direction ($z$), due to the stronger product buoyancy effect, the droplets first sediment to the bottom. Then the droplets gradually float up as the concentration between the wall and the droplets accumulate. As $Ra$ further increases, the terminal distance between the neighbouring droplets decreases. For $Ra=245$, the product buoyancy effect becomes even stronger. The droplets are more attractive to each other and collide in the end.

From the results, we find that different strengths of the product buoyancy effect lead to different terminal distances between the droplets along the horizontal axis. Therefore, we define the terminal distances $\ell _\infty$ to quantify that strength:

(4.1)\begin{equation} \ell_\infty =\underset{t \to \infty}{\text{lim}} \ell(t) = \underset{t \to \infty}{\text{lim}} (\text{min}(\ell_1(t), \ell_2(t))). \end{equation}

The dependence of $\ell _\infty$ on the Rayleigh number $Ra$ is shown in figure 4(a). We identify two different types of interaction according to $\ell _\infty$. First, $Ra < 50$: without collision, where the droplets remain at an equilibrium distance without colliding with each other. The distance $\ell _\infty$ between the droplets reduces as $Ra$ increases. Second, $Ra \geq 50$, with collision, where the droplets collide due to the sufficiently strong attraction driven by the product buoyancy effect.

Figure 4. (a) Terminal distance $\ell _\infty$ between the nearest droplets with $Ga=0$, $Sc=100$, $Pe=5$ and different $Ra$ from $0.1$ to $245$. Two interaction modes are identified, marked with different colours: $Ra \leq 50$, the droplets remain at an equilibrium distance (without collision: blue), $Ra \geq 50$, the droplets collide with each other due to the strong attraction (with collision: red). (b) The interaction modes for $Ga=0$, $Sc=100$, $0.5 \leq Pe \leq 10$, $0.1 \leq Ra \leq 245$. The blue circles represent the cases without collision, while the red triangles those with collision. The results indicate that a higher Pe results in a higher Ra threshold, above which the collision occurs.

We also simulate cases for different $Pe$ and $Ra$. The results can be classified into the two mentioned interaction modes, which are presented by different symbols in figure 4(b). From the plot, we find that a higher $Ra$ leads to collision while a higher $Pe$ leads to repulsion. The result in figure 4(b) indicates that there is competition between repulsion by diffusiophoresis and attraction by the product buoyancy effect. The two effects can be balanced, leading to the equilibrium distance for $0.5 \leq Ra \leq 50$ in figure 4(a). This complies with the experimental results of Krueger et al. (2016), who found that as the surfactant concentration ($Pe$) was increased, higher density differences ($Ra$) were needed for the collective behaviour to occur.

In summary, we investigate the influence of the diffusiophoretic effect and the product buoyancy effect to the droplet interaction. Specifically, the diffusiophoretic effect (characterized by $Pe$) results in repulsion between droplets, and the product buoyancy effect (characterized by $Ra$) leads to their attraction. The mechanism of the interaction will further be quantitatively investigated in §§ 6 and 7.

5. Effect of increasing Galileo number

In this section we analyse the influence of the droplet buoyancy effect for product-buoyancy-effect-dominant cases (high $Ra$), as quantified by the Galileo number $Ga$. We numerically investigate the interactions between a pair of droplets with $Pe = 5$, $Ra = 245$ and $0 \leq Ga \leq 0.19$.

Snapshots at different times are plotted in figure 5. For all examined cases with different $Ga$, the droplets collide in the end. To further analyse the interaction, we examine the temporal change of the horizontal distance $\ell$ and the vertical height of the droplets.

Figure 5. Concentration contours for a pair of droplets with $Sc=100$, $Pe=5$, $Ra=245$ and two different $Ga$: (a) $Ga=0.11$ (b) $Ga=0.19$.

We first plot the horizontal distance $\ell$ versus time $t$ in figure 6(a). The plot indicates that as $Ga$ increases, the waiting time for the collision to occur is longer. Next, we have a close inspection of the movement of droplets near the moment of collision through plotting $\ell -\ell _c$ as a function of $t_c - t$ in log-log scale in figure 6(b), where $\ell _c$ is the collision distance and $t_c$ is the collision time. Remarkably, all curves collapse on each other near the collision point. It suggests that the attracting behaviour of the droplets are mainly determined by $Ra$ and $Pe$, while the change of $Ga$ only leads to a delayed collision and does not influence the horizontal interaction.

Figure 6. The plot of distance $\ell$, $\ell -\ell _c$ and height $h$ versus time $t$ or $t_c-t$ with $Sc=100$, $Pe=5$, $Ra=245$ and different $Ga$. Here $t_c$ and $\ell _c$ are the collision time and distance. (a) The distance between the two droplets $\ell$ as a function of time for different $Ga$. Subfigure (b) shows $\ell -\ell _c$ and (c) $h$ along time $t_c-t$, where $\ell _c$ is the distance at the collision point and $t_c$ is the collision time.

We also plot the height $h$ (right $y$ axis) along $t_c-t$ in figure 6(c). As $Ga$ increases, the droplets wait for longer time before the occurrence of the approaching stage, and the rising velocity is also smaller for larger $Ga$. This is because a longer time is needed to build up a sufficiently large vertical concentration gradient to lift up a heavier droplet.

In summary, based on the numerical results, we find that $Ga$ will not change the horizontal interaction between droplets, but large $Ga$ delays the collision time by reducing vertical rising velocity.

6. Attraction model with buoyancy

In this section we further investigate the origin of attraction and develop a model to estimate the attracting velocity in the buoyancy-dominant cases. Employing the point heat source model, we first derive a scaling law for the horizontal velocity around a droplet, and then calculate the attractive velocity, using the methods of reflections and Faxen's law. Since the droplet buoyancy effect only leads to delayed collision, we neglect it, i.e. we assume $Ga=0$ throughout this section.

6.1. The velocity field near a single droplet

We start by simulating a single droplet to investigate the flow around it. Figure 7(a) shows the concentration and horizontal velocity ($v_y$) around a single droplet at $Ra=245$. From the fields, we observe a strong downwards plume, which leads to a higher concentration underneath the droplet. In the meantime, a horizontal flow is induced sidewards of the droplet. This inward flow drives the attraction between two nearby droplets.

Figure 7. (a) Concentration (left half) and velocity (right half) fields near a single droplet at $Ra=245$. The streamlines are shown by the white curves. The red dashed line is at the same height as the droplet. (b) The symmetric model is plotted in cylindrical coordinates $(r, z)$ to describes the flow near the droplet with buoyancy. The buoyancy induces a strong downwards flow under the droplet and a horizontal flow near the droplet. The width of the downwards flow is $h_1$ and the horizontal one $h_2$. In the simulation we define the width $h_1, h_2$ of each flow branch by the width between $10\,\%$ of the maximum vertical and horizontal velocity.

We represent the buoyancy-driven flow near a single droplet by the schematics in figure 7(b). Since the flow around a single droplet is axisymmetric, it is illustrated in cylindrical coordinates $(r, z)$. There is horizontal inward flow sidewards of the droplet and vertical downwards flow under the droplet. In the following analysis, we will see that this horizontal inward flow is the origin of the attractive motion between nearby droplets.

A similar case that has been well studied is the natural convection near a heat source or dissolutions source (Fujii 1963; Moses, Zocchi & Libchaberii 1993; Dietrich et al. 2016). Fujii (1963) theoretically studied the buoyancy-driven convection near a fixed heat source in the fluid, and quantitatively obtained the buoyancy-driven velocity. The theoretical results were later verified in experiments with a heating sphere in a fluid by Moses et al. (1993). Both the velocity and the width of the plume scale with the Rayleigh number $Ra$ (Fujii 1963; Moses et al. 1993; Dietrich et al. 2016):

(6.1)$$\begin{gather} v_z\sim \frac{D}{R}Ra^{1/2}, \end{gather}$$

(6.2)$$\begin{gather}h/R \sim Ra^{{-}1/4}. \end{gather}$$

We define the width $h_1, h_2$ of each flow branch as the distance between $10\,\%$ of the maximum vertical and horizontal velocity. Due to the limited domain size, the vertical and horizontal flow cannot attain the asymptotic velocity of 0, preventing us from using a smaller threshold for the $h_1, h_2$ definition. The normalized values $h_1/h_1(Ra=245)$ and $h_2/h_2(Ra=245)$ vs $Ra$ for different domain sizes ($16\times 16 \times 24$ and $30\times 30 \times 24$) obtained in simulations are plotted in figure 8. When $Ra$ is large enough ($Ra \geq 100$), the width of the channel well agrees with (6.2). For $Ra < 100$, a larger domain leads to a better agreement with theory. However, there is still deviation, which attributes to the finite-size effect and the existence of a strong enough diffusiophoretic effect.

Figure 8. The width of the downwards flow ($h_1$) and the horizontal flow ($h_2$), normalized by the corresponding height at $Ra=245$ for different $Ra$ at different domain sizes. We simulate the cases at different domain size $L_x \times L_y \times L_z =16\times 16\times 24$ denoted by $L_y=16$ and $30\times 30\times 24$ denoted by $L_y=30$. The blue and red symbols are correspondingly the numerical results for $h_1$ and $h_2$. The solid curve represents (6.2).

Given the width of both the horizontal and vertical flow, by continuity, we can further derive the relationship between the strength of the two velocities ($v_y$ for horizontal and $v_z$ for vertical), namely

(6.3) \begin{equation} v_y(r) \times 2{\rm \pi} rh_2\sim v_z \times {\rm \pi}h_1 ^2/4, \end{equation}

where $r$ refers to the horizontal distance from the droplet centre (along the red dashed curve in figure 7a). Then we define the local Reynolds number $Re_y (r)$ using the horizontal velocity $v_y (r)$. With (6.1) and (6.3), we obtain

(6.4)\begin{equation} Re_y (r)=\frac{v_y(r) R}{\nu}\sim \frac{1}{Sc} \frac{Ra^{1/4}}{r/R} . \end{equation}

We verify (6.4) with the numerical simulations of a single droplet in the domain. Note that due to the periodic boundary condition, the horizontal velocity is also influenced by the neighbouring droplets outside the domain,

(6.5)\begin{equation} \frac{Re_y (r)}{Ra^{1/4}}\sim \frac{R}{r}-\frac{R}{L_y-r}. \end{equation}

The results for different $Ra$ and $Pe$ are shown in figure 9. The numerical results agree well with the theory equation (6.5) for $r/R>4$. The results deviate near the droplet surface $r/R<4$, because the horizontal velocity reduces to zero approaching the droplet surface. A similar flow field was also observed and interpreted as a monopole component in the flow field around the active droplet by de Blois et al. (2019).

Figure 9. Plot of $Re_y(r)$ normalized by $Ra^{1/4}$ along the red dashed line in figure 7(a) for various distances $r$ to the droplet centre normalized by the radius $R$ of the droplet. The markers are the numerical results. The solid black curve represents relationship (6.5) with a fitted prefactor $0.021 \pm 0.003$ (Everitt & Skrondal 2010).

6.2. Droplet velocity using method of reflections and Faxen's law

In this section we apply Faxen's law and the method of reflections to account for the interactions between multiple droplets. The principle of the method of reflections is to perform successive approximations for the interaction of droplets within the fluid (Guazzelli & Morris 2011; Varma, Montenegro-Johnson & Michelin 2018). The velocity of the droplet is calculated iteratively, and in each step, the velocity of the droplet is updated with the disturbance from other droplets using Faxen's law (Guazzelli & Morris 2011). Despite the far-field assumption of the method, even for close distance $\ell /R\sim O(1)$, it reaches a surprisingly accurate result (Ishikawa et al. 2006; Spagnolie & Lauga 2012).

First, we consider a pair of active droplets (droplet 1 and 2) far apart. Since there is no external force and the droplet is isotropic, the droplets are stationary:

(6.6)\begin{equation} U_1^0=U_2^0=0. \end{equation}

Here $U_i^j$ represents the velocity of droplet $i$ after the $j$th reflection process.

Then in first reflection, we suppose that the droplets are only moderately far apart, and each droplet makes a disturbance at the velocity of the other. From (6.4), droplet 1 causes a fluid velocity disturbance at droplet 2:

(6.7)\begin{equation} u_2^0\sim \frac{D Ra^{1/4}}{\ell}. \end{equation}

Here $u_i^j$ is the fluid velocity disturbance at the centre of droplet $i$ caused by the other droplet after the $j$th reflection. According to Faxen's law, the velocity of droplet 2 due to the velocity disturbance caused by droplet 1 is (Guazzelli & Morris 2011, p. 87)

(6.8)\begin{equation} U_2^1= \left(1+\frac{R^2}{6}\nabla ^2 \right)u_2^0. \end{equation}

Since the two droplets are equivalent, the same velocity is obtained for droplet 1 after the first reflection.

Then we start with the second reflection, the velocity of the droplet obtained in the first reflection will cause disturbance to the other one. The fluid velocity caused by droplet 1 at the centre of droplet 2 is (Lamb 1932, p. 599)

(6.9)\begin{equation} u_2^1= \left(\frac{3R}{2\ell}-\frac{R^3}{2\ell^3} \right)U_1^1. \end{equation}

Again with Faxen's law, the velocity disturbance of droplet 2 after the second reflection is given by

(6.10)\begin{equation} U_2^2= \left(1+\frac{R^2}{6}\nabla ^2 \right)u_2^1. \end{equation}

For higher-order reflection, it is found that the velocity disturbance after reflection

(6.11)\begin{equation} U_2^n\sim O \left(\left(\frac{R}{\ell} \right)^{n-1} \right). \end{equation}

Therefore, we neglect the higher-order small terms, and the velocity of the droplet is approximated as

(6.12)\begin{equation} U(\ell)=U_2=U_2^0+U_2^1+U_2^2+O \left(\frac{R}{\ell} \right)=u_2^0+O \left(\frac{R}{\ell} \right)\sim \frac{DRa^{1/4}}{\ell}. \end{equation}

We define the Reynolds number of the droplet $Re_d$ by the droplet velocity $U$:

(6.13)\begin{equation} Re_d(\ell)=\frac{UR}{\nu}\sim\frac{D}{\nu}\frac{Ra^{1/4}}{\ell/R}=\frac{1}{Sc}\frac{Ra^{1/4}}{\ell/R}. \end{equation}

Here $Re_d (\ell )$ is different from $Re_y (r)$, where $Re_d(\ell )$ expresses the velocity of a droplet influenced by the other droplet at distance $\ell$, while $Re_y(r)$ corresponds to the fluid velocity at distance $r$ away from a single droplet.

Equation (6.13) indicates that the inward horizontal flow results in the attractive motion of a nearby droplet, but it considers the influence from only one neighbouring droplet. Note that the lateral boundaries are periodic. We consider the influence from the two neighbouring droplets and obtain

(6.14)\begin{equation} \frac{Re_d}{Ra^{1/4}}\sim \frac{R}{\ell}-\frac{R}{L_y-\ell}. \end{equation}

The equation can also be expressed as

(6.15)\begin{equation} \frac{Re_d}{Ra^{1/4}}\frac{L_y}{R}\sim \frac{1}{\ell/L_y}-\frac{1}{1-\ell/L_y}, \end{equation}

where $({Re_d}/{Ra^{1/4}})({L_y}/{R})$ is written as a function of $\ell /L_y$, facilitating the generalization of the equation to cases of different domain sizes.

6.3. Model validation

To validate the model, we calculate $({Re_d}/{Ra^{1/4}})({L_y}/{R})$ as a function $\ell /L_y$ for the droplet interaction discussed in § 4 and plot them in figure 10. The numerical results agree well with the theoretical model of (6.15). For low $Ra$ ($Ra\leq 50$), there is a deviation between the numerical results and (6.15) near the droplet, which can be explained by the influence of diffusiophoretic flow near the droplet.

Figure 10. Plot of $Re_d(\ell )L_y$ normalized by $Ra^{1/4}R$ versus the normalized distance $\ell /L_y$ between the pair of droplets of different $Pe$, $Ra$ and domain sizes. Specifically, the case where the domain size $20 \times 20 \times 24$ is denoted by $L_y=20$, the case with a domain size of $30 \times 30 \times 24$ is indicated by $L_y=30$ and the domain size of $16 \times 16 \times 24$ is indicated without explicit reference to $L_y$. The markers are for numerical results and the black solid line for relationship (6.15) with the fitted prefactor $0.012 \pm 0.002$ (Everitt & Skrondal 2010).

To test the domain size effect of the theory, we simulate cases of domain size $20 \times 20 \times 24$ and $30 \times 30 \times 24$, and the results are also plotted in figure 10, which collapse well with theory, confirming the applicability of the theory across different domain sizes.

To further test our theory, we simulate the case of three droplets initially located at the centre of the domain with $Ra=245$, $Sc=100$ and $Pe=5$, where the snapshots are given in figure 11(a). The horizontal velocity of the middle droplet remains at zero due to the symmetry about the middle axis. With our model of § 6.2, $Re_d$ follows

(6.16)\begin{equation} \frac{Re_d}{Ra^{1/4}}\sim \frac{R}{\ell}-\frac{R}{L_y-2\ell}. \end{equation}

Indeed, in figure 11(b), for large $\ell /R$, again an excellent agreement is seen between the numerical results and the prediction of (6.16). The fitted coefficient $0.013 \pm 0.001$ is also similar to that of the case for a pair of droplets in figure 10 $0.012 \pm 0.002$, which indicates the theory is consistent for both cases. However, it is worth noting that they are different from the fitted coefficient in figure 9, because the latter case is for a different set-up: the horizontal velocity around a single droplet.

Figure 11. (a) Snapshots at different times of concentration fields emerging from three neighbouring droplets. Here $Sc=100$, $Pe=5$ and $Ra=245$. (b) Plot of $Re_d$ normalized by $Ra^{1/4}$ versus the normalized smallest distance $\ell /R$. The symbols show the numerical results and the solid line shows relationship (6.16) with a fitted prefactor $0.013 \pm 0.001$.

7. Repulsive effect by diffusiophoresis

Given the good agreement between the attraction model with numerical results, now we further study the repulsive velocity from simulations. In the last section we found that the product buoyancy effect results in a solute convection near it, which leads to attraction between droplets. However, this solute convection also changes the flow field near the droplet, which changes the repulsive interaction due to diffusiophoretic effects. However, the concentration field with the product buoyancy effect is challenging to obtain theoretically. Therefore, in this section we analyse the repulsive effect of diffusiophoresis with numerical simulations.

We simulate a pair of droplets fixed at the centre of the domain (figure 1a) with a horizontal distance $\ell =3$ in the $y$ direction. Since the repulsive diffusiophoretic motion mainly comes from the slip velocity induced by the concentration field, we estimate the repulsive velocity by the integral of slip velocity at the surface (Stone & Samuel 1996):

(7.1)\begin{equation} U_{rep}=\frac{1}{4{\rm \pi} R^2}\int_S{{\boldsymbol{u}}_s\,{\rm d}S}. \end{equation}

We define $Re_{rep}=U_{rep}R/\nu$ to represent the repulsive interaction. We simulate the cases of different $Pe$ and $Ra$, and the resulting concentration field is shown in figure 12. The relationship between $Re_{rep}$ and $Pe$, $Ra$ is shown in figure 13. Figure 13(a) shows the relationship between the normalized $Re_{rep}$ and $Pe$. The results indicate that the diffusiophoretic effect leads to repulsive motion, which agrees with our conclusions in § 4, and $Re_{rep}$ is proportional to $Pe$. Figure 13(b) shows $Re_{rep}$ for different $Ra$, and we find that a stronger buoyancy effect reduces the repulsive velocity between droplets. We fit the results with a power law ansatz and get

(7.2)\begin{equation} Re_{rep}\sim PeRa^{{-}0.38}. \end{equation}

Figure 12. The concentration field for a pair of fixed droplets at distance $3$ for $Pe=5$, $Sc=100$ and three different $Ra$ numbers: (a) $Ra=1$, (b) $Ra=10$, (c) $Ra=200$. Here $\theta$ is the angle between the bottom point and the maximum concentration point to represents the plume position.

Figure 13. (a) Normalized droplet repulsive Reynolds number $Re_{rep}/Re_{rep}(Pe=1)$ for different $Pe$. The plot shows that the $Re_{rep}/Re_{rep}(Pe=1)$ is proportional to $Pe$. The fitting result is $Re_{rep} \sim Pe^{0.98\pm 0.03}$, which indicates a linear relationship with fitting exponent standard error $0.03$. (b) Plot of $Re_{rep}/Pe$ vs $Ra$. The solid line represents the fitted function, which shows that $Re_{rep}/Pe$ is proportional to $Ra^{-0.38 \pm 0.02}$ (Everitt & Skrondal 2010).

To better understand the decrease of $Re_{rep}$ with increasing $Ra$, we first study the influence of $Ra$ on its surrounding concentration field. We plot the concentration gradient near a single droplet of different $Ra$ at $Pe=5$ from simulations in figure 14(a). It indicates that the concentration gradient has a significant drop as $Ra$ increases. This can be rationalized as follows: as $Ra$ increases, buoyancy-driven convection reduces the thickness of the concentration boundary layer (Fujii 1963; Dietrich et al. 2016). As the surface concentration gradient remains constant (2.11), a reduction in the boundary layer thickness leads to a lower local concentration gradient near the droplet.

Figure 14. (a) Concentration gradient $|\partial c/ \partial r|$ at normalized distance $r/R$ for different $Ra$ near a single droplet, which indicates that the concentration gradient decreases as $Ra$ increases. The inset shows the normalized concentration profile. (b) The plume position $\theta$ vs $Ra$, which reflects the fact that the plume is pulled more towards the droplet bottom as $Ra$ increases.

Moreover, we find that buoyancy also influences the position of the plume at the droplet surface. Through the concentration field in figure 12, as $Ra$ increases, the plume moves closer to the bottom of the droplet. To evaluate the effect, we define $\theta$ as the angle between the maximum concentration point and the droplet bottom point to represent the position of the plume as indicated in figure 12. Figure 14(b) shows the change of the plume position for different $Ra$ at $Pe=5$. This finding thus suggests that a stronger buoyancy effect (higher $Ra$) pulls the plume towards the bottom point and this can reduce the horizontal component of the repulsive diffusiophoretic velocity.

We acknowledge that the arguments above are handwaving and qualitative. The complex system dynamics resulting from the coupling between convection and the concentration field makes a theoretical derivation of the relationship between $Re_{rep}$ and $Ra$ too challenging.

However, if we combine the equations for the attractive (6.13) and repulsive velocities (7.2), we obtain

(7.3)\begin{equation} Pe\sim Ra^{0.63}, \end{equation}

which perfectly describes the transition between the attracting and the repelling regimes; see figure 4(b). This plot nicely reflects that the mechanism behind the interaction between droplets is the competition between the attractive force by buoyancy and the repulsive force by diffusiophoresis.

8. Discussion

In the previous studies without the buoyancy effect (Michelin et al. 2013), the Péclet number is found to dominate the droplet motion: if the Péclet number is low ($Pe<4$), the droplet remains stable and stationary, and at high Pe ($Pe>4$), the droplet breaks the symmetry and undergoes self-propulsion. A thorough study on onset of propulsion for the case without buoyancy has been conducted in our previous work Chen et al. (2021). However, in the current work, we find that the same dynamics applies to droplet behaviour across a range of Péclet numbers ($0.5< Pe<10$) and the onset of self-propulsion ($Pe=4$) does not change the behaviour of the droplet.

To elucidate the phenomena and explain the influence of buoyancy effect to the onset of self-propulsion, we conduct simulations of a single droplet positioned in the centre of the domain. To minimize interaction of the droplet with the boundaries of the computational domain, we employ a domain size of $L_x \times L_y \times L_z= 30\times 30 \times 24$. Specifically, we simulate the cases for $Pe=0.5$ and $10$, with varying Rayleigh numbers ($Ra=0$, $0.001$, $0.1$, $1$ and $10$), and initiate the droplets with an initial horizontal velocity perturbation $\Delta u_y=0.001$.

The results are presented in figure 15. We observe that at low Péclet number ($Pe=0.5$), if $Ra \leq 0.001$, the droplet remains stable (figure 15a), consistent with the findings of Michelin et al. (2013). However, if $Ra \geq 0.1$, the droplet initially settles down to the bottom and then rises up from the wall (figure 15b). This suggests that at low $Pe$, a sufficiently strong product buoyancy effect induces downward motion.

Figure 15. Concentration contours for a single droplet settlement with parameters $Ga=0$, $Sc=100$ and (a) $Pe=0.5$, $Ra=0.001$; (b) $Pe=0.5$, $Ra=1$; (c) $Pe=10$, $Ra=0$; (d) $Pe=10$, $Ra=10$. The inset of figure (d) shows velocity vectors around the droplet (the white dashed square region) and the vertical velocity field ($v_z$).

At high Péclet number ($Pe=10$), when $Ra=0$, the droplet spontaneously breaks symmetry and propels itself in the same direction as the perturbation ($\kern 1.5pt y$ direction), as shown in figure 15(c). However, as $Ra$ increases, as evident from the droplet trajectory in figure 16, horizontal motion is suppressed, and the droplet tends to move vertically (in the $z$ direction) rather than horizontally (in the $y$ direction). If $Ra \geq 0.1$, the droplet settles down without much horizontal motion. Two factors contribute to the phenomena: firstly, as shown in the inset of figure 15(d), the product buoyancy effect leads to a downward flow near the droplet, which continuously drives the droplet in a downward motion, and secondly, as depicted in figure 14(b), higher Ra numbers draw the plume closer to the bottom of the droplet, limiting the horizontal diffusiophoretic motion.

Figure 16. Numerical results of the droplet trajectory for cases with parameters $Ga=0$, $Pe=10$ and varying Rayleigh numbers: $Ra=0$, $0.001$, $0.1$, $1$, $10$, $100$. The results suggest that at high $Ra$, horizontal motion is suppressed and the droplet motion is in the vertical direction.

In summary, the involvement of the product buoyancy effect significantly alters the motion and interaction of droplets. We find that the droplet self-propulsion in the horizontal direction is restrained at Rayleigh numbers $Ra\geq 0.1$ for both high and low Péclet numbers. Consequently, there is no discernible distinction in the interaction mechanism between low and high Péclet number scenarios within this Pe and Ra range.

Moreover, while previous studies have extensively examined droplet motion near a solid wall with the influence of the diffusiophoretic effect (Lippera et al. 2020b; Desai & Michelin 2022), our findings indicate that the involvement of the product buoyancy effect significantly alters the droplet motion. A comprehensive investigation of the buoyancy effect to droplet motion requires further research in the future.

9. Summary and conclusions

We have studied the interaction between droplets influenced by diffusiophoretic, droplet buoyancy and product buoyancy effects, characterized by the Péclet number ($Pe$), Galileo number ($Ga$) and Rayleigh number ($Ra$) correspondingly. We have simulated the cases over a range of $Pe$, $Ra$ and $Ga$, while keeping $Sc$ constant at $100$.

Through numerical simulations of a pair of droplets, we have found that the product buoyancy effect leads to the attractive motion between droplets, while the Marangoni/diffusiophoretic effect results in repulsion. A larger $Pe$ corresponds to a higher $Ra$ threshold for droplet collision. When $Ra$ is relatively small, the distance between droplets reaches an equi-distance equilibrium, and with increasing $Ra$, the closest balanced distance decreases, indicating that the product buoyancy weakens the Marangoni/diffusiophoretic repulsion. For sufficiently high Rayleigh numbers ($Ra\geq 50$), the droplets collide with each other. Additionally, we investigate the droplet buoyancy effect and reveal that the attracting behaviour is similar for different $Ga$, and the change of $Ga$ only leads to a delayed collision.

The simulation of a single droplet indicates that the attraction stems from convective flow induced by the density difference between the dissolving product and ambient fluid. This insight leads to a model describing the horizontal velocity near the droplet, elucidating that the local Reynolds number is inversely proportional to the distance from the droplet, as shown in (6.4).

Building on this model, we derive the equation for the attracting velocity of the droplet at high $Ra$ using Faxen's law and the method of reflections. The attracting velocity scales proportionally with $Ra^{1/4}$ and inversely with the distance between droplets. The relationship is validated through simulations involving two or three droplets.

Then we investigate the repulsive effect by simulating a pair of fixed droplets and the repulsive velocity is approximated by the integral of the slip velocity (7.1). We find that $Re_{rep}$, which represents the repulsive velocity, is proportional to $PeRa^{-0.38}$. The linear dependence of $Re_{rep}$ on $Pe$ results from a larger diffusiophoretic repulsive force for larger $Pe$, while the $Ra$ dependence indicates that higher $Ra$ leads to reduced repulsion through a reducing horizontal concentration gradient and altering plume positioning as shown in figure 14.

Combining the scaling relations of the attractive and repulsive velocity, we obtain $Pe\sim Ra^{0.63}$, aligning perfectly with the transition curve between the attractive and repulsive regime in figure 4(b). This indicates that the interplay between attractive buoyancy and repulsive diffusiophoretic forces governs droplet interaction.

We also find that for a single droplet, the sufficiently strong product buoyancy effect ($Ra \geq 0.1$) restricts the diffusiophoretic horizontal motion because the product buoyancy effect drives a downward flow near the droplet, and also drives the concentration plume closer to the bottom of the droplet, limiting the horizontal motion due to diffusiophoretic effect. This restriction also results in no discernible distinction in the interaction mechanism between low and high Péclet number with $Ra \geq 0.1$, because in both cases, the horizontal self-propulsion is restricted.

This work contributes to the understanding of the interaction between active droplets, and specifically reveals the significant role of the dissolving product buoyancy. The insights obtained clarify the attraction observed in experiments by Krueger et al. (2016), and we also offer a model to predict the velocity of interacting active droplets. Our results provide a framework to understand the droplet attraction induced by the convective flow and suggest a mechanism to manipulate collective behaviours by tuning the buoyancy.

Many questions remain open. For example, how to determine the cluster size for multiple droplets? How does the flow field change if droplets are near a fluid–air interface? How to quantitatively determine the threshold Rayleigh number above which the droplets collide with each other and show collective behaviours? With the obtained insights into the attraction here, we hope it is seen as worthwhile to further explore the formation and motion of a cluster of active particles.