## 1. Introduction

The fundamental principles of microorganisms propulsion have gained attention across disciplines over the last few decades (Brennen & Winet Reference Brennen and Winet1977; Stone & Samuel Reference Stone and Samuel1996; Lauga & Powers Reference Lauga and Powers2009; Marchetti *et al.* Reference Marchetti, Joanny, Ramaswamy, Liverpool, Prost, Rao and Simha2013; Li & Ardekani Reference Li and Ardekani2016; Blackiston *et al.* Reference Blackiston, Lederer, Kriegman, Garnier, Bongard and Levin2021). Given the abundance of such microorganisms like bacteria and plankton within ecosystems (Hays, Richardson & Robinson Reference Hays, Richardson and Robinson2005), studying both their individual and collective motion is critical for understanding ecosystem dynamics (Guasto, Rusconi & Stocker Reference Guasto, Rusconi and Stocker2012). Interactions between microorganisms can stem from either physical factors, like hydrodynamics (Ramia, Tullock & Phan-Thien Reference Ramia, Tullock and Phan-Thien1993; Ishikawa, Simmonds & Pedley Reference Ishikawa, Simmonds and Pedley2006) or biological factors, such as visual (Trushin Reference Trushin2004) or chemical signals (Adler Reference Adler1975). 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.* Reference Maass, Krüger, Herminghaus and Bahr2016; Pedley Reference Pedley2016; Datt & Elfring Reference Datt and Elfring2019; Hokmabad *et al.* Reference Hokmabad, Baldwin, Krüger, Bahr and Maass2019; Gompper *et al.* Reference Gompper2020; Chen *et al.* Reference Chen, Chong, Liu, Verzicco and Lohse2021; Li Reference Li2022). Artificial microswimmers are designed to propel themselves by converting free energy from the environment into kinetic energy (Ogrin, Petrov & Winlove Reference Ogrin, Petrov and Winlove2008), showing similar interactions observed among living microorganisms, such as chemotaxis, collective entrainment and cluster formation (Lozano *et al.* Reference Lozano, Ten Hagen, Löwen and Bechinger2016; Maass *et al.* Reference Maass, Krüger, Herminghaus and Bahr2016; Jin, Krüger & Maass Reference Jin, Krüger and Maass2017; Lohse & Zhang Reference Lohse and Zhang2020; Jin *et al.* Reference Jin, Chen, Maass and Mathijssen2021).

One extensively investigated type of artificial microswimmers is a dissolving active oil droplet suspended in water (Maass *et al.* Reference Maass, Krüger, Herminghaus and Bahr2016; Birrer, Cheon & Zarzar Reference Birrer, Cheon and Zarzar2022; Dwivedi, Pillai & Mangal Reference Dwivedi, Pillai and Mangal2022). 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.* Reference Herminghaus, Maass, Krüger, Thutupalli, Goehring and Bahr2014; Morozov & Michelin Reference Morozov and Michelin2019*a*,Reference Morozov and Michelin*b*; Michelin Reference Michelin2023). This phenomenon extends to other coupled systems such as particles with catalytic surfaces, producing flow termed diffusiophoretic flow (Anderson Reference Anderson1989). 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.* (Reference Moerman, Moyses, Van Der Wee, Grier, Van Blaaderen, Kegel, Groenewold and Brujic2017) performed a comprehensive experimental study, quantifying repulsive velocity and analysing the force-distance relationship between active droplets. Lippera, Benzaquen & Michelin (Reference Lippera, Benzaquen and Michelin2020*a*) subsequently provided a theoretical analysis of the repulsive interaction, categorizing distinct motion modes. Further exploration by Lippera, Benzaquen & Michelin (Reference Lippera, Benzaquen and Michelin2021) focused on obliquely colliding droplets, and identified whether the droplets interact directly or through their chemical wake. Besides the direct interactions, Jin *et al.* (Reference Jin, Krüger and Maass2017) 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 (Reference Daftari and Newhall2022), 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.* (Reference Li, Diddens, Prosperetti, Chong, Zhang and Lohse2019, Reference Li, Diddens, Prosperetti and Lohse2021) and Li, Meijer & Lohse (Reference Li, Meijer and Lohse2022), 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.* Reference Edwards, Atkinson, Cheung, Liang, Fairhurst and Ouali2018; Li *et al.* Reference Li, Lv, Diddens, Tan, Wijshoff, Versluis and Lohse2018, Reference Li, Diddens, Segers, Wijshoff, Versluis and Lohse2020; Diddens, Li & Lohse Reference Diddens, Li and Lohse2021). Finally, for a pair of fixed droplets, Lopez de la Cruz *et al.* (Reference Lopez de la Cruz, Diddens, Zhang and Lohse2022) 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.* (Reference Krueger, Bahr, Herminghaus and Maass2016) 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.* (Reference Hokmabad, Nishide, Ramesh, Krüger and Maass2022) 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 (Reference Théry, Maaß and Lauga2023) 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.* Reference Jin, Krüger and Maass2017). 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.

Considering similarities between diffusiophoresis and the Marangoni effect (Desai & Michelin Reference Desai and Michelin2021), 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

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.* Reference Hokmabad, Dey, Jalaal, Mohanty, Almukambetova, Baldwin, Lohse and Maass2021; Michelin Reference Michelin2023).

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

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.* (Reference Krueger, Bahr, Herminghaus and Maass2016), 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

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 Reference Michelin, Lauga and Bartolo2013; Chen *et al.* Reference Chen, Chong, Liu, Verzicco and Lohse2021), 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

*a*,

*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

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,

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

the Schmidt number $Sc$

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

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

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

We apply periodic boundary conditions along the horizontal directions ($x$ and $y$) of the domain (figure 1*b*), 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,

and

## 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 Reference Verzicco and Orlandi1996; Ostilla-Mónico *et al.* Reference Ostilla-Mónico, Yang, Van Der Poel, Lohse and Verzicco2015; van der Poel *et al.* Reference van der Poel, Ostilla-Mónico, Donners and Verzicco2015; Spandan *et al.* Reference Spandan, Lohse, de Tullio and Verzicco2018). The model for the droplet–droplet and droplet–wall collision is based on the spring-dashpot model (Costa *et al.* Reference Costa, Boersma, Westerweel and Breugem2015).

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.* (Reference Krueger, Bahr, Herminghaus and Maass2016): 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.* Reference Chen, Chong, Liu, Verzicco and Lohse2021). 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.* Reference Majlesara, Abouali, Kamali, Ardekani and Brandt2020). 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.* (Reference Majlesara, Abouali, Kamali, Ardekani and Brandt2020), 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.* (Reference Majlesara, Abouali, Kamali, Ardekani and Brandt2020).

## 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.

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:

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.

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.* (Reference Krueger, Bahr, Herminghaus and Maass2016), 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.

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.

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.

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 Reference Fujii1963; Moses, Zocchi & Libchaberii Reference Moses, Zocchi and Libchaberii1993; Dietrich *et al.* Reference Dietrich, Wildeman, Visser, Hofhuis, Kooij, Zandvliet and Lohse2016). Fujii (Reference Fujii1963) 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.* (Reference Moses, Zocchi and Libchaberii1993). Both the velocity and the width of the plume scale with the Rayleigh number $Ra$ (Fujii Reference Fujii1963; Moses *et al.* Reference Moses, Zocchi and Libchaberii1993; Dietrich *et al.* Reference Dietrich, Wildeman, Visser, Hofhuis, Kooij, Zandvliet and Lohse2016):

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.

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

where $r$ refers to the horizontal distance from the droplet centre (along the red dashed curve in figure 7*a*). 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

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,

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.* (Reference de Blois, Reyssat, Michelin and Dauchot2019).

### 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 Reference Guazzelli and Morris2011; Varma, Montenegro-Johnson & Michelin Reference Varma, Montenegro-Johnson and Michelin2018). 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 Reference Guazzelli and Morris2011). 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.* Reference Ishikawa, Simmonds and Pedley2006; Spagnolie & Lauga Reference Spagnolie and Lauga2012).

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:

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:

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 Reference Guazzelli and Morris2011, p. 87)

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 Reference Lamb1932, p. 599)

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

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

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

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

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

The equation can also be expressed as

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.

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

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.

## 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 1*a*) 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 Reference Stone and Samuel1996):

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

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 Reference Fujii1963; Dietrich *et al.* Reference Dietrich, Wildeman, Visser, Hofhuis, Kooij, Zandvliet and Lohse2016). 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.

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

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.* Reference Michelin, Lauga and Bartolo2013), 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.* (Reference Chen, Chong, Liu, Verzicco and Lohse2021). 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 15*a*), consistent with the findings of Michelin *et al.* (Reference Michelin, Lauga and Bartolo2013). However, if $Ra \geq 0.1$, the droplet initially settles down to the bottom and then rises up from the wall (figure 15*b*). This suggests that at low $Pe$, a sufficiently strong product buoyancy effect induces downward motion.

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.

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.* Reference Lippera, Morozov, Benzaquen and Michelin2020*b*; Desai & Michelin Reference Desai and Michelin2022), 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.* (Reference Krueger, Bahr, Herminghaus and Maass2016), 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.

## Funding

We greatly appreciate valuable discussions with M. Assen, U. Jain and C. Maass. We acknowledge the support from the Netherlands Center for Multiscale Catalytic Energy Conversion (MCEC), an NWO Gravitation program funded by the Ministry of Education and support from the ERC-Advanced Grant ‘DDD’ under the project number 740479. The simulations in this work were carried out on the national e-infrastructure of SURFsara, a subsidiary of SURF cooperation, the collaborative ICT organization for Dutch education and research, MareNostrum 4 that is based in Spain at the Barcelona Computing Center (BSC) under PRACE projects 2018194742, 2020225335 and 2020235589, on Irene at Trés Grand Centre de calcul du CEA (TGCC) under PRACE project 2019215098 and on Marconi successor at CINECA, Italy under PRACE project 2019204979. K.L.C. is supported by the Natural Science Foundation of China under grant no. 92052201.

## Declaration of interests

The authors report no conflict of interest.