1. Introduction
Stable density stratification in liquids can readily be generated from temperature or compositional variations due to the presence of gravity on Earth. Particles, microorganisms, drops and bubbles suspended in stratified fluids play a significant role in a variety of natural phenomena and industrial applications, such as the transport of pollutants in the atmosphere (Pereira et al. Reference Pereira, Freitas, Moraes, Ferreira, Shimabukuro, Rao and Longo2009), algal blooms in ponds and lakes (Anderson et al. Reference Anderson, Cembella and Hallegraeff2012; Paerl et al. Reference Paerl, Gardner, McCarthy, Peierls and Wilhelm2014), and encapsulation and coating in food science and cosmetics engineering (Dubey et al. Reference Dubey, Shami and Rao2009). Suspended objects and a background stratification are strongly coupled: the settling or rising of objects is altered by the stratification, whereas the former can influence the stratification itself by entraining nearby fluid in the form of a column or a tail. Accurate prediction of the motion of suspended objects in stratified environments is crucial for modelling and predicting the transport and distribution of swimming organisms, and for optimising industrial processes involving mixing and dispersion (Durham et al. Reference Durham, Kessler and Stocker2009; Mathijssen et al. Reference Mathijssen, Lisicki, Prakash and Mossige2023; Shaik & Elfring Reference Shaik and Elfring2024).
Much effort has been devoted to understanding the dynamics of suspended objects in stratified fluids, ranging from rigid spheres to deformable bubbles or droplets (Hanazaki et al. Reference Hanazaki, Konishi and Okamura2009b
; Doostmohammadi et al. Reference Doostmohammadi, Dabiri and Ardekani2014; Magnaudet & Mercier Reference Magnaudet and Mercier2020; More & Ardekani Reference More and Ardekani2020, Reference More and Ardekani2023). Previous work on rigid spheres demonstrates that the drag on a sinking (rising) sphere is enhanced by the stratification, because the moving sphere entrains low-density (high-density) liquid that reduces the net buoyancy force (Srdić-Mitrović et al. Reference Srdić-Mitrović, Mohamed and Fernando1999; Torres et al. Reference Torres, Hanazaki, Ochoa, Castillo and Van Woert2000). Recent reviews survey the physics underlying drag enhancement in stratified fluids with a linear density gradient, as well as two miscible fluids separated by a more or less sharp interface (Magnaudet & Mercier Reference Magnaudet and Mercier2020; More & Ardekani Reference More and Ardekani2023). Theoretical and experimental studies have shown how drag enhancement is controlled by key parameters such as the Reynolds number
$Re$
, the Froude number
$Fr$
, and the Prandtl number
$\textit{Pr}$
(Hanazaki et al. Reference Hanazaki, Kashimoto and Okamura2009a
; Zhang et al. Reference Zhang, Mercier and Magnaudet2019; More & Ardekani Reference More and Ardekani2020).
Recently, there has been much interest in the intricate dynamics of moving droplets or bubbles in stratifications due to their deformability and tangential surface motion (Blanchette & Shapiro Reference Blanchette and Shapiro2012; Bayareh et al. Reference Bayareh, Doostmohammadi, Dabiri and Ardekani2013; Mandel et al. Reference Mandel, Zhou, Waldrop, Theillard, Kleckner and Khatri2020; Shaik & Ardekani Reference Shaik and Ardekani2020; Zeng et al. Reference Zeng, Chong, Wang, Diddens, Li, Detert, Zandvliet and Lohse2021). In particular, the presence of a concentration or temperature gradient around a droplet generates differences in interfacial tension, leading to Marangoni flow (Chang et al. Reference Chang, Wilcox and Lefever1979; Schwabe et al. Reference Schwabe, Preisser and Scharmann1982; Balasubramaniam & Chai Reference Balasubramaniam and Chai1987; Rednikov et al. Reference Rednikov, Ryazantsev and Velarde1994; Haj-Hariri et al. Reference Haj-Hariri, Shi and Borhan1997; Zhang et al. Reference Zhang, Subramanian and Balasubramaniam2001). Such flow can drive self-propulsion of droplets, enabling diverse potential applications such as targeted cargo delivery, decontamination, and chemical analysis and sensing (Maass et al. Reference Maass, Krüger, Herminghaus and Bahr2016; Yi et al. Reference Yi, Sanchez, Gao and Yu2016; Jurado-Sánchez & Wang Reference Jurado-Sánchez and Wang2018). Early studies investigated droplet motions by focusing separately on either the surface tension gradient (Young et al. Reference Young, Goldstein and Block1959; Leven & Newman Reference Leven and Newman1976; Chang et al. Reference Chang, Wilcox and Lefever1979; Chun & Wuest Reference Chun and Wuest1979; Schwabe et al. Reference Schwabe, Preisser and Scharmann1982; Kamotani et al. Reference Kamotani, Ostrach and Vargas1984; Chen & Stebe Reference Chen and Stebe1996) or the density gradient in the surrounding medium (Bayareh et al. Reference Bayareh, Doostmohammadi, Dabiri and Ardekani2013; Shaik & Ardekani Reference Shaik and Ardekani2020). However, more intricate dynamics arises from the coupling of these two effects (Blanchette & Shapiro Reference Blanchette and Shapiro2012; Mandel et al. Reference Mandel, Zhou, Waldrop, Theillard, Kleckner and Khatri2020). For example, in a linear density stratification with a length scale much larger than the droplet size, the Marangoni force can be comparable to the buoyancy of (say) a rising lighter droplet, which decreases the time required for the droplet to rise to its maximum height (Mandel et al. Reference Mandel, Zhou, Waldrop, Theillard, Kleckner and Khatri2020). Drops falling through a sharp density and surface tension stratification with a length scale smaller than the drop size can bounce back up due to Marangoni flow (Blanchette & Shapiro Reference Blanchette and Shapiro2012). Recent studies show that the coupling between Marangoni flow and gravity-driven flow can make the drop levitate or bounce continuously (Li et al. Reference Li, Diddens, Prosperetti, Chong, Zhang and Lohse2019, Reference Li, Diddens, Prosperetti and Lohse2021). The transition from levitation to bouncing is due to an oscillatory Marangoni instability. While most studies emphasise vertical motion along the concentration gradient, earlier theoretical and experimental investigations of oceanic mixing reveal that stably stratified diffusive fluids with inclined boundaries can also generate spontaneous lateral flows (Phillips Reference Phillips1970; Wunsch Reference Wunsch1970). In such systems, fluid streams along the sloping boundaries within a boundary layer whose Reynolds number scales inversely with the Prandtl number. To the best of our knowledge, there has been little investigation of the lateral motion of droplets due to the interplay of solutal Marangoni flow, natural convection and stratification.
Here we report a novel behaviour of silicone oil droplets in a continuously stratified ethanol–water mixture, namely lateral self-propulsion with or without vertical oscillation. To explore how physical properties such as droplet size and concentration gradient influence the dynamics, we performed experiments with droplet radii ranging from
$\sim 400$
${\unicode{x03BC}}$
m to 800
${\unicode{x03BC}}$
m in stratified ethanol–water mixtures with an ethanol concentration gradient
$\mathrm{d} w_e/\mathrm{d} z$
in the range
$\sim 0.16{-}0.28$
mm
$^{-1}$
, where
$w_e$
is the mass fraction of ethanol in the mixture. We quantify the lateral self-propulsion modes of the droplets via their mean square displacement (MSD) as a function of time. Theoretical scaling analysis is used to derive the critical condition for the onset of lateral motion and the steady lateral speed far beyond onset, which are in good agreement with experimental observations. The phenomenon described here is an interesting example of Marangoni-driven self-propulsion, akin to active matter systems in which autonomous droplet motions are fuelled by energy released from vertical interfacial Marangoni effects.
2. Methods
2.1. Preparation of stratified ethanol–water mixtures
Before preparing the stratified ethanol–water mixtures, ethanol and deionised water were degassed in a vacuum chamber at 0.1 bar for 30 min to prevent bubble generation during mixing. The temperature of the laboratory and the liquids used was maintained at
$20\,^{\circ}$
C during experiments. As illustrated in figure 1(a), 20 ml of water were injected into the bottom of a cuvette (50
$\times$
50
$\times$
50 mm
$^3$
) containing 30 ml ethanol to create a density stratification in the transition region between the ethanol and water layers. The injection rate was less than 2 ml min−1 to minimise flow disturbance to the stratification. To estimate the ethanol concentration gradient, we used methylene blue dissolved in water (40 mg l−1) as an indicator. The diffusion coefficient of methylene blue,
$(0.67\pm0.13)\times 10^{-9}$
m
$^2$
s−1, closely matches that of ethanol in water
$(0.5{-}1.2) \times 10^{-9}$
m
$^2$
s−1 (see figure 11 below), enabling it to serve as an effective proxy for ethanol diffusion (Selifonov et al. Reference Selifonov, Shapoval, Mikerov and Tuchin2019). As shown in figure 1(c), stable stratification of the ethanol–water mixture was observed after the water injection. These stable layers, which can persist for tens of minutes or even several hours, exhibit distinct colour differences due to varying ethanol concentrations.
(a) Schematic of the generation of stratified ethanol–water mixtures by injecting water at the bottom of an ethanol layer in a tank. (b) Schematic diagram of the optical set-up for droplet imaging. (c) Evolution over 30 min of the stratified ethanol–water mixtures after injection of dyed water (20 ml) beneath pure ethanol (30 ml) at rate 2 ml min−1. (d) Calibration of the ethanol concentration using the colour intensity in the digital images. Vertical stripes measuring 100 pixels wide and 1000 pixels high with different ethanol concentrations
$w_e$
were used to determine the greyscale intensity in images of stratified ethanol–water mixtures. (e) Evolution of ethanol mass fraction
$w_e(z)$
over 30 min. The individual profiles are approximately linear for
$w_e \approx 10{-}90\,\%$
.

Figure 1. Long description
Panel A: A schematic diagram shows the setup for generating stratified ethanol-water mixtures by injecting water at the bottom of an ethanol layer in a tank. Panel B: A schematic diagram illustrates the optical setup for droplet imaging, including side and top views with light sources and a CCD camera. Panel C: A series of images displays the evolution of stratified ethanol-water mixtures over 30 minutes after injecting dyed water beneath pure ethanol. Each image is labeled with a specific time point ranging from 5 to 30 minutes. Panel D: A calibration graph shows the relationship between ethanol concentration and color intensity in digital images, with vertical stripes of different ethanol concentrations used to determine greyscale intensity. The horizontal axis is labeled with ethanol concentration in percent, and the vertical axis is labeled with greyscale intensity. Panel E: A line graph depicts the evolution of ethanol mass fraction over 30 minutes. The horizontal axis is labeled with distance in millimeters, and the vertical axis is labeled with ethanol mass fraction in percent. Different colored lines represent data at various time points, ranging from 5 to 30 minutes.
To calibrate the relation between the intensity of the dye and the ethanol mass fraction
$w_e$
, we followed Xue et al. (Reference Xue, Khodaparast, Zhu, Nunes, Kim and Stone2017). Digital images of water–ethanol mixtures with different ethanol concentrations were captured. Then an intensity–concentration calibration map was established, as shown in figure 1(d). Figure 1(e) shows concentration profiles at six different times estimated in this way. The concentration gradient
$\mathrm{d} w_e/\mathrm{d} z\approx 0.28$
mm
$^{-1}$
is nearly linear, and remains stable for approximately 30 min. This stability is verified by simulations based on Fick’s second law (see Appendix A for details). By replacing the initial pure ethanol with ethanol–water mixtures with
$w_e =55{-}90\,\%$
, we generated smaller gradients
$\mathrm{d} w_e/\mathrm{d} z\approx 0.14{-}0.25$
mm
$^{-1}$
, respectively. The deviation of the linear gradients is within 0.01 mm
$^{-1}$
over a 30 min period. Since droplet motions are observed for less than 5 min, the associated gradient error is estimated to be less than 0.01 mm
$^{-1}$
. Quantifying the concentration gradient with high spatial precision remains challenging due to its millimetre scale. To ensure reliability, we used the identical cuvette with illumination conditions identical to those applied during the experiments with droplets.
2.2. Experimental procedures
Following gradient generation, the mixture was left undisturbed for 15 min to allow convective flows to subside. During each experiment, a droplet of silicone oil (10 cSt, Sigma–Aldrich) was released into the cuvette using a micropipette (0.178 mm diameter) at a depth corresponding to the centre of the concentration gradient. By changing the injected volume, the droplet radius could be varied over a range
$\sim 400{-}800$
${\unicode{x03BC}}$
m. Only a single droplet was injected at a time to avoid interactions between droplets. Observations began immediately after droplet injection. A new gradient was established prior to each new trial. The droplet’s motion was recorded with a CCD camera (Zongyuan, ZY-PHC4K, 1920
$\times$
1080 pixel
$^2$
) in either side view or top view at 30 frames per second (fps). The experimental set-up for droplet imaging is shown in figure 1(b). Two CCD cameras equipped with long working distance lenses with magnification
$5\times$
were used, and two light sources provided illumination.
We observed the motion of each droplet for
$30{-}180$
s until it moved out of view of the camera. The trajectories of the droplets and their sizes were analysed using in-house MATLAB programs. Note that the standard deviation of the measured droplet radius (
$R\sim 400{-}800$
${\unicode{x03BC}}$
m) ranges from 2.8 to 4.7
${\unicode{x03BC}}$
m. Given the smallness of these deviations, error bars are not visibly distinguishable in the figures depicting the radius values. Therefore, we use throughout the time-averaged value of the droplet radius to represent the radius
$R$
.
2.3. Particle image velocimetry measurements
Polystyrene microparticles (Baseline Chrom. Tech., China) with diameter 2.5
${\unicode{x03BC}}$
m and concentration 0.2 mg ml−1 were added to the ethanol and water to visualise the flow fields. The particle-laden gradient was prepared in the same manner as described in § 2.1, and remained still for 10 min after the injection. Flow fields were then recorded using the CCD camera at 30 fps. No noticeable random fluid motion was observed, confirming that droplet dynamics occurred in an undisturbed background. By releasing a droplet into the gradient, we further captured the flow field surrounding the droplet. To further probe the internal flow field of the droplet, particles of the same diameter and concentration were dispersed in silicone oil. After generating the ethanol–water gradient without particles, a particle-laden oil droplet was released, and the internal flow field was recorded. The low tracer concentration used here made velocity extraction from particle image velocimetry (PIV) recordings challenging. Increasing the concentration improves signal quality but causes particles to accumulate at the droplet interface, obstructing interfacial flows and preventing lateral self-propulsion.
(a–f) Schematic illustration of MSD definition for a droplet’s lateral self-propulsion trajectory.

2.4. Characterisation of lateral self-propulsion
2.4.1. Definition of the MSD for time lags
To quantify droplet trajectories in the lateral direction, we analysed the MSD. Unlike the conventional definition based on displacements from the initial position, we adopted an alternative formulation in which the MSD is averaged over all time intervals within a single trajectory (Qian et al. Reference Qian, Sheetz and Elson1991; Michalet Reference Michalet2010).
As illustrated in figure 2(a), the trajectory is sampled at times
$\Delta t, 2\Delta t, \ldots , N\Delta t$
, where
$\Delta t$
is the fixed time interval, and
$N$
is the total number of frames. This yields
$N(N-1)/2$
forward displacements,
corresponding to time lags
$\Delta t_{ij} = (j-i)\Delta t$
, where
$\boldsymbol{r}_{\!j}$
and
$\boldsymbol{r}_i$
are the droplet displacements relative to the original position. The MSD is then defined as
\begin{align} \textit{MSD}(n\Delta t) = \frac {1}{N-n} \sum _{i=1}^{N-n} |\boldsymbol{d}_{i,i+n}|^2=\frac {1}{N-n} \sum _{i=1}^{N-n} |\boldsymbol{r}_{i\Delta t+n\Delta t} - \boldsymbol{r}_{i\Delta t}|^2, \quad n = 1, \ldots , N-1. \end{align}
As shown in figures 2(b–f), many displacement measurements are available for small time lags, but only a few for large time lags. Consequently, only the early portion of the MSD curve is statistically reliable and representative of typical droplet motion. In the following analysis, MSD data are presented up to time lag
$N\Delta t/2$
, ensuring sufficient averaging statistics.
2.4.2. Definition of lateral speed
The instantaneous droplet speed is defined as
where
$|\boldsymbol{d}_{i,i+1}|$
is the displacement between successive frames, and
$\Delta t$
is the corresponding time interval. The time-averaged speed magnitude is then given by
\begin{equation} U = \frac {1}{(N-1)\Delta t} \sum _{i=1}^{N-1} |\boldsymbol{d}_{i,i+1}|, \end{equation}
which is independent of instantaneous motion direction. In the subsequent analysis, the lateral self-propulsion speed
$U$
refers to this time-averaged speed. A sampling interval of
$(4/30)$
s was selected to balance temporal resolution and noise, thereby providing a representative measure of self-propulsion speed while avoiding excessive noise error. This definition emphasises overall propulsion strength rather than directional persistence, as Marangoni-driven droplet motion involves continuous reorientation.
3. Characteristics of lateral self-propulsion of droplets in stratified liquids
We start with an ethanol gradient
$\mathrm{d} w_e/\mathrm{d} z = 0.28$
mm
$^{-1}$
, and observe two typical behaviours. As seen in figures 3(a,b) (see also supplementary movies 1 and 2), when the droplet is small (
$R \sim 400$
${\unicode{x03BC}}$
m), it bounces vertically (figure 3
a), and undergoes lateral Brownian-like motion with a maximum displacement of a few hundred
${\unicode{x03BC}}$
m (figure 3
b). As
$R$
increases (
$R \sim 800$
${\unicode{x03BC}}$
m), the droplet becomes capable of propelling itself laterally to more distant regions (figures 3
c,d); supplementary movies 3 and 4). The green and red dashed circles in figures 3(a) and 3(c) show the initial positions of the droplets. The red dashed arrow in figure 3(c) indicates the droplet’s lateral displacement. The trajectories in figures 3(b) and 3(d) are coloured according to the instantaneous speed of the droplets. Figure 3(e) shows continuous bouncing trajectories
$z(t)$
of droplets with
$R = 378{-}760$
${\unicode{x03BC}}$
m in an ethanol gradient
$\mathrm{d} w_e/\mathrm{d} z = 0.28$
mm
$^{-1}$
. The trajectories are aligned with the maximum height
$z_{\textit{max}}$
to enhance visibility. Smaller droplets have larger amplitude and smaller frequency.
The lateral motion of the droplets can be characterised by the MSD in (2.2), shown in figure 3(f) as a function of time. The green and red curves correspond to the small and large droplets in figures 3(b) and 3(d), respectively. For the small droplets with radii 402 and 525
${\unicode{x03BC}}$
m (figure 3
b), the MSD increases approximately linearly with time
$t$
, i.e.
$r^2 \sim t$
. However, for the largest droplet,
$r^2 \sim t^2$
over 50 s (figure 3
d), indicating superdiffusive motion. We note that effective diffusivities of the small droplets are of the order of 10
$^{-8}$
m
$^2$
s−1, which exceeds the values for particles of comparable size undergoing thermal Brownian motion. This enhancement arises from weak Marangoni flows that provide the energy source to amplify stochastic displacements, yet remain insufficient to sustain directional self-propulsion. To distinguish this regime from cases of persistent self-propulsion for large droplets, we call it ‘Brownian-like motion’.
Dynamics of droplets in an ethanol–water stratification
$\mathrm{d}w_e/\mathrm{d}z = 0.28$
mm
$^{-1}$
from side view and top view. (a,b) Vertical and lateral motions of droplets with small radii (
$R = 378$
, 402
${\unicode{x03BC}}$
m). (c,d) Vertical and lateral motions of droplets with large radii (
$R = 760$
, 837
${\unicode{x03BC}}$
m). The small droplet exhibits vertical bouncing, and lateral Brownian-like motion, whereas larger droplets move laterally towards more distant regions. (e) Vertical droplet displacement versus time
$t$
for different droplet sizes. (f) The MSD of the lateral motions of droplets, with green and red curves for representive small and large droplets, respectively. (g) Schematic of Marangoni flow around an oil droplet.

The vertical bouncing due to an oscillatory Marangoni instability has been elucidated in previous studies (Li et al. Reference Li, Diddens, Prosperetti, Chong, Zhang and Lohse2019, Reference Li, Diddens, Prosperetti and Lohse2021). Because the silicone oil droplet (density 930 kg m−
$^{3}$
at
$25\,^{\circ}$
C) is less dense than water (997 kg m−
$^{3}$
at
$25\,^{\circ}$
C), but denser than ethanol (785 kg m−
$^{3}$
), after injection it first sinks towards its density-matched height. Meanwhile, the ethanol concentration gradient produces a surface tension gradient on the droplet interface because the interfacial tension of oil with the ethanol-rich liquid at the apex of the droplet is smaller than that with the ethanol-poor liquid at the bottom. Consequently, a downward Marangoni flow is induced that tends to lift the droplet (figure 3
g). When the Marangoni flow is stable, the droplet levitates at a fixed height governed by the balance of buoyancy and the Marangoni force. Droplets above a critical size bounce periodically due to the aforementioned oscillatory Marangoni instability (figure 3
a). The dependence of the bouncing height and the rising and sinking times on the droplet radius, the stratification strength and the drop viscosity have been determined by Meijer et al. (Reference Meijer, Li, Diddens and Lohse2023).
Here, we focus on lateral self-propulsion of droplets. Self-propulsion can be attributed to the balance between solutal advection and diffusion, with the characteristic ratio quantified by the Marangoni number
$\textit{Ma}$
, defined as
where
$D$
is the diffusivity of ethanol in water,
$V_M$
is the magnitude of the (downward) equatorial velocity of the Marangoni flow in the limit
$D\to \infty$
(Young et al. Reference Young, Goldstein and Block1959),
$\sigma$
is the surface tension coefficient, and
$\eta$
and
$\hat {\eta }$
are the dynamic viscosities of the ethanol–water mixture and the droplet, respectively. At low
$\textit{Ma}$
, diffusion dominates, resulting in the smoothing out of concentration gradients. However, above a critical value of
$\textit{Ma}$
, spontaneous symmetry breaking occurs in the concentration field surrounding the droplets, leading to self-propulsion. Symmetry-breaking-induced self-propulsion has been reported for droplets in surfactant solutions, as well as for catalytic particles with diffusiophoretic effects (Michelin et al. Reference Michelin, Lauga and Bartolo2013; Izri et al. Reference Izri, van der Linden, Michelin and Dauchot2014; Suga et al. Reference Suga, Suda, Ichikawa and Kimura2018; Chen et al. Reference Chen, Chong, Liu, Verzicco and Lohse2021; Michelin Reference Michelin2023). Suga et al. (Reference Suga, Suda, Ichikawa and Kimura2018) reported the self-propulsion of micron-sized nematic liquid crystal droplets in an aqueous surfactant solution, and observed a transition from random to helical to straight motion of the droplet over time. Michelin et al. (Reference Michelin, Lauga and Bartolo2013) showed that catalytic particles spontaneously break symmetry when the Péclet number (characterising the ratio between advection driven by diffusiophoresis and mass diffusion) exceeds a critical value. However, Michelin et al. (Reference Michelin, Lauga and Bartolo2013) involves diffusiophoretic flow, whereas we consider Marangoni flow induced by a surface tension gradient.
To determine which factors govern the lateral self-propulsion, we varied the concentration gradient by replacing the upper ethanol layer with an ethanol–water mixture with
$w_e = 60{-}90$
wt %, corresponding to stratification strengths
$\mathrm{d} w_e/\mathrm{d} z=0.16{-}0.25$
mm
$^{-1}$
, respectively. Figures 4(a–c) show representative cases for
$R = 773$
, 785 and 784
${\unicode{x03BC}}$
m and
$\mathrm{d} w_e/\mathrm{d} z = 0.28$
, 0.22 and 0.16 mm
$^{-1}$
. The simultaneous lateral self-propulsion and vertical bouncing when
$\mathrm{d} w_e/\mathrm{d} z = 0.28$
mm
$^{-1}$
is similar to that shown in figure 3(c). However, the droplets in gradients
$\mathrm{d} w_e/\mathrm{d} z = 0.22$
and 0.16 mm
$^{-1}$
display only lateral self-propulsion without vertical bouncing (figures 4
b,c).
Droplet motion in stratified ethanol–water mixtures with different ethanol concentration gradients
$\mathrm{d} w_e/\mathrm{d} z$
. (a) Simultaneous lateral self-propulsion and vertical bouncing of a droplet with
$R = 773$
${\unicode{x03BC}}$
m and
$\mathrm{d} w_e/\mathrm{d} z = 0.28$
mm
$^{-1}$
. (b,c) Lateral self-propulsion of droplets with
$R= 785$
${\unicode{x03BC}}$
m,
$\mathrm{d} w_e/\mathrm{d} z = 0.22$
mm
$^{-1}$
, and
$R = 784$
${\unicode{x03BC}}$
m,
$\mathrm{d} w_e/\mathrm{d} z = 0.16$
mm
$^{-1}$
, respectively. Red dashed circles mark the initial positions of the droplets. The red arrows indicate the direction of lateral displacement. (d) Vertical trajectories of droplets with different radii as functions of time
$t$
in three linear stratifications. (e) The MSD of the lateral motion of droplets with different radii and in different concentration gradients. The insets show the trajectories for the large droplets at different
$\mathrm{d} w_e/\mathrm{d} z$
. (f) Time-dependent speed
$U ^ {\prime }$
of droplets in a gradient
$\mathrm{d} w_e/\mathrm{d} z = 0.28$
mm
$^{-1}$
. (g) Lateral self-propulsion speed
$U$
of droplets as a function of droplet radius for different
$\mathrm{d} w_e/\mathrm{d} z$
. A selection of typical error bars is shown.

Figure 4. Long description
Panel A: A series of images showing the simultaneous lateral self-propulsion and vertical bouncing of a droplet with a radius of 378 micrometers and a concentration gradient of 0.28 millimeters inverse. Panel B: A series of images showing the lateral self-propulsion of a droplet with a radius of 402 micrometers and a concentration gradient of 0.22 millimeters inverse. Panel C: A series of images showing the lateral self-propulsion of a droplet with a radius of 396 micrometers and a concentration gradient of 0.16 millimeters inverse. Red dashed circles mark the initial positions of the droplets, and red arrows indicate the direction of lateral displacement. Panel D: A line graph showing the vertical trajectories of droplets with different radii as functions of time in three linear stratifications. The x-axis represents time in seconds, and the y-axis represents the vertical position in millimeters. Panel E: A log-log plot showing the mean squared displacement (MSD) of the lateral motion of droplets with different radii and in different concentration gradients. The x-axis represents the product of the number of steps and the time interval in seconds, and the y-axis represents the MSD in square millimeters. The insets show the trajectories for the large droplets at different time intervals. Panel F: A line graph showing the time-dependent speed of droplets in a gradient of 0.28 millimeters inverse. The x-axis represents time in seconds, and the y-axis represents the speed in millimeters per second. Panel G: A scatter plot showing the lateral self-propulsion speed of droplets as a function of droplet radius for different concentration gradients. The x-axis represents the droplet radius in micrometers, and the y-axis represents the speed in millimeters per second.
Overall, as depicted in figure 4(d), pronounced vertical bouncing occurs simultaneously for all droplets (
$R \approx 400{-}800$
${\unicode{x03BC}}$
m,
$\textit{Ma}\approx0.8\times 10{^5}$
–
$1.0\times 10{^6}$
) only for the largest gradient (
$\mathrm{d} w_e/\mathrm{d} z = 0.28$
mm
$^{-1}$
). For medium gradients (
$\mathrm{d} w_e/\mathrm{d} z = 0.25$
and 0.22 mm
$^{-1}$
), only the small (
$R \lesssim 450$
${\unicode{x03BC}}$
m,
$\textit{Ma}\lesssim 1.5 \times 10{^5}$
) droplets show obvious vertical bouncing. For large self-propelling droplets, the vertical bouncing is either of very small amplitude (less than the droplet radius) or non-existent (see supplementary movies 5 and 6). Under even lower gradients (
$\mathrm{d} w_e/\mathrm{d} z = 0.19$
and 0.16 mm
$^{-1}$
), no bouncing occurs for any droplet, i.e. small droplets levitate (
$\textit{Ma}\lesssim 1.3 \times 10{^5}$
), and large droplets self-propel horizontally without bouncing (
$\textit{Ma}\gtrsim 1.3 \times 10{^5}$
; see supplementary movies 7 and 8). For an even lower gradient
$\mathrm{d} w_e/\mathrm{d} z = 0.14$
mm
$^{-1}$
, no significant lateral self-propulsion was observed for any droplet size. We attribute this to a Marangoni flow that is too weak to drive lateral motion (see details in Appendix B).
The temporal MSD for the lateral motion with different
$\mathrm{d} w_e/\mathrm{d} z$
is shown in figure 4(e). Note that the horizontal axis,
$n\,\Delta t$
(
$\Delta t = 1/30\ \text{s}$
), represents the time-window size along the trajectory, rather than the starting time of the droplet motion. For small droplets with
$R \sim 400$
${\unicode{x03BC}}$
m, Brownian-like motion occurs for all gradients
$\mathrm{d} w_e/\mathrm{d} z$
, and the MSDs increase linearly with time. In contrast, large droplets (
$R\sim800$
${\unicode{x03BC}}$
m) exhibit self-propulsion with
$\textit{MSD}\sim t^2$
. Moreover, for
$\mathrm{d} w_e/\mathrm{d} z = 0.16$
mm
$^{-1}$
, the large droplet describes helical rather than random trajectories when viewed from above (figure 3
d). Similar helical paths of Marangoni-driven droplets have been reported previously (Suga et al. Reference Suga, Suda, Ichikawa and Kimura2018). The transition between helical and straight motion is governed by the self-propelling velocity, which is linked to the Ericksen number. The emergence of such complex trajectories in stratified environments highlights a promising avenue for future research, which we will explore in subsequent studies.
We also extracted the instantaneous lateral speed from these trajectories. Figure 4(f) shows representative time-varying speed
$U ^{\prime }$
for droplets in a gradient
$\mathrm{d} w_e/\mathrm{d} z = 0.28$
mm
$^{-1}$
, with standard deviation of
$U ^{\prime }$
ranging from 0.03 to 0.08 mm s−1. The time averaged
$U ^{\prime }$
is taken as the self-propulsion speed
$U$
in subsequent analysis.
As shown in figure 4(g),
$U$
increases with increasing droplet radius for all gradients tested. The lateral self-propulsion speed
$U$
normalised by
$R$
as a function of
$R$
for
$\mathrm{d} w_e/\mathrm{d} z = 0.16$
, 0.22 and 0.28 mm
$^{-1}$
is shown in Appendix C. The ratio
$U/R$
for each value of
$\mathrm{d} w_e/\mathrm{d} z$
, while scattered, shows only a small trend as a function of
$R$
, which suggests that
$U$
scales as
$R$
. Moreover, we find that the droplet radius remains approximately constant throughout the observation period, indicating that the mass transfer from the droplet to the ethanol–water mixture is negligible. These sustained motions are thereby driven by the vertical ethanol concentration gradient. Notably, the velocity shows no measurable decay along the trajectory. This constancy is attributed to the short recording duration before droplets exit the camera’s field of view, as well as their migration into previously unperturbed regions of the gradient where they continually access fresh driving potential.
4. Scaling analysis for the onset and steady-state amplitude of lateral motion
We now use theoretical scaling analysis to gain further insight into the onset of lateral self-propulsion, focusing on the case of negligible vertical bouncing. Our analysis starts from the recognition of Li et al. (Reference Li, Diddens, Prosperetti and Lohse2021) that departures of the ethanol concentration from the far-field linear gradient are confined to a boundary layer adjoining the droplet, whose thickness
$\delta$
is governed by a balance of diffusion (which increases
$\delta$
) and advection (which decreases
$\delta$
) of ethanol.
Before beginning our scaling analysis, we need to demonstrate that the Marangoni velocity scale
$V_M$
, which was derived in the limit of infinite diffusivity (
$\textit{Ma} = 0$
) by Young et al. (Reference Young, Goldstein and Block1959), remains the appropriate velocity scale for the case
$\textit{Ma}\gg 1$
that interests us. This is indeed the case. The reason is that the gradient of the ethanol concentration along the droplet’s surface (which produces the Marangoni stress that drives flow in the outer fluid) differs by only 50 % between the cases
$\textit{Ma}\gg 1$
and
$\textit{Ma} = 0$
. Let
$(\mathrm{d} w_e/\mathrm{d}z)_{\infty }\equiv \beta$
be the ethanol concentration gradient far from the droplet. For the case
$\textit{Ma}\gg 1$
, the far-field region extends right up to the edge of the thin boundary layer, so the vertical gradient of ethanol concentration that drives Marangoni flow is approximately
$\beta$
. For the case
$\textit{Ma}=0$
, by contrast, the ethanol concentration field satisfies Laplace’s equation. The solution that satisfies the boundary conditions (matching to the far field and zero flux across the droplet’s surface) is
$w_e = \beta (r + R^3/2 r^2)\cos \theta$
, which implies that the vertical ethanol concentration gradient at the droplet’s surface is
$3\beta /2$
– just 50 % larger than for the case
$\textit{Ma}\gg 1$
. The fact that the ethanol concentration gradients at the droplet’s surface differ by only an
$O(1)$
numerical factor between the cases
$\textit{Ma}\gg 1$
and
$\textit{Ma} =0$
means that
$V_M$
remains the appropriate velocity scale when
$\textit{Ma}\gg 1$
.
Now that the appropriateness of the velocity scale
$V_M$
has been demonstrated, we recall the simple scaling analysis by which Li et al. (Reference Li, Diddens, Prosperetti and Lohse2021) determined the boundary-layer thickness
$\delta$
in the limit
$\textit{Ma}\gg 1$
. The balance of advection and diffusion in this limit requires
which confirms that the boundary layer is thin if
$\textit{Ma} \gg 1$
.
4.1. Onset of the lateral motion
The PIV measurements reveal that instability originates from breaking of the axisymmetry of the concentration field. The external interfacial flow (see supplementary movie 9;
$\mathrm{d} w_e/\mathrm{d} z = 0.22$
mm
$^{-1}$
,
$R = 657$
${\unicode{x03BC}}$
m) shows stronger Marangoni stresses on the left-hand side of the droplet than on the right-hand side, driving the droplet’s rightward motion. No rotational motion of the droplet is observed from its internal flow field (see supplementary movie S10;
$\mathrm{d} w_e/\mathrm{d} z = 0.22$
mm
$^{-1}$
,
$R = 810$
${\unicode{x03BC}}$
m), suggesting that droplet self-propulsion arises from lateral net Marangoni stress. To quantify this instability, we consider a stationary droplet that is perturbed by moving it laterally to the right (say) with an infinitesimal velocity
$\hat U$
. Two characteristic time scales can be identified. The advective time scale is the time
$\tau _a\sim \delta /\hat U$
for the droplet to move a distance equal to the boundary-layer thickness
$\delta$
. The diffusive time scale
$\tau _d\sim \delta ^2/D$
is the time required for the perturbed boundary layer to regain its fore–aft symmetry by diffusion. Instability sets in if
$\tau _d \gt c_1 \tau _a$
, where
$c_1$
is a constant. Thus we have
Now the Marangoni shear stress in the boundary layer is
$\sigma _{r\theta }\sim \eta V_M/\delta$
(Li et al. Reference Li, Diddens, Prosperetti and Lohse2021). By a well-known property of isoviscous Stokes flow around a sphere, the normal stress
$\sigma _{rr}$
is of the same order of magnitude as
$\sigma _{r\theta }$
. The lateral Marangoni force on the droplet must be balanced by the Stokes drag, which requires
Eliminating the product
$\hat U\delta$
between (4.2) and (4.3), we find the instability criterion
Using the definition (3.1) of
$V_M$
, we rewrite (4.4) to obtain the critical droplet radius
$R_{\textit{crit}}$
above which lateral self-propulsion occurs:
\begin{equation} R_{\textit{crit}} = \sqrt {\frac {2 c_1 D (\eta + \hat \eta )}{-(\mathrm{d}\sigma /\mathrm{d} w_e)(\mathrm{d} w_e/\mathrm{d} z)}}. \end{equation}
Figure 5 shows the radius
$R$
of droplets exhibiting lateral motion (red triangles) and no lateral motion (black squares) for five ethanol gradients
$\mathrm{d} w_e/\mathrm{d} z$
in the range 0.16–0.28 mm
$^{-1}$
. For the lowest four of these, least squares fitting gives
$R_{\textit{crit}} \propto (\mathrm{d} w_e/\mathrm{d} z)^{-0.45 \pm 0.11}$
, which is consistent with the prediction in (4.5). However, this scaling fails for
$\mathrm{d} w_e/\mathrm{d} z = 0.28$
mm
$^{-1}$
. We attribute this to the large-amplitude (
$\gtrsim R$
) vertical bouncing that accompanies the lateral motion for this gradient, which is not consistent with the assumption of the scaling analysis that no bouncing occurs. We speculate that for large gradients, the large-amplitude oscillations entrain and drag a portion of the surrounding liquid, thereby modifying the original local concentration field, which may substantially influence the droplet dynamics. A full investigation of the complex hydrodynamics at larger gradients is beyond the scope of the present study. Note that to clearly distinguish the criteria for ‘no lateral motion’ and ‘lateral motion’, the vertical axis in figure 5 is limited to 600
${\unicode{x03BC}}$
m. In our experiments, droplets with radii beyond this threshold, up to
$\sim 900$
${\unicode{x03BC}}$
m, all exhibit lateral motion.
Radius
$R$
of droplets exhibiting lateral motion (red triangles) and no lateral motion (black squares) for ethanol concentration gradient
$\mathrm{d} w_e/\mathrm{d} z = 0.16{-}0.28$
mm
$^{-1}$
(double logarithmic plot).

4.2. Steady-state amplitude
The scaling law for the steady-state lateral speed
$U$
beyond onset can be obtained directly from the force balance (4.3). For a Péclet number
$\textit{Pe} = \textit{UR}/D \gg 1$
(e.g.
$\textit{Pe} = 100$
for a droplet of radius 500
${\unicode{x03BC}}$
m with speed
$U = 0.2$
mm s−1), the boundary layer on the back (downstream) side of the droplet, which drives the droplet’s motion, has thickness
$\delta \sim R$
. Substituting this expression into (4.3), we obtain
where
$c_2$
is a constant. Only the back of the boundary layer need be considered here, because the front exerts a normal force on the droplet that is
$\textit{Pe}^{-1/2}$
times that exerted by the back (see Appendix D).
Figure 6(a) shows
$U$
versus
$V_M$
for the experiments with
$\mathrm{d} w_e/\mathrm{d} z = 0.16{-}0.25$
mm
$^{-1}$
, where
$V_M$
was calculated as described in Appendix F. The experiments with
$\mathrm{d} w_e/\mathrm{d} z = 0.28$
mm
$^{-1}$
, as well as a subset at
$\mathrm{d} w_e/\mathrm{d} z = 0.19{-}25$
mm
$^{-1}$
, are omitted because of the aforementioned large-amplitude bouncing that occurs in that case. The solid line represents the best-fit power law
$U \propto V _M ^\beta$
with exponent
$\beta = 1.11 \pm 0.23$
obtained from least squares fitting to the logarithmic data, and is consistent with the prediction of (4.6). The reported error corresponds to the standard error of the slope with a 95 % confidence interval. The ratio
$U/V_M$
as a function of droplet radius
$R$
exhibits noticeable scatter and a slight decrease at larger
$R$
(figure 6
b). However, this dependence is considerably weaker than the scaling of
$U$
with
$V_M$
, thereby reinforcing the dominant relation
$U\propto V_M$
.
(a) Lateral self-propulsion speed
$U$
of the droplet as a function of the Marangoni velocity scale
$V_M$
(3.1) for ethanol concentration gradients
$\mathrm{d} w_e/\mathrm{d} z = 0.16{-}0.25$
mm
$^{-1}$
on a double logarithmic plot. The grey line is
$U \propto V_M ^{1.11 \pm 0.23}$
, obtained by a least squares fit to the data. (b) Plot of
$U/V_M$
as a function of droplet radius
$R$
.

Figure 6. Long description
Two scatter plots depict the relationship between droplet speed, Marangoni velocity, and droplet radius. Panel A: A scatter plot shows the lateral self-propulsion speed of the droplet as a function of the Marangoni velocity scale for ethanol concentration gradients on a double logarithmic plot. The x-axis represents the Marangoni velocity scale in mm per second, and the y-axis represents the lateral self-propulsion speed in mm per second. Different colored circles represent various ethanol concentration gradients. A grey line obtained by a least squares fit to the data is also shown. Panel B: A scatter plot shows the ratio of the lateral self-propulsion speed to the Marangoni velocity scale as a function of droplet radius. The x-axis represents the droplet radius in micrometers, and the y-axis represents the ratio of the lateral self-propulsion speed to the Marangoni velocity scale. Different colored squares represent various ethanol concentration gradients.
Phase diagram for droplet motion in the space of
$\textit{Ma}/Ra ^{1/2}$
versus
$Ra$
. Black triangles denote stable droplets, and red circles from Li et al. (Reference Li, Diddens, Prosperetti and Lohse2021) represent bouncing droplets. The horizontal line between them indicates the critical value of
$\textit{Ma}/Ra ^{1/2}$
for bouncing. The yellow rectangles, green diamonds and blue triangles correspond to droplets exhibiting only bouncing, bouncing with lateral motion, and only lateral motion, respectively, in our experiments. The black solid line with slope
$-1/2$
represents the critical value of
$\textit{Ma}$
for lateral self-propulsion.

4.3. Regime diagram
By combining our data with that of Li et al. (Reference Li, Diddens, Prosperetti and Lohse2021), we can construct a regime diagram (figure 7) of the droplet dynamics in the space of
$\textit{Ma}/Ra ^{1/2}$
versus
$Ra$
, where
is the Rayleigh number. Our method for calculating
$\textit{Ma}$
and
$Ra$
for each data point is described in Appendix F. As shown by Li et al. (Reference Li, Diddens, Prosperetti and Lohse2021), stable levitating drops undergo an oscillatory Marangoni instability to periodic bouncing when
$\textit{Ma}/Ra^{1/2} \gt 275$
. In addition, our experiments and scaling analysis show that lateral motion occurs when
$\textit{Ma} \gtrsim (1.3 \pm 0.2) \times 10^5$
(line with slope
$-1/2$
in figure 7). Note that our data for bouncing droplets without lateral motion (yellow rectangles) overlap with the data of Li et al. (Reference Li, Diddens, Prosperetti and Lohse2021) for the same behaviour (red circles). This agreement highlights the consistency between the present observations and prior studies, while our work extends this to the regime of lateral self-propulsion.
5. Discussion
The lateral self-propulsion of silicone oil droplets in stratified ethanol–water mixtures has been investigated by tracking the droplet dynamics from both side and top views. The characteristic droplet motion changes from vertical periodic bouncing to lateral self-propulsion as the droplet size increases, with an additional dependence on the magnitude of the ethanol concentration gradient. In some cases, low-amplitude bouncing and lateral translation occur simultaneously. When the concentration gradient is relatively small (0.16 mm
$^{-1}$
), large droplets display only lateral self-propulsion without bouncing. Scaling analysis shows that the transition from Brownian-like motion to self-propulsion occurs when the Marangoni number
$\textit{Ma}$
exceeds a critical value
$\textit{Ma}_{\textit{crit}}\approx 10^5$
. The onset of self-propulsion reflects spontaneous symmetry breaking of axisymmetric Marangoni convection around the droplet. The criterion
$\textit{Ma} \gt \textit{Ma}_{\textit{crit}}$
implies that the critical radius for lateral motion decreases as the
$-1/2$
power of the concentration gradient, a prediction that is consistent with our experiments for four of the five ethanol gradients used. Furthermore, those experiments also show that the speed of steady lateral self-propulsion increases approximately linearly with the characteristic Marangoni velocity scale, supporting a further prediction of the scaling analysis. Analysis of the evolution of the complex droplet trajectories beyond onset is a promising direction for research, which we will explore in future.
In the phase diagram of
$\textit{Ma}/Ra^{1/2}$
versus
$Ra$
(figure 7), the blue triangles denoting purely lateral motion overlap with the green diamonds corresponding to lateral motion with bouncing. The reason for this overlap remains unclear and warrants further investigation in subsequent studies. One possibility is that it represents the projection onto the
$Ra$
–
$\textit{Ma}/Ra^{1/2}$
plane of a more complicated three-dimensional structure involving a third dimensionless parameter along an axis perpendicular to the plane. Two possible such parameters are the viscosity ratio and the aspect ratio
$R/L$
, where
$L$
is the container size. However, the viscosity ratio was nearly constant in our experiments, and it seems unlikely that variations in the aspect ratio by a factor
$\approx 2$
could have much effect given that
$R/L\ll 1$
in all cases. Another possibility is that the degree of coexistence of lateral self-propulsion and vertical bouncing depends on the relative magnitudes of initial perturbations in the droplet’s horizontal and vertical positions. A promising way to resolve this question is via direct numerical simulations, which we have already begun to build.
The reader may also have noted that a few red circles in figure 7, representing pure bouncing, fall below the line separating bouncing and stable phases at large
$Ra \approx 100$
. We attribute this to the confinement effect of the container. In prior work of Li et al. (Reference Li, Diddens, Prosperetti and Lohse2021), the instability threshold
$c \approx 275$
was validated experimentally for
$6 \times 10^{-3} \lesssim Ra \lesssim 3$
. Experiments at
$Ra \lt 6 \times 10^{-3}$
are inaccessible because droplets with
$R \lt 20$
${\unicode{x03BC}}$
m are too small to observe, while for
$Ra \gt 3$
, the finite size of the container comes into play. It is speculated that
$c \approx 273$
remain valid only if the container size
$L$
is large enough, with
$L \gt \mathcal{L}$
, where
$\mathcal{L} \sim (-(\mathrm{d} \rho / \mathrm{d} z)\,\mu D / g)^{-1 / 4}$
(Phillips Reference Phillips1970). The presence of points below the dashed line thus indicates an altered instability criterion induced by confinement. Again, numerical simulation is a promising way to extend the parameter space beyond current experimental limitations to study these rich dynamics.
The instability that we have studied is not limited to concentration gradients, but can also be extended to other stratified systems involving e.g. temperature or surfactant gradients (Izzet et al. Reference Izzet, Moerman, Gross, Groenewold, Hollingsworth, Bibette and Brujic2020; Hu et al. Reference Hu, Li, Chen, Wu, Huo, Gu, Wang and Deng2024, Reference Hu, Wang, Li, Chen, Wu, Huo and Deng2025), and vertically vibrating fluids (Bush Reference Bush2015; Guan et al. Reference Guan, Tamim, Magoon, Stone and Sáenz2025). For instance, it is found that a bubble can exhibit spontaneous ‘hovering’ along horizontal surfaces under stationary laser heating (Hu et al. Reference Hu, Li, Chen, Wu, Huo, Gu, Wang and Deng2024). Self-propulsion due to spontaneous breaking of axisymmetry is thus robust and generic. This rich dynamics suggests promising applications in microfluidic transport, sensing and cleaning, and drug delivery as well as in biological systems, where slow but persistent directional motion plays a crucial role.
Supplementary movies
Supplementary movies are available at https://doi.org/10.1017/jfm.2026.11828.
Acknowledgements
The authors thank Professor D. Lohse, Professor H. Wyss, Professor Y. Wang and Dr B. Zeng for valuable discussions. We are grateful to three anonymous referees whose suggestions led to substantial improvements in the original manuscript.
Funding
This work was supported in part by the Research Impact Fund (R4015-21) and the Health@InnoHK programme of the Innovation and Technology Commission of the Hong Kong SAR Government. H.C.S. was supported in part by the National Natural Science Foundation of China (no. 22525805) and the RGC Senior Research Fellowship (SRFS2425-7S04). N.M.R. was supported by the Horizon Europe ERC grant 101098375 SOFT-PLANET. C.W. was supported in part by the National Natural Science Foundation of China (no. 22405226).
Declaration of interests
The authors declare the following competing interests: Ho Cheung Shum holds shares in, acts as Director of, or acts as advisor of MicroDiagnostics Limited, PharmaEase Tech Limited, Upgrade Biopolymers Limited, Monexus Innovation Limited, Multera Inc., and EN Technology Limited, and he is the Founding Centre Director and Co-Director of Advanced Biomedical Instrumentation Centre Limited. He is also a member of the Board of Directors of Advanced Biomedical Instrumentation Centre Limited and Hong Kong Institution of Sciences Limited. The work in this paper is, however, not directly related to the works of these entities, as far as we know. The authors declare no other competing interests.
Appendix A. Simulation of evolution of the ethanol–water gradient
The mass transport of ethanol in water is described by Fick’s second law extended for concentration-dependent diffusivity:
where
$w(z,t)$
represents the ethanol mass fraction,
$t$
is time,
$z$
is the spatial coordinate, and
$D(w)$
is the concentration-dependent mutual diffusion coefficient. To capture the temporal evolution of the ethanol–water concentration gradient, we performed numerical simulations based on (A1). Experimental mutual diffusivity data for ethanol–water at
$25\,^{\circ}$
C were compiled from the literature (Par et al. Reference Par2013). A shape-preserving piecewise cubic Hermite interpolating polynomial was used to construct a continuous
$D(w)$
function. The diffusion process was simulated using a one-dimensional finite difference scheme based on Fick’s second law (A1). The system was initialised with a sharp interface, defined by the initial condition
\begin{equation} w(z, 0)=\begin{cases} 1, & z\lt 0 \text{ (pure ethanol)}, \\ 0, & z \geqslant 0\text{ (pure water)}, \end{cases} \end{equation}
and boundary conditions
with domain of length
$L = 20\rm\, mm$
. The spatial domain
$z \in [-L/2, L/2]$
is discretised using
$N = 501$
equally spaced points (
$\Delta z = 40$
${\unicode{x03BC}}$
m). Time integration is performed with step
$\Delta t = 0.1$
s over a total duration of 300–1800 s. Equation (A1) is discretised using the second-order Crank–Nicolson scheme.
The simulated evolution of ethanol mass fraction due to diffusion reproduces the experimental trends (figure 8). Exact agreement between experiment and simulation is not expected, since convection during injection is unavoidable even at low flow velocities, whereas the model assumes purely diffusive transport.
Experimental and numerical profiles of ethanol mass fraction
$w_e(z)$
at 5, 15 and 25 min. The profiles remain approximately linear for
$w_e \approx 10{-}90\,\%$
over 25 min.

Appendix B. Termination of lateral motion for a low concentration gradient (
$\mathrm{d} w_e/\mathrm{d} z = 0.14$
mm
$^{-1}$
)
In the smallest ethanol concentration gradient
$\mathrm{d} w_e/\mathrm{d} z = 0.14$
mm
$^{-1}$
, the lateral displacement of droplets of different sizes is smaller than their diameter over 5 min (figure 9), indicating no significant lateral self-propulsion. We attribute this to the Marangoni flow being too weak to drive lateral motion.
Droplet trajectories in ethanol–water gradient with
$\mathrm{d} w_e/\mathrm{d} z = 0.14$
mm
$^{-1}$
. The displacement of droplets is less than
$2R$
over 5 min, exhibiting no obvious lateral self-propulsion.

Appendix C. Lateral self-propulsion speed normalised by droplet radius in different linear stratifications
Figure 10 shows the lateral speed
$U$
of the droplets normalised by
$R$
as a function of
$R$
for
$\mathrm{d} w_e/\mathrm{d} z = 0.16$
, 0.22 and 0.28 mm
$^{-1}$
. The ratio
$U/R$
increases with
$\mathrm{d} w_e/\mathrm{d} z$
, and exhibits no obvious dependence on
$R$
at 0.16 and 0.22 mm
$^{-1}$
. At
$\mathrm{d} w_e/\mathrm{d} z = 0.28$
mm
$^{-1}$
,
$U/R$
shows a slight increase with
$R$
, which is far weaker than that of
$U$
with
$R$
in figure 4(g).
Lateral self-propulsion speed
$U$
normalised by radius
$R$
as a function of
$R$
for
$\mathrm{d} w_e/\mathrm{d} z = 0.16$
, 0.22 and 0.28 mm
$^{-1}$
.

Appendix D. Normal stresses exerted by the front and back boundary layers
In this appendix, we show that for large
$Pe = RU/D$
, the back (downstream) side of the boundary layer provides the dominant contribution to the normal stress that balances the Stokes drag on the moving droplet. Let
$\delta _-\sim R$
and
$\delta _+\sim R\textit{Pe}^{-1/2}$
be the characteristic thicknesses of the boundary layer on the back and front of the droplet, respectively. Further, let
$V_M^-$
and
$V_M^+$
be the characteristic Marangoni velocities in the back and front of the boundary layer, respectively. The balance of the shear stress and the integrated buoyancy inside the boundary layer requires (Li et al. Reference Li, Diddens, Prosperetti and Lohse2021)
where
$\Delta \rho$
is the density difference between the boundary layer and the far field. With the definitions of
$\delta _-$
and
$\delta _+$
given above, (D1) implies
Now, the lateral forces applied by the back and front of the boundary layer are
$F_{\pm }\sim R^2\eta V_M^\pm /\delta _\pm$
. Thus
showing that the contribution from the back of the boundary layer is dominant.
Appendix E. Physical properties of ethanol–water mixtures and silicone oil
The physical properties of the ethanol–water mixture vary depending on the weight fraction of ethanol
$w_e$
. The density and dynamic viscosity data for the mixture are obtained from Khattab et al. (Reference Khattab, Bandarkar, Fakhree and Jouyban2012). The diffusivity of ethanol within the mixture is taken from Par et al. (Reference Par2013). Interfacial tension
$\sigma$
between 10 cSt silicone oil and the ethanol–water mixture is taken from Li et al. (Reference Li, Diddens, Prosperetti and Lohse2021). The density, dynamic viscosity, diffusivity and surface tension are shown as functions of
$w_e$
in figure 11.
Physical properties of ethanol–water mixtures as functions of the ethanol weight fraction
$w_e$
. (a,b) Density and dynamic viscosity at
$20\,^{\circ}$
C, respectively (Khattab et al. Reference Khattab, Bandarkar, Fakhree and Jouyban2012). (c) Plot of diffusivity
$D$
. Blue data points are from Par et al. (Reference Par2013), and the black solid line is a polynomial fit. (d) Interfacial tension
$\sigma$
between silicone oil and the ethanol–water mixture (Li et al. Reference Li, Diddens, Prosperetti and Lohse2021).

Appendix F. Estimation of Marangoni velocity
$V_M$
and Marangoni number
$\textit{Ma}$
To calculate the Marangoni velocity
$V_M$
and the Marangoni number
$\textit{Ma}$
, we need a value of the ethanol concentration
$w_e$
for each droplet. We estimate
$w_e$
based on the linear stratification and the sizes of the droplets. The experimental images show that small droplets are located at greater heights
$z$
than larger droplets (
$z$
values corresponding to the lowest positions are used for bouncing drops), consistent with previous studies (Li et al. Reference Li, Diddens, Prosperetti and Lohse2021). We take the position of droplets with
$R\gtrsim800$
${\unicode{x03BC}}$
m as the density-matched height where
$w_{e0} = 0.4$
. The ethanol concentration is then calculated as
$w_e = w_{e0} + (\mathrm{d} w_e/\mathrm{d} z)(\mathrm{d} z/\mathrm{d} R) (800-R)$
, with
$\mathrm{d} z/\mathrm{d} R = 1.67$
obtained by fitting a line to
$z$
versus
$R$
from experimental measurements. According to this equation, the ethanol concentration values for drops with radius
$R = 400{-}800$
${\unicode{x03BC}}$
m are
$w_e$
0.40–0.50, 0.40–0.54 and 0.40–0.59 for experiments with
$\mathrm{d} w_e/\mathrm{d} z = 0.16$
, 0.22 and 0.28 mm
$^{-1}$
, respectively. Finally, using the physical properties corresponding to the calculated value of
$w_e$
(figure 11), we determine the interfacial tension
$\sigma$
, the viscosity
$\eta$
, and the diffusivity
$D$
that are needed to calculate
$V_M$
and
$\textit{Ma}$
.





we
we(z)
we≈10−90%
dwe/dz=0.28
−1
R=378
μ
R=760
μ
t
dwe/dz
R=773
μ
dwe/dz=0.28
−1
R=785
μ
dwe/dz=0.22
−1
R=784
μ
dwe/dz=0.16
−1
t
dwe/dz
U′
dwe/dz=0.28
−1
U
dwe/dz
R
dwe/dz=0.16−0.28
−1
U
VM
dwe/dz=0.16−0.25
−1
U∝VM1.11±0.23
U/VM
R
Ma/Ra1/2
Ra
Ma/Ra1/2
−1/2
Ma
we(z)
we≈10−90%
dwe/dz=0.14
−1
2R
U
R
R
dwe/dz=0.16
−1
we
20∘
D
σ