3.2.1. Governing equations and force–velocity relation

We develop a theoretical framework to relate the translational motion of an evaporating sessile droplet to the balance of capillary force arising under vapour-mediated interactions and viscous drag force. The analysis is restricted to the quasistatic regime relevant to the present experiments, in which inertial effects are negligible and the droplet shape adjusts rapidly compared with the time scale of translation. Under these conditions, the flow inside the droplet is governed by incompressible Stokes equations. The velocity field $\boldsymbol{u}(\boldsymbol{x})$ and pressure $p(\boldsymbol{x})$ satisfy $\boldsymbol{\nabla }\boldsymbol{\cdot }\boldsymbol{\sigma } = \boldsymbol{0}$ and $\boldsymbol{\nabla }\boldsymbol{\cdot }\boldsymbol{u} = 0$ , where the Newtonian stress tensor is given by $\boldsymbol{\sigma } = -p \boldsymbol{I} + \mu [ \boldsymbol{\nabla }\boldsymbol{u} + (\boldsymbol{\nabla }\boldsymbol{u})^{T} ]$ .

The liquid–solid interface is modelled using a Navier-slip boundary condition, $\boldsymbol{u}_s = l_s {\boldsymbol{\tau }_s}/{\mu }$ where $l_s$ is the slip length and $\boldsymbol{\tau }_s$ is the tangential stress at the liquid–solid interface $S_{ls}$ . On the liquid–air interface $S_{lg}$ , the normal stress is balanced by capillarity $\boldsymbol{\sigma } \boldsymbol{\cdot }\boldsymbol{n} = \gamma _d \kappa \boldsymbol{n}$ , where $\boldsymbol{n}$ is the outward unit normal, $\gamma _d$ is the surface tension of the bulk droplet and $\kappa$ is the mean curvature. Numerical simulations of evaporating binary droplets show that surface tension gradients are strongly localized within a narrow region near the contact line (Karpitschka et al. Reference Karpitschka, Liebig and Riegler2017; Diddens et al. Reference Diddens, Kuerten, Van der Geld and Wijshoff2017a ; Dekker et al. Reference Dekker, Rocha, Lohse and Diddens2026). In this region, the surface tension gradient follows a power law dependence on the distance $s$ to the contact line, scaling as $\sim s^{-3/2}$ (Ramírez-Soto & Karpitschka Reference Ramírez-Soto and Karpitschka2022). Following previous analyses of vapour-mediated droplet interactions, the system is described using a simplified picture of a droplet consisting of a bulk phase with a single, spatially uniform composition, surrounded by a nanometric precursor film whose composition may vary along the substrate. Within the precursor film, local phase equilibrium at the liquid–air interface constrains the liquid composition to adjust to the vapour partial pressure above the surface. Spatial variations in vapour concentration therefore alter the local vapour partial pressure, which shifts the equilibrium liquid composition and hence the surface tension within the film (Wang et al. Reference Wang, Karapetsas, Valluri, Sefiane, Williams and Takata2021; Chao et al. Reference Chao, Ramírez-Soto, Bahr and Karpitschka2022; Zhao et al. Reference Zhao, Orejon, Sefiane and Shanahan2025). Because the surface stresses are confined to a region small compared with the droplet size, their integrated effect at macroscopic scales can be effectively represented by a line force acting at the apparent contact line, consistent with Cira et al. (Reference Cira, Benusiglio and Prakash2015) and Zhao et al. (Reference Zhao, Orejon, Sefiane and Shanahan2025). When this effective line force becomes azimuthally unbalanced, a net driving force arises and the droplet undergoes motion.

To relate the translational velocity of the droplet to this net force, we employ the Lorentz reciprocal theorem for incompressible Stokes flows (Happel & Brenner Reference Happel and Brenner2012; Masoud & Stone Reference Masoud and Stone2014, Reference Masoud and Stone2019). We define an auxiliary flow $(\hat {\boldsymbol{u}},\hat {\boldsymbol{\sigma }})$ in the same instantaneous droplet geometry, corresponding to the rigid-body translation of the droplet with prescribed velocity $\hat {\boldsymbol{U}}$ along the substrate, sharing the same direction with the actual droplet motion. The boundary conditions are chosen to satisfy the following principles: on the liquid–solid interface $S_{ls}$ , the auxiliary flow obeys the same Navier-slip condition as the actual flow, namely $\hat {\boldsymbol{u}}_s = l_s {\hat {\boldsymbol{\tau }}_s}/{\mu }$ , to ensure that stress integrals at the contact line converge; on the liquid–air interface $S_{lg}$ , the auxiliary velocity is prescribed as a rigid-body translation, $\hat {\boldsymbol{u}} = \hat {\boldsymbol{U}}$ , thereby isolating the translational degree of freedom in the Lorentz reciprocal theorem analysis; finally, the auxiliary flow vanishes in the far field, $\hat {\boldsymbol{u}} \to 0$ as $ r \to \infty$ , so that the surface integrals form a closed system.

In the auxiliary problem, the surface tension is uniform and the interface shape is held fixed. This choice is justified by the quasistatic nature of the dynamics, for which shape relaxation occurs on time scales much shorter than translational motion. The auxiliary problem therefore isolates the translational degree of freedom while preserving the hydrodynamic boundary conditions.

The Lorentz reciprocal theorem states that

(3.1) \begin{align} \int _S \hat {\boldsymbol{u}} \boldsymbol{\cdot }\left ( \boldsymbol{\sigma } \boldsymbol{\cdot }\boldsymbol{n} \right ) \, \mathrm{d}S = \int _S \boldsymbol{u} \boldsymbol{\cdot }\left ( \hat {\boldsymbol{\sigma }} \boldsymbol{\cdot }\boldsymbol{n} \right ) \, \mathrm{d}S, \end{align}

where $S=S_{lg}\cup S_{ls}$ is a closed surface enclosing the fluid domain. On the liquid–air interface $S_{lg}$ , $\int _{S_{lg}} \boldsymbol{u} \boldsymbol{\cdot }( \hat {\boldsymbol{\sigma }} \boldsymbol{\cdot }\boldsymbol{n} ) \, \mathrm{d}S=\boldsymbol{U}\boldsymbol{\cdot }\int _{S_{lg}} ( \hat {\boldsymbol{\sigma }} \boldsymbol{\cdot }\boldsymbol{n} ) \, \mathrm{d}S$ , where $\boldsymbol{U}$ represents the translational velocity of the droplet’s centre of mass. On the liquid–solid interface $S_{ls}$ , both $\boldsymbol{u}$ and $\hat {\boldsymbol{u}}$ satisfy the same Navier-slip condition, so contributions to the surface integrals are balanced, $\int _{S_{ls}} \boldsymbol{u} \boldsymbol{\cdot }\hat {\tau }_s \text{d}S = \int _{S_{ls}} \boldsymbol{u} \boldsymbol{\cdot }({\mu } \hat {\boldsymbol{u}}/{l_s} ) \text{d}S=\int _{S_{ls}} \hat {\boldsymbol{u}} \boldsymbol{\cdot }\tau _s \text{d}S$ .

The integral reduces to

(3.2) \begin{align} \hat {\boldsymbol{U}} \boldsymbol{\cdot }\boldsymbol{F} = \boldsymbol{U} \boldsymbol{\cdot }\hat {\boldsymbol{F}}_d, \end{align}

where $\boldsymbol{F} = \int _{S_{lg}} (\boldsymbol{\sigma } \boldsymbol{\cdot }\boldsymbol{n}) \text{d}S=\int _C \gamma \boldsymbol{t} \, \text{d}l$ is the net capillary driving force concentrated along the contact line, $\boldsymbol{t}$ being the unit vector tangent to the substrate.

Under the small contact angle condition, the velocity of the auxiliary flow can be described by lubrication function (Snoeijer & Andreotti Reference Snoeijer and Andreotti2013),

(3.3) \begin{align} \boldsymbol{\nabla }\boldsymbol{\cdot }\big ( h^{3}\boldsymbol{\nabla }\hat {p} \big ) = 6\mu \boldsymbol{\hat {U}} \boldsymbol{\cdot }\boldsymbol{\nabla }h. \end{align}

The total hydrodynamic resistance $\hat {\boldsymbol{F}}_d$ in the auxiliary field is dominated by logarithmic singularity of the wetting perimeter. This yields the characteristic scaling, $\hat {{F}}_d = \xi \mu R \hat {U}l_n/\theta$ . Equation (3.2) provides a rigorous connection between droplet velocity and capillary forcing, independent of the detailed internal flow structure, yields

(3.4) \begin{align} \frac {\xi R \mu l_n}{\theta }\boldsymbol{U}= \int _C \gamma \boldsymbol{t} \, \text{d}l, \end{align}

where $l_n$ denotes the logarithmic cutoff factor arising from the viscous stress singularity near the contact line, and $\xi$ is a dimensionless coefficient of order unity reflecting details of the microscopic flow structure. Equation (3.4) can be rewritten as $F_d=F$ , where the functional form of drag force reads $F_d = \xi \mu R Ul_n/\theta$ . The droplet velocity is determined once the effective surface tension in the contact-line region is related to the ambient vapour distribution, which is addressed in the following sections.

3.2.2. Humidity distribution around sessile droplets

When analysing the motion of a target droplet, we take the vapour concentration field imposed by a neighbouring ‘source’ droplet as the relevant driving field. Evaporation of a sessile droplet under ambient conditions is governed by diffusion in the surrounding gas phase (Picknett & Bexon Reference Picknett and Bexon1977; Deegan et al. Reference Deegan, Bakajin, Dupont, Huber, Nagel and Witten1997, Reference Deegan, Bakajin, Dupont, Huber, Nagel and Witten2000; Hu & Larson Reference Hu and Larson2002), so the local vapour concentration $\phi$ satisfies

(3.5) \begin{align} \frac {\partial \phi }{\partial t}=D_v {\nabla} ^2\phi , \end{align}

where $D_v$ is the diffusion coefficient of water vapour in air. In the present water–PG system, evaporation is dominated by water because the saturation vapour pressure of water ( $p_{\textit{sat}}^{\mathit{H_2O}}\approx 2.3\,\mathrm{kPa}$ ) is orders of magnitude larger than that of PG ( $p_{\textit{sat}}^{\textit{PG}}\approx 10\,\mathrm{Pa}$ ) under ambient conditions (Stull Reference Stull2002; Wagner & Pruß Reference Wagner and Pruß2002). Although vapour-phase PG adsorption can affect droplet motion in sufficiently dilute PG solutions (Zhao et al. Reference Zhao, Orejon, Sefiane and Shanahan2025), the concentrations considered here lie outside that regime (figure S6 in the Supplementary material). We therefore neglect PG vapour adsorption and treat the vapour phase as a water–air mixture, with PG taken as effectively non-volatile. For simplicity, we express the vapour field in terms of the relative humidity $\phi (r)$ , where $r$ is the distance from the centre of the ‘source’ droplet. The relevant vapour Péclet number is $ \textit{Pe}_v= \textit{UR}/D_v$ . Using $D_v=2.6\times 10^{-5}\,\mathrm{m^2\,s^{-1}}$ at $T=293\,\mathrm{K}$ and $p=1\,\mathrm{atm}$ (Hu & Larson Reference Hu and Larson2002), we obtain $ \textit{Pe}_v\sim O(10^{-3})$ , indicating diffusion-dominated and quasisteady transport. The humidity field therefore satisfies

(3.6) \begin{align} \boldsymbol{\nabla} ^2\phi =0, \end{align}

which forms the basis of the vapour-field analysis in this study.

The droplets in the experiments are sessile spherical caps of contact angle $\theta$ on a solid substrate. Since the vapour field around an evaporating sessile droplet depends on geometric confinement (Benusiglio et al. Reference Benusiglio, Cira, Lai and Prakash2017), the chamber dimensions are chosen to be at least two orders of magnitude larger than the droplet size, and the droplets are placed sufficiently far from the walls so that confinement effects are negligible. Following the treatment by Cira et al., the droplets considered here are predominantly composed of water (10 wt % PG by volume, corresponding to 2.6 % PG in mole fraction), and the water vapour at the liquid–air interface is therefore assumed to be locally saturated (Raoult Reference Raoult1887; Saxena & Hildemann Reference Saxena and Hildemann1997; Cira et al. Reference Cira, Benusiglio and Prakash2015; Zhao et al. Reference Zhao, Orejon, Sefiane and Shanahan2025). A no-flux condition is imposed on the dry substrate outside the droplet footprint, and the ambient humidity $\phi _\infty$ is recovered in the far field. Under these boundary conditions, Laplace’s equation admits an analytical solution, derived in toroidal coordinates (Picknett & Bexon Reference Picknett and Bexon1977; Hu & Larson Reference Hu and Larson2002),

(3.7) \begin{align} &\frac {\phi -\phi _\infty }{1-\phi _\infty }=\sqrt {2\cosh \alpha -2\cos \beta }\notag \\[5pt]& \qquad \times \int _0^\infty \frac {\cosh {(\theta \omega )}}{\cosh {(\pi \omega )}} \frac {\cosh {[(2\pi -\beta )\omega ]}}{\cosh {[(\pi -\theta )\omega ]}}P_{-1/2+i\pi }(\cosh \alpha )\text{d}\omega , \end{align}

where $\alpha$ and $\beta$ are toroidal coordinates related to the cylindrical coordinates $(r,z)$ , and $P_{-1/2+i\pi }(\cosh \alpha )$ denotes the Legendre function of the first kind.

Of particular relevance for droplet–droplet interactions is the vapour concentration distribution along the substrate outside the droplet footprint. Since the vapour concentration converges within the small experimental contact angles ranging $ 5^\circ \leqslant \theta \leqslant 15^\circ$ , the analytical solution can be reasonably simplified and reduced to the classical problem of diffusion from a circular disk into a half-space. Along the substrate ( $z = 0$ ), the vapour concentration exhibits an inverse-distance dependence on the horizontal distance $r$ ,

(3.8) \begin{align} \frac {\phi - \phi _\infty }{1 - \phi _\infty } = \frac {1}{2} - \frac {1}{\pi } \arctan {\left [ \frac {r^2 - 2R^2}{2R\sqrt {r^2-R^2}} \right ]},\, r \gt R , \\[-28pt] \nonumber \end{align}

(3.9) \begin{align} \frac {\partial \phi }{\partial r}\bigg |_{d} = - (1-\phi _\infty ) \,\frac {2}{\pi d}\, \frac {R}{\sqrt {d^2-R^2}} = -(1-\phi _\infty ) \,\frac {\lambda (d)R}{d^2}, \\[0pt] \nonumber \end{align}

where $d=L+2R \in [2R,+\infty )$ is the centre-to-centre distance of two droplets, and $\lambda (d) = 2/( \pi \sqrt {1-R^2/d^2})$ . Notably, $\lambda (d)$ is a slowly varying, order-unity function of the separation, bounded between its near-field and far-field limits, $\lambda (2R) = {4}/{(\sqrt {3}\pi )}, \lambda (\infty ) = {2}/{\pi }$ . In the following analysis, we take $\lambda =\lambda (2R) \simeq 0.735$ , and the humidity distribution can be expressed as

(3.10) \begin{align} \phi (r) = \phi _\infty + \lambda (1-\phi _\infty )\frac {R}{r}. \end{align}

In this representation, geometric effects are fully absorbed into $\lambda$ , while the humidity field retains a universal functional form (figure 3 a).

To assess the validity of this asymptotic representation, we performed finite-element simulations of Laplace’s equation. The numerical results confirm that, for $\theta \lesssim 30^\circ$ , the vapour concentration profiles along the substrate collapse onto a nearly universal curve that is weakly dependent on the contact angle (figure S1 in the Supplementary material), thereby supporting the validity of (3.10) over the experimentally relevant range.

3.2.3. Dependence of contact angle on humidity

To quantify the effect of ambient humidity on the effective surface tension at the contact line, we consider a sessile droplet in a spatially uniform environment with prescribed relative humidity $\phi$ . In this setting, the droplet remains stationary and the humidity field is homogeneous, the variations in the surface tension are reflected in changes in the apparent contact angle $\theta$ , as given by the horizontal force balance at the intersection of the thin film and the bulk droplet, $\gamma = \gamma _d \cos \theta$ . The distribution of $\phi$ produces a gradient in $\gamma$ , the surface tension of the precursor liquid, and subsequently the driving force for droplet motion (Cira et al. Reference Cira, Benusiglio and Prakash2015). Here, $\theta$ should be the static apparent contact angle measured under spatially uniform humidity conditions, rather than the dynamic angle during droplet motion.

Figure 3.

(a) The local relative humidity $\phi$ along the substrate. Blue curve, theoretical results (3.8), for a thin film with $\theta$ = 0; dots, numerically solved results of $\phi$ , for the ambient humidity in the far field $\phi _{\infty } = 30.0\,\%$ and $\theta = 10^{\circ }$ ; red curve, $\phi (r)=\lambda (1-\phi _\infty )R/r+\phi _\infty$ , an asymptotic solution of (3.8) with $\lambda$ = 0.735; black curve, centrosymmetric solution of $\phi (r)$ for the evaporation from a drop floating freely in the air (Cira et al. Reference Cira, Benusiglio and Prakash2015). (b) Apparent contact angle $\theta$ as a function of ambient relative humidity $\phi$ for a 10 wt % PG droplet on a hydrophilic silicon substrate. The fitting coefficients are $m_l=0.036$ and $b_l=0.96$ , $m_t=0.28$ and $\phi _{\textit{eq}}=0.95$ , for two models, respectively.

Several empirical and theoretical relationships have been proposed to characterize the dependence of the apparent contact angle on the ambient humidity. Reported scalings range from $\theta \sim (\phi _{\textit{eq}}-\phi )^{1/3}$ to $\theta \sim (\phi _{\textit{eq}}-\phi )^{1/2}$ , depending on whether evaporation-induced Marangoni flows are treated explicitly within a lubrication-theory framework (Karpitschka et al. Reference Karpitschka, Liebig and Riegler2017; Wu et al. Reference Wu, Doi and Man2021) or incorporated implicitly through effective interfacial relations (Cira et al. Reference Cira, Benusiglio and Prakash2015). Careful experimental measurements typically yield effective exponents lying between these limiting values, reflecting the combined influence of thermodynamic and hydrodynamic effects near the contact line. Figure 3(b) shows our experimental data together with both scalings discussed. They capture the observed trends over the accessible humidity range, with minor quantitative differences that do not affect the droplet dynamics analysis.

Overall, these $\theta -\phi$ relationships supply surface tension variation combining evaporation, interfacial transport, compositional rearrangements and possible Marangoni-induced flows within the contact-line region, without the need to resolve the exact form of surface tension. They provide a basis for modelling the droplet dynamics addressed in this work.

3.2.4. Inverse-square force and cubic kinematics

Building on the humidity analysis and the experimental calibration of the humidity-dependent contact angle established above, we now quantify how the resulting surface tension heterogeneity generates a net driving force on a sessile droplet and governs its subsequent motion. As shown in (3.4), the capillary driving force acting on the sessile droplets originates from the circumferential variation of $\gamma$ along the contact line. Projected onto the direction of motion, the force reads $F =\boldsymbol{F} \boldsymbol{\cdot }\boldsymbol{e_x}=\int _{0}^{2 \pi } \gamma \cos \varphi R \text{d} \varphi$ , where $\varphi$ is the azimuthal angle, $\boldsymbol{e_x}$ is the unit vector in the x-direction, as illustrated schematically in figure 1(b). Under the quasistatic assumption, the humidity and the associated surface tension of the precursor film $\gamma$ are assumed to remain slaved to the instantaneous droplet configuration. Accordingly, $\gamma$ adjusts instantaneously to its local equilibrium value set by the humidity with $\gamma (\phi )=\gamma _d\cos \theta (\phi )$ . Here $\theta$ parametrizes the macroscopic response of the contact line to humidity, rather than representing a dynamic contact angle during droplet motion. As mentioned above, contact angle is a function of humidity $\cos {\theta }=f(\phi )$ . To simplify the analysis, $\cos {\theta }$ is expanded linearly,

(3.11) \begin{align} f(\phi ) \equiv \cos \theta = f^{\prime }\left (\phi _{\infty }\right )\left (\phi -\phi _{\infty }\right )+f\left (\phi _{\infty }\right )=m \phi +b, \end{align}

with $m=f^{\prime }(\phi _{\infty })$ and $b=-f^{\prime }(\phi _{\infty })\phi _{\infty }+f(\phi _{\infty })$ . Here

(3.12) \begin{align} F=\int _{0}^{2 \pi } \gamma _{d} \cos \theta \cos \varphi R \mathrm{d} \varphi = \gamma _d f'(\phi _\infty ) \int _0^{2\pi } \phi R \cos \varphi \, \text{d}\varphi . \end{align}

Performing a local Taylor expansion of $\phi$ about the ‘target’ droplet centre gives

(3.13) \begin{align} \int _0^{2\pi }\phi R \cos \varphi \, \text{d}\varphi =\pi R^2 \frac {\partial \phi }{\partial d}, \end{align}

where the zeroth-order contribution vanishes by angular symmetry. This expression shows that the driving force may be interpreted as the product of the wet area of the droplet and the vapour-induced surface-tension gradient evaluated along the direction of motion, with the latter given, to leading order in the humidity perturbation, by the directional derivative of the humidity field at the droplet centre. Therefore, the distribution of the humidity field contributes significantly to the functional form of the driving force. The resulting driving force reads

(3.14) \begin{align} F=-\frac {\gamma _{d} m \lambda \pi R^{3}\left (1-\phi _{\infty }\right )}{d^{2}}. \end{align}

The higher-order terms contributing a relative correction can be neglected, which is discussed in the Supplementary material (§ II).

Figure 4.

Inclined-plane sliding experiments used to assess the value of $l_n$ by a moving droplet. (a) Schematic of a droplet on an inclined silicon wafer (15 $^{\circ }$ ). (b) Generation of a 10 wt % PG aqueous droplet of known volume $\varOmega$ from a conical capillary tube (scale bar is 1 mm), ensuring reproducible gravitational force along the incline. (c) Droplet sliding on the inclined wafer at its terminal velocity (scale bar is 2 mm). (d) Comparison of the measured force–velocity relation with the linear scaling and the Cox–Voinov asymptotic form. The drag is evaluated against the component of gravity acting parallel to the substrate.

To close the dynamical problem, it is necessary to determine the coefficient appearing in (3.4). According to lubrication theory for droplets with small contact angles and low capillary number ( $Ca \ll 1$ ), viscous dissipation is dominated by flow near the moving contact line, yielding a linear force–velocity relation $F_d = \xi \, \mu R U l_n/{\theta }$ , where the dimensionless prefactor $\xi = 3\pi$ arises from the radial and circumferential integration of shear dissipation in the wedge-shaped lubrication flow (de Gennes Reference de Gennes1985). Dedicated sliding experiments were then performed to quantify the logarithmic cutoff factor $l_n$ associated with the motion of the same binary droplet, thereby closing the model. For comparison, we consider the Cox–Voinov asymptotic scaling (Voinov Reference Voinov1976; Cox Reference Cox1986; de Gennes et al. Reference de Gennes, Brochard-Wyart and Quéré2003; Reyssat Reference Reyssat2014), which predicts $F_d \sim \gamma R \, Ca^{2/3}$ , reflecting the contribution of dynamic contact-angle effects to the viscous resistance (see the Supplementary material § IV).

Figure 5.

(a) Temporal evolution of the separation distance $L$ as a function of $\tau -t$ for different droplet radii $R$ and ambient relative humidities $\phi _\infty$ . (b) The same data plotted in the rescaled representation collapse onto straight lines. The experiments span droplet radii in the range 0.8–1.8 mm and relative humidities from 19.9 % to 81.1 %. (c) Comparison between the prefactors $k_t$ calculated theoretically from (3.15), and $k_e$ measured experimentally from (b). Data are selected for $\phi _{\infty } \sim$ 60 %. (d) The ratio of $k_e/k_t$ versus environment humidity $\phi _{\infty }$ . Fitting coefficients are $m_l = 0.036$ and $b_l = 0.96$ , $m_t = 0.28$ and $\phi _{\textit{eq}} = 0.95$ .

Figure 5. Long description

Panel A: A line graph shows the temporal evolution of the separation distance as a function of time for different droplet radii and ambient relative humidities. The x-axis represents time in seconds, and the y-axis represents the separation distance in millimeters. Different colors and symbols represent various droplet radii and relative humidities. Panel B: The same data is plotted in a rescaled representation, collapsing onto straight lines. The x-axis represents the rescaled time, and the y-axis represents the rescaled separation distance. Panel C: A scatter plot compares the prefactors calculated theoretically from a specific equation and measured experimentally. The x-axis represents the theoretical prefactor, and the y-axis represents the experimental prefactor. Data points are selected for a specific relative humidity. Panel D: A scatter plot shows the ratio of coefficients versus environment humidity. The x-axis represents the environment humidity, and the y-axis represents the ratio of coefficients. Fitting coefficients are provided in the caption.

As shown in figure 4(d), the measured force–velocity relation follows a linear scaling, with no distinct $Ca^{2/3}$ trend observed over the explored range. It may be attributed to Marangoni-induced contraction in binary droplets, as previously discussed in similar systems in Karpitschka et al. (Reference Karpitschka, Liebig and Riegler2017) and Shiri et al. (Reference Shiri, Sinha, Baumgartner and Cira2021).

Combining the driving force (3.14) with the measured drag $F_d$ into (3.4) yields the droplet velocity

(3.15) \begin{align} U = - k/(6d^2)= - k/[6(L+2R)^2], \end{align}

where $ k = 2 \gamma _d m \lambda R^2 \theta (1 - \phi _\infty ) / \mu l_n$ . Integrating $2U = \text{d}L/\text{d}t$ yields a cubic kinematic scaling for the centre-to-centre distance $d=L+2R$ ,

(3.16) \begin{align} d^3 - (2R)^3 = k (\tau - t), \end{align}

where $\tau$ denotes the time of contact.

Applying this transformation to the experimental data reveals a striking linear relation (figures 5 a and 5 b) that all curves collapse onto lines passing through the origin. This confirms that the observed droplet dynamics are consistent with the inverse-square force scaling of (3.14). These experiments covered droplet radii ranging from 0.8 to 1.8 mm and $\phi$ from 19.9 % to 81.1 %. Then, we compare the theoretical prefactor $k_t$ computed from (3.15) with the experimentally extracted value $k_e$ , as shown in figures 5(c) and 5(d). Good agreement is observed for $\phi _\infty \sim 60\,\%$ across both calibrations of $f(\phi )$ . For more extreme experimental conditions with $\phi _\infty \lt$ 25 % or $\gt$ 75 %, $k_e / k_t$ calculated from both models deviates from 1, but still spreads mainly between 0.5 and 1.5. The deviation possibly originates from the use of $\theta = f(\phi _\infty )$ to predict the influence of local relative humidity on $\theta$ (figure 3 b). A more precise model of $f(\phi )$ may be worth pursuing, nevertheless, it is shown from figure 5(b) that (3.16) holds for all tested experimental conditions with $\phi _\infty \in [19.1\,\%, 81.1\,\%]$ , indicating the inverse-square nature of $F$ is beyond reasonable doubt.

Figure 6.

Comparison of the droplet kinematics under different experimental conditions. Here (a) 10 wt % PG droplet pairs with different initial radii. (b) Droplet pairs with identical PG concentration in each experiment, where the concentration is varied across different values. (c) Comparison between the theoretical prediction and experimental data reported by Cira et al. (Reference Cira, Benusiglio and Prakash2015), where the theoretical curves are evaluated using the parameters extracted from their measurements. The slope of the experimental fit, $k_e = 15.89$ , is consistent with the theoretical value $k_t = 16.22$ .

Figure 6. Long description

Panel A: A line graph shows the relationship between the cube of the difference in radii and time for droplet pairs with different initial radii. The x-axis represents time in seconds, and the y-axis represents the cube of the difference in radii in cubic millimeters. Three datasets are plotted: red circles for a radius ratio of 1.21 and a rate constant ratio of 1.13, black squares for a radius ratio of 2.02 and a rate constant ratio of 1.06, and blue diamonds for a radius ratio of 4.27 and a rate constant ratio of 1.70. The lines indicate theoretical predictions, with the red line showing a steeper slope than the black and blue lines. Panel B: Another line graph depicts the relationship between the cube of twice the radius and time for droplet pairs with identical concentrations but varying radii. The x-axis represents time in seconds, and the y-axis represents the cube of twice the radius in cubic millimeters. Three datasets are shown: red circles for 20 weight percent propylene glycol with a radius of 1.04 millimeters, black squares for 40 weight percent with a radius of 1.23 millimeters, and blue diamonds for 60 weight percent with a radius of 1.54 millimeters. The lines represent theoretical predictions, with the red line having the steepest slope, followed by the black and blue lines. Panel C: A line graph compares theoretical predictions with experimental data from Cira et al. The x-axis represents time in seconds, and the y-axis represents the cube of twice the radius in cubic millimeters. Red circles represent experimental data, the dashed line shows a linear fit to the experimental data, and the solid blue line represents the theoretical prediction. The experimental data and the theoretical prediction show a similar trend, with the linear fit closely matching the theoretical prediction.

The effect of transient droplet deformation near the final stage of approach ( $\tau -t\sim 0$ ) is shown to be negligible within the present parameter range (see the Supplementary material § V for details). Accordingly, the droplets may be treated as maintaining a fixed shape in the dynamical analysis. In addition, the maximum surface temperature difference induced by evaporative cooling is of order $10^{-3}\,\mathrm{K}$ (Hu & Larson Reference Hu and Larson2002). The associated thermally induced Marangoni force ( $\sim 10^{-9}\,\mathrm{N}$ ) is at least two orders of magnitude smaller than the force arising from concentration-induced surface-tension gradients ( $\sim 10^{-7}\,\mathrm{N}$ ). Evaporative thermal effects can therefore be safely neglected in the present dynamical analysis.