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.


Introduction
The fundamental principles of microorganisms propulsion have gained attention across disciplines over the past 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 such as bacteria and plankton in our ecosystem (Hays et al. 2005), studying their individual and collective motion is critical for understanding the dynamics of the entire ecosystem (Guasto et al. 2012).The interactions between microorganisms can be purely physical, i.e. based on hydrodynamics (Ramia et al. 1993;Ishikawa et al. 2006), or biological, i.e. based on visual signals (Trushin 2004) or by chemical signals (Adler 1975).Disentangling these effects makes it difficult to analyze the various interactions in real microorganism colonies.To reduce the complexity, in recent years artificial microswimmers as a simplified model have been investigated in order to understand the interactions between living microorganisms (Pedley 2016;Maass et al. 2016;Datt & Elfring 2019;Hokmabad et al. 2019;Gompper et al. 2020;Chen et al. 2021;Li 2022).Such artificial microswimmers are designed to propel themselves by converting free energy from the environment into kinetic energy (Ogrin et al. 2008).Similar interactions as those between living microorganisms are observed, such as chemotaxis, collective entrainment, and cluster formation (Maass et al. 2016;Lozano et al. 2016;Jin et al. 2017;Lohse & Zhang 2020;Jin et al. 2021).
One extensively studied type of artificial microswimmer is a dissolving active oil droplet floating in water (Maass et al. 2016).The driving mechanism behind the propulsion of such active droplets is the Marangoni effect.The basic feature is that whenever there is an inhomogeneity of surfactant concentration at the surface of the droplet, the consequent surface tension difference generates a tangential Marangoni flow adjacent to the surface, which leads to the self-propulsion of the droplet (Herminghaus et al. 2014;Morozov & Michelin 2019a,b;Michelin 2022).This effect can also be generalised to other coupled systems such as the particles with catalytic surfaces.The resulting flows are then referred to as diffusio-phoretic flow (Anderson 1989).With the Marangoni effect or diffusiophoresis be present, such active droplets become repulsive.This simply happens because the Marangoni flow or the diffusiophoretic flow will propel the active droplet or particles towards higher surfactant concentration direction (where the surface tension is lower), and the concentration of surfactant molecules is lower between two close-by droplets or particles than that in the periphery.
Repulsive interactions induced by Marangoni effects between active droplets have been well studied in numerous experimental and theoretical works.A clear experimental observation of the repulsion was conducted by Moerman et al. (2017), who quantitatively measured the repulsive velocity for a pair of active droplets and analyzed the relations between the repulsive force and their distance.Later on, Lippera et al. (2020) theoretically analyzed the repulsive interaction between a pair of droplets and identified different motion modes.In a further study, Lippera et al. (2021) investigated the repulsive interactions for a pair of obliquely-colliding droplets, and identified whether the droplets interact directly or through their chemical wake.Besides the direct interaction between active droplets, Jin et al. (2017) found the active droplets also show trail avoidance behavior.They reported that the active droplet emitted filled micelles in the wake play a role as chemical repellents and cause trajectory avoidance.Based on such observation, Daftari & Newhall (2022) developed a mathematical model to mimic the trail-avoidance with multiple active droplets.They observed that the active droplets could trap themselves due to the trail-avoidance, a feature which has been called transit self-caging behavior.
However, besides the Marangoni effect caused by the difference in surface tension, also buoyancy effect caused by the inhomogeneity in the density field can affect the flow, namely by natural convection.One interesting example of a droplet for which the interplay between Marangoni and buoyancy forces leads to rich dynamics is the phenomenon of the "jumping droplet" (Li et al. 2019(Li et al. , 2021(Li et al. , 2022)).In that case the droplet repeatedly jumps up due to the Marangoni effect, and then slowly sinks due to buoyancy (Li et al. 2019(Li et al. , 2022)).Another example is the evaporation of a binary micro-droplet, where buoyancy competes with Marangoni forces and can drive the convection inside the droplet, playing a crucial role in the evaporation process (Edwards et al. 2018;Li et al. 2018Li et al. , 2020;;Diddens et al. 2021).Finally for a pair of fixed droplets, Lopez de la Cruz et al. (2022) reported an oscillatory flow near the droplet, which again is triggered by the competition between Marangoni and buoyancy effects.
Coming back to dissolving active droplet, also here, despite the repulsive interaction by the Marangoni effect, Krueger et al. (2016) observed opposing collective behaviours when buoyancy is significant in the collective droplet system.They reported that the active droplets attract each other and form clusters hovering in the fluid and ascribed this finding to buoyancy effect because the collective behavior only occurs when the density difference between the solvent and the droplet is above a certain threshold.In a further study, Hokmabad et al. (2022) investigated the spontaneous rotation of the cluster formed by the attracted droplets.However, an explanation of the detailed mechanism of the attraction is still missing, especially on how the buoyancy effect drives the attraction and overtakes the repulsion driven by the Marangoni effect.
Inspired by the above-mentioned studies, we focus on the collective behavior of active droplets with buoyancy effects.Note that there are two different types of buoyancy effects, either by the density difference between the droplet and ambient fluid, or by the solutal density difference between the dissolving product (i.e.filled micelles) and the ambient fluid.Hereafter, we call them the "droplet buoyancy effect" and the "product buoyancy effect", respectively.In this work, we will quantitatively analyze 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.The numerical simulations allow us to capture the flow field around the droplet and how it induces the droplet interaction.Then we develop a model to predict the attracting velocity based on the method of reflections and on Faxen's law.We test our model with systems of two and three interacting droplets.Then we analyze the repulsive velocity (based on the Marangoni flow) from simulations of a pair of fixed droplets.Finally, by comparing these results with those from the attractive velocity model due to buoyancy effect, we obtain a good prediction for the regime transition of droplet collision.
The paper is organized as follows: We first describe the problem setup in Section 2. The numerical method and validations of the numerical scheme are provided in Section 3. We first qualitatively analyze the role of the diffusiophoretic effect (characterized by ) and the product buoyancy effect in Section 4 and then we analyze the role of the droplet buoyancy effect in Section 5. Next, we develop a model to explain the attraction and calculate the attracting velocity and how it scales in section 6.1 and 6.2 .The model is then further tested with the system of a pair of droplets and three droplets in section 6.3.Then the repulsive effect is analyzed with cases of a pair of fixed droplets in section 7. Finally, concluding remarks are given in section 8. F 1. The setup 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 characteristic length.The domain size expressed in this length is then 16 × 16 × 24.The numerical grid resolution is 161 × 161 × 241.The top and bottom boundary conditions are set as solid wall (marked by gray plane) and the boundary conditions at  and  directions are periodic.(b) Because of the periodic boundary condition in  direction, the two droplets in the domain align with a series of droplets.In  direction, the periodic boundary condition results in a balanced force.The distance between the two neighboring droplets inside the domain is ℓ 1 and the distance of neighboring particles between inside and outside of the domain is ℓ 2 .

Setup 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 causes a repulsive interaction via the Marangoni effect (Jin et al. 2017).Simultaneously, the oil-filled micelles are generated near the droplet surface.Besides that, the droplet buoyancy effect and product buoyancy effect are taken into consideration in the simulations.
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 to describe the system are the solutal concentration ĉ and the velocity û.Note that all dimensional physical fields are marked with a hat (e.g.ĉ, û), while the dimensionless ones are without a hat (e.g., u).
The droplets emit a solute (filled micelles) at a rate  > 0. The concentration boundary condition at the droplet surface is given by where  is the diffusion coefficient of the dissolution product,  the dissolution rate at the surface, and  ĉ/ n the concentration gradient normal to the surface.The tangential concentration gradient at the surface induces a slip velocity, which is the so-called diffusiophoretic flow.The magnitude of the slip velocity   is proportional to the local tangential concentration gradient, given by û = ∇  ĉ, (2.2)

Focus on Fluids articles must not exceed this page length
where  is the mobility and ∇  represents the tangential gradient.Since the filled micelles are emitted as a chemical repellent, the case of  > 0 is considered.We define ρ0 as the density of the surrounding fluid without any dissolved product and ρ as the density of the droplet itself.Note that in the experiments by Krueger et al. (2016), the density difference among ρ0 , ρ and the density of the dissolving product (filled micelles) ρ is lower than 3%.Therefore, we consider the density difference within the Boussinesq approximation, i.e. the density of the fluid ρ is assumed to be linearly proportional to the filled micelle concentration ρ( ĉ) = ρ0 (1 +  ĉ), (2.3) where  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 nondimensionalized by  for lengths,  0 = / for concentrations, and / for velocities.Then the non-dimensional governing equations can be written as and the velocity of the droplet   satisfies where ∫ ( • ) is the force integrated over the surface of the droplet and   is the unit vector of the  axis.
The dimensionless control parameters of these equations are the Rayleigh number , which represents the strength of the product buoyancy effect, the Péclet number , which indicates the strength of the diffusiophoretic effect, which is kept as a constant in our study, and the Galileo number , which represents the strength of the droplet buoyancy effect, (2.10) The concentration boundary condition (equation (2.1)) at the droplet surface reads in nondimensional form   = −1. (2.11) The non-dimensional version of the velocity boundary condition at the droplet surface (equation (2.2)) is (2.12) We apply periodic boundary conditions along the horizontal directions ( and ) of the domain (figure 1 (b)), and 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,   = 0, (2.13) and  = 0. (2.14)

Numerical methods and validation
The Navier-Stokes equations and advection-diffusion equation are solved using direct numerical simulation (DNS) in Cartesian coordinates.We spatially discretize the equations using the central second-order finite difference scheme.Uniform staggered grids are used in all directions.The time integration is accomplished by using a fractional-step method.The non-linear terms are computed explicitly by a low-storage third-order Runge-Kutta scheme and the viscous and diffusion terms by a Crank-Nicolson scheme (Verzicco & Orlandi 1996;van der Poel et al. 2015;Ostilla-Mónico et al. 2015;Spandan et al. 2018).The model for the droplet-droplet and droplet-wall collision is based on the spring-dashpot model by 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.Therefore, when the gap width is below 2 grid spacings, we assume that the droplets are in contact.
The numerical setup is shown in figure 1.The radius  of the droplet is the characteristic length of the system, and the domain size is   ×   ×   = 16 × 16 × 24.Two droplets of unit radius are initially aligned along the  axis at the center of the domain with an initial distance  0 = 4.We use uniform grids   ×   ×   = 161 × 161 × 241.Since there is the periodic boundary condition along the  axis, the simulations with two droplets are actually a part of a series of droplets aligning along the  axis.In  direction, the periodic boundary conditions result in a balanced force.We use ℓ 1 () and ℓ 2 () to denote the distance between the droplet and its two neighboring droplets along the  axis, and we further define ℓ() = min(ℓ 1 (), ℓ 2 ()) as the droplet's distance to its nearest neighbor as shown in figure 1  (b).
We will take a range of parameters based on the data in the experiment by Krueger et al. (2016): Péclet number 0.5  10, Rayleigh number over the range 0.1  245 and Galileo number 0  0.19.The Schmidt number in the experiments is at the order of 10 4 .However, the simulations at such high  are challenging due to the very small diffusivity compared to viscosity.Therefore, in our study, we set the Schmidt number  = 100.
Our code has been used to simulate diffusiophoretic particles.For the corresponding code validation, we refer the readers to our previous work (Chen et al. 2021).As a further validation, we test our code by simulating particle-laden flow with both droplet buoyancy effect and product buoyancy effect, and comparing with the existing results from the literature (Majlesara et al. 2020).These authors consider the cases about sedimenting cold/hot (fixed temperature) spherical particles in a long vertical fluid channel and study their terminal velocity.In that work, the product buoyancy effect is induced by the temperature variation, characterized by the Grashof number  = /.The particle buoyancy effect is characterized by the Reynolds number , which carries the same information as the Galileo number,  = 2 √ 3 3 .We apply our code to simulate the same cases as Majlesara et al.  (2020), and compare the normalized terminal velocity   / 0 for various  and , where  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).

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  = 100, 0.5  10, and 0.1  245.The droplet buoyancy effect is absent in this subsection ( = 0), and will be analyzed in subsection 5.
To demonstrate different interaction modes, we first focus on the cases  = 5 and  = 0.1, 2 and 245 in figure 3, all for  = 100 as throughout in this paper.For  = 0.1, the diffusiophoretic effect is dominant.The mutual repulsion drives the droplets to the horizontally balanced positions (1/4  and 3/4  ).The repulsion of the neighbouring droplets acts as 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 horizontal direction around the balanced positions.In vertical direction, the droplets rebound from the walls because of the concentration accumulation in between.
As  increases to 2, the droplets approach their neighoring 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 ℓ <   /2.Along the vertical direction (), 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  further increases, the terminal distance between the neighbouring droplets decreases.For  = 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

Without collision
With collision F 4. (a) Terminal distance ℓ ∞ between the nearest droplets with  = 0,  = 100,  = 5 and different  from 0.1 to 245.The dashed curve is a guide to the eyes.Two interaction modes are identified, marked with different colors:  50, the droplets remain at an equilibrium distance (without collision: blue),  50, the droplets collide with each other due to the strong attraction (with collision: red).(b) The interaction modes for  = 0,  = 100, 0.5  10, 0.1  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.define the terminal distances ℓ ∞ to quantify that strength: The dependence of ℓ ∞ on the Rayleigh number  is shown in figure 4 (a).We identify two different types of interaction according to ℓ ∞ : (a)  < 50: without collision, where the droplets remain at an equilibrium distance without colliding with each other.The distance ℓ ∞ between the droplets reduces as  increases.(b)  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  and .The results can be classified into the two mentioned interaction modes, which are presented by different symbols in figure 4 (b).The results indicate that there is competition between repulsion by diffusiophoresis and attraction by the product buoyancy effect.Higher  results in a higher  threshold, above which the collision occurs.This complies with the experimental results by Krueger et al. (2016), who find that the surfactant concentration () is increased, higher density differences () are needed for collective behavior to occur.
In summary, we numerically observe the interaction between droplets.We find very similar features as in the experiments by Krueger et al. (2016).While the diffusiophoretic effect (characterised by ) results in repulsion between droplets, the product buoyancy effect (characterised by ) leads to their attraction.We identify two different interaction modes: when the diffusiophoretic effect is dominant, the droplets reach a finite distance without collision; when the product buoyancy effect is dominant, the droplets collide in the end.In the next section, we will further investigate the role of the droplet buoyancy effect.

Effect of increasing Galileo number
In this section, we analyze the influence of the droplet buoyancy effect, as quantified by the Galileo number .We numerically investigate the interactions between a pair of droplets with  = 5,  = 245 and 0  0.19.
Snapshots at different times are plotted in figure 5.For all examined cases with different  , the droplets collide in the end.To further analyze the interaction, we will examine the temporal change of the horizontal distance ℓ and the vertical height of the droplets.
We first plot the horizontal distance ℓ versus time  in figure 6 (a).The plot indicates that as  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 ℓ − ℓ  as a function of   −  in log-log scale in figure 6 (b), where ℓ  is the collision distance and   is the collision time.Remarkably, all curves collapse on each other near the collision point.It suggests that the attracting behavior of the droplets are mainly determined by  and , while the change of  only leads to a delayed collision.
We also plot the height ℎ (right y axis) along   −  in figure 6 (c).As  increases, the droplets wait for longer time before the occurrence of the approaching stage, and the rising velocity is also smaller for larger .This is because a longer time is needed to build up a sufficiently large vertical concentration gradient to lift up a heavier droplet.

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  = 0 throughout this section.

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 (  ) around a single droplet at  = 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 (, ).There is horizontal inward flow sidewards of the droplet, and vertical downwards flow under the droplet.
A similar case that has been well studied is the natural convection near a heat source or dissolutions source (Fujii 1963;Moses et al. 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  (Fujii 1963;Moses et al. 1993;Dietrich et al. 2016): ℎ/ ∼  −1/4 .(6.2) We define the width ℎ 1 , ℎ 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 Pe=1 Ra=100 Pe=3 Ra=100 Pe=10 Ra=245 Pe=5 Ra=10 Pe=5 Ra=50 Pe=5 Ra=100 Pe=5 Ra=150 Pe=5 Ra=200 Pe=5 Ra=245 Theory F 9.   () normalized by  1/4 along the red dashed line in figure 7(a) for various distances  to the droplet center normalized by the radius  of the droplet.The markers are the numerical results and the solid lines are a guide to the eye.The solid black curve represents relationship (6.5) with a fitted prefactor 0.021.smaller threshold for the ℎ 1 , ℎ 2 definition.The normalized values ℎ 1 /ℎ 1 ( = 245) and ℎ 2 /ℎ 2 ( = 245) versus  obtained in simulations are plotted in figure 8.When  is large enough ( 100), the width of the channel well agrees with (6.2).For  < 100, the numerical results deviate, because of the existence of a strong enough diffusiophoretic effect.
Given the width of both the horizontal and the vertical flow, by continuity, we can further derive the relationship between the strength of the two velocities (  for horizontal and   for vertical), namely   () × 2 ℎ 2 ∼   × ℎ 2 1 /4, (6.3)where  refers to the horizontal distance from the droplet center (along the red dashed curve in figure 7 (a)).Then we define the local Reynolds number   () using the horizontal velocity   ().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 neighboring droplets outside the domain, (6.5) The results for different  and  are shown in figure 9.The numerical results agree well with the theory equation (6.5) for / > 4. The results deviate near the droplet surface / < 4, because the horizontal velocity reduces to zero approaching the droplet surface.

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 et al. 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 Theory F 10.   (ℓ) normalized by  1/4 versus the normalized distance ℓ/ between the pair of droplets.The markers are for numerical results with lines to guide the eye and the black solid line for relationship (6.14) with the fitted prefactor 0.012.using Faxen's law (Guazzelli & Morris 2011).Despite the far-field assumption of the method, even for close distance ℓ/ ∼  (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: where    represents the velocity of droplet  after the 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: where    is the fluid velocity disturbance at the center of droplet  caused by the other droplet after the th reflection.According to the Faxen's law, the velocity of droplet 2 due to the velocity disturbance caused by droplet 1 is (Guazzelli & Morris (2011), 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 center of droplet 2 is (Lamb (1932), p. 599): Again with Faxen's law, the velocity disturbance of the droplet 2 after the second reflection is given by: (6.10) 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   by the droplet velocity : , where   (ℓ) expresses the velocity of a droplet influenced by the other droplet at distance ℓ, while   () corresponds to the fluid velocity at distance  away from a single droplet.Equation (6.13) considers the influence from only one neighboring droplet.Note that the lateral boundaries are periodic.We consider the influence from the two neighboring droplets and obtain (6.14)

Model validation
The   of the droplets at different distances ℓ of different  and  obtained from simulations are shown in figure 10.The numerical results collapse for sufficiently large droplet separation ℓ.For large distances ℓ/, the numerical results can be described by equation (6.14) which excellently agrees with the data for the velocity of the droplet especially for high .For low  ( 50), there is a deviation between the numerical results and the relationship (6.14) near the droplet, which can be explained by the influence of diffusiophoretic flow near the droplet.
To further test our theory, we simulate the case of three droplets initially located at the center of the domain with  = 245,  = 100, and  = 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 subsection 6.2,   follows: Indeed, in figure 11 (b), for large ℓ/ again an excellent agreement is seen between the numerical results and the prediction of equation (6.15).

Repulsive effect by diffusiophoresis
Given the good agreement between the attraction model with numerical results, now we further study the repulsive velocity from simulations.
We simulate a pair of droplets fixed at the center of the domain (figure 1 (a)) with a horizontal distance ℓ = 3 at  direction.Since the repulsive diffusiophoretic motion mainly comes from the slip velocity induced by concentration field, we estimate the repulsive velocity by the integral of slip velocity at the surface (Stone & Samuel 1996): (7.1)We define  rep =  rep / to represent the repulsive interaction.We simulate the cases of different  and , and the resulting concentration field is shown in figure 12.The relationship between  rep and ,  is shown in figure 13. Figure 13 (a) shows the relationship between the normalized  rep and .The results indicate that the diffusiophoretic effect leads to repulsive motion, which agrees with our conclusions in Section 4, and  rep is proportional to . Figure 13 (b) shows  rep for different , 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  rep ∼  −0.38 .
(7.2) To better understand the decrease of  rep with increasing , we first study the influence of  on its surrounding concentration field.We plot the concentration gradient near a single droplet of different  at  = 5 from simulations in Figure 14 (a).It indicates that the concentration gradient has a significant drop as  increases.This can be rationalized as follows: As  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 (equation (2.11)), a reduction in the boundary layer thickness leads to a lower local 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  increases, the plume moves closer to the bottom of the droplet.To evaluate the effect, we define  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  at  = 5.This finding thus suggests that a stronger buoyancy effect (higher ) 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  rep and  too challenging.
However, if we combine the equations for the attractive (6.13) and repulsive velocities (7.2), we obtain  ∼  0.63 , (7.3) 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.

Summary & Conclusions
We have studied the interaction between droplets with diffusiophoretic effect, droplet buoyancy effect and product buoyancy effect.The corresponding parameters are Péclet number (), Galileo number () and Rayleigh number ().We have simulated the cases over a range of , , and , with  being fixed at 100.
For a pair of droplets, using numerical simulations, 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  results in a larger  threshold, above which droplet collision occurs.If the Rayleigh number is sufficiently small, the distance between droplets reaches an equi-distance equilibrium, and as  increases, the closest balanced distance decreases, which indicates that the product buoyancy weakens the repulsion caused by the Marangoni/diffusiophoretic effect.For sufficiently high Rayleigh numbers ( 50), the droplets collide with each other.Then we investigated the influence of droplet buoyancy effect and found that the attracting behavior is similar for different , and the change of  only leads to a delayed collision.
With the simulation of a single droplet, we have found that the attraction originates from convective flow induced by the density difference between the dissolving product and ambient fluid.Based on this, we have created a simple model which well describes the horizontal velocity near the droplet.The local Reynolds number is inversely proportional to the distance from the droplet as shown in equation (6.4).
With the above model as a starting point, we have obtained the equation for the attracting velocity of the droplet at high  with Faxen's law and the method of reflections.The attracting velocity is proportional to  1/4 and inversely proportional to the distance between the droplets.The results have been verified by the simulation results for cases with two and three droplets.
Then we have investigated the repulsive effect by simulating the case of a pair of fixed droplets and the repulsive velocity was approximated by the integral of the slip velocity (7.1).We have found that  rep , which represents the repulsive velocity, is proportional to  −0.38 .The linear dependence of  rep on  is simply due to a larger diffusiophoretic repulsive force for larger .In contrast, the -dependence of  rep is more complicated.It reflects that an increasing  leads to a smaller horizontal concentration gradient and favours the plume to be closer to the bottom point of droplet, which reduces the repulsive velocity.
Combining the scaling relations of the attractive and repulsive velocity, we obtain  ∼  0.63 , which perfectly describes the transition curve between the attractive and repulsive regime in figure 4 (b).This indicates that the mechanism behind the interaction between droplets are the competition between attractive buoyancy force and repulsive diffusiophoretic force.
The present work contributes to the understanding of the interaction between active droplets, and specifically reveals the significant role played by the dissolving product buoyancy.It shows that product buoyancy can lead to attractive motion between active particles, which helps us understand the attraction of active droplet in the experiments of Krueger et al. (2016).We have proposed a simple model to predict the velocity of the interacting active droplets.The results provide a framework to understand the droplet attraction induced by the convective flow.Moreover, the present work reveals a possible way to change the collective behaviors by tuning the buoyancy.
In our simulation, the propulsion of the active droplet is simply modelled as diffusiophoreis.Alternatively, we could have taken Marangoni flow.Until now, buoyancy-driven attractive motions are only observed in the cases of active droplets but scarcely in phoretic particles, possibly due to the difficulties to generate large enough density difference between the product and ambient fluid by phoretic particles.
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.

F 2 .
Code validation for a settling particle at fixed temperature in a long vertical channel.The terminal velocity   , normalized by the reference velocity  0 , versus the Reynolds number , which is linearly correlated to Galileo number,  = 2 √ 3 3 .We show results for three Grashof numbers  = /.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.
contours (left panels) and distance between droplets (right panels) as funtion of time for a pair of droplets in the domain with parameters  = 0,  = 100,  = 5 and various  = 0.1 (a),  = 2 (b) ,  = 245 (c).At the right we plot the distances ℓ 1 and ℓ 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.
plot of distance ℓ, ℓ − ℓ  and height ℎ versus time  or   −  with  = 100,  = 5,  = 245 and different .  and ℓ  are the collision time and distance.(a) The distance between the two droplets ℓ as a function of time for different .(b) ℓ − ℓ  and (c) ℎ along time   − , where ℓ  is the distance at collision point and   is the collision time.
a) Concentration (left half) and velocity (right half) fields near a single droplet at  = 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 coordinate (, ) 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 ℎ 1 and the horizontal one ℎ 2 .In the simulation, we define the width ℎ 1 , ℎ 2 of each flow branch by the width between 10% of the maximum vertical and horizontal velocity.
width of the downwards flow (ℎ 1 ) and the horizontal flow (ℎ 2 ) normalized by the corresponding height at  = 245 for different .The blue and red symbols are correspondingly the numerical results for ℎ 1 and ℎ 2 .The solid curve represents (6.2).
a) Snapshots at different times of concentration fields emerging from three neighboring droplets.Here  = 100,  = and  = 245.(b)   normalized by  1/4 versus the normalized smallest distance ℓ/.The symbols show the numerical results and the solid line shows relationship (6.15) with a fitted prefactor 0.013.
concentration field for a pair of fixed droplets at distance 3 for  = 5,  = 100 and three different  numbers: (a)  = 1, (b)  = 10, (c)  = 200. is the angle between the bottom point and the maximum concentration point to represents the plume position.
a) Normalized droplet repulsive Reynolds number  rep / rep ( = 1) for different .The plot shows that the rep / rep ( = 1) is proportional to .(b)  rep / versus .The solid line represents the fitted function, which shows that  rep / is proportional to  −0.38 .
a) Concentration gradient |/ | at normlized distance / for different  near a single droplet, which indicates the concentration gradient decreases as  increases.The inset shows the normalized concentration profile.(b) The plume position  versus , which reflects that the plume is more pulled towards the droplet bottom as  increases.