Dynamics of evaporating, interconnected droplets

Abstract

We report on the dynamics of a pair of sessile droplets that are connected by a microchannel, yet open to the atmosphere and hence free to evaporate. Our results reveal that fluid exchange between droplets occurs via a pumping flow driven by differences in hydrostatic and Laplace pressure between the two droplets. Evaporation causes the droplets to slowly lose volume and change shape, which subsequently affects the fluid transport between them. We observe that, for equal contact areas, a larger droplet typically feeds a smaller droplet during evaporation and the flow in the connecting channel is unidirectional. However, for unequal contact areas, the flow can reverse in the connecting channel following a sudden switch in droplet shape that occurs during evaporation. A stability analysis reveals that the dynamics of the exchange flow are underpinned by a supercritical pitchfork bifurcation. Evaporative volume loss permits the droplet pair to step through a sequence of quasi-stationary states determined by the instantaneous volume of the system. Enforcing unequal contact area unfolds the bifurcation such that droplet-shape switching and the associated flow reversal can be understood in terms of a jump from the disconnected to the connected branch of the bifurcation. This establishes symmetry breaking as a mechanism to induce evaporation-driven flow reversal in connected droplets.

---

1 Introduction

Microfluidic devices that are passively driven, rather than requiring external pumping, have great potential as point-of-care diagnostic tools,¹ due to their compactness and given that they do not require tethering to hardware, e.g., a pressure controller or syringe pump. Capillary action is one mechanism that provides a means of passively transporting fluids through stand-alone devices and hence rendering them portable. Capillary pumping, based on gradients in Laplace pressure, can be achieved in closed systems with ‘fluid walls’,^2,3 by confining liquid circuitry beneath an immiscible fluid layer to prevent evaporation of the working fluid. However, open microfluidics have the advantage of being accessible and straight-forward to manufacture^4–6 and have found application as diagnostic tools and in biochemical synthesis.^4,7

Droplet driven capillary micropumps^8–12 rely on flows between droplets of different sizes i.e., a ‘pumping droplet’ and a ‘reservoir droplet’. Liquid transport is primarily driven by internal pressure differences between droplets and flow control can be engineered by considering the material properties (hydrophobic vs. hydrophilic) and geometry of the device and the droplets. Such microfluidic platforms can reliably deliver microlitre volumes of fluid without the need for external energy input,¹³ making them particularly suitable for applications requiring precise handling of small liquid volumes in resource-limited or portable settings.

On short time scales, droplet evaporation can typically be neglected. However, for volatile liquids, evaporation is inevitable. Governed by the diffusion of liquid molecules into the surrounding gas phase, evaporation typically occurs over long timescales and can therefore be harnessed to induce sustained microfluidic flows which can persist for minutes to hours.¹⁴ This has been demonstrated in sweat-sensing devices, which imbibe sweat through a filter–paper interface and sustain continuous fluid flow without any external pump.^15,16 In open capillary micropumps, however, evaporation leads to a reduction in droplet volume that alters droplet shape, leading to changes in the pressure difference between the droplets, which in turn affects the pumping flow. Yet the influence of evaporation on open droplet-based microfluidic flows – including the potential loss of working fluid – is typically overlooked.

Evaporative effects have been considered in other systems comprising multiple, sessile droplets that are in close proximity.^17–21 The vapour emitted by one droplet increases the local ambient vapour concentration around its neighbours, thereby reducing the local vapour concentration gradient that drives evaporation. As a result, droplets in close proximity exhibit reduced evaporation rates and prolonged lifetimes compared to their isolated counterparts – a phenomenon known as ‘shielding’. Here, we minimise vapour-mediated interactions between droplets, by separating the droplets from one another, and instead focus on how fluid is exchanged between droplets via a connecting channel.

We investigate a model capillary micropump based on two connected droplets and simultaneously examine both the fluid transport between two droplets and the volume loss due to evaporation. Our results reveal that, under certain conditions, the direction of capillary flow between droplets can reverse, rendering ephemeral the distinction between the ‘pumping’ and ‘reservoir’ droplets. Such a flow reversal has been previously observed in droplet-based micropumps,^22,23 but was attributed to inertial effects. We model the dynamics of interconnected droplets in the absence of inertia and, in doing so, demonstrate that evaporation-induced flow reversal can instead be a direct consequence of device and droplet geometry.

2 Experimental method

In experiments, a pair of droplets was positioned upon two pillars connected via a channel. As the droplets evaporated, the evolving profile of each droplet was monitored and the flow field in the connecting channel was measured. A schematic diagram of the experimental set up is shown in Fig. 1. Experiments were performed using droplets of deionized water, containing 0.04% (w/v) polystyrene particles of diameter 5 μm (TSI Inc. KOBO, SP-500, refractive index 1.53, density ρ = 1.01 g mL⁻¹), in ambient laboratory environments for which the temperature and humidity were measured to be 23 ± 0.1 °C and 47 ± 0.5%, respectively. The surface tension of the working fluid was measured to be γ = 71.5 ± 0.5 × 10⁻³ N m⁻¹ using a Goniometer (Attension Theta Flex, Biolin Scientific), while we use literature values for viscosity μ = 0.95 mPas and density ρ = 1000 kg m⁻³. We terminate our experiments once we are no longer able to resolve droplets of finite volume or a droplet contact line unpins from the edge of the pedestal on which they rest.

Fig. 1 Experimental system. Schematic diagram of the experiment, comprising two sessile droplets resting on pillars of circular cross-section with radii a₁ and a₂, respectively. Each pillar has a circular hole at its center, denoted by O₁ and O₂, that forms the opening for a vertical, cylindrical channel. The two cylindrical channels are connected by a horizontal microchannel of rectangular cross-section. The contact angles of the droplets are represented by θ₁ and θ₂, while their heights are indicated by H₁ and H₂. Since the droplets resemble spherical caps, their radii of curvature, r₁ and r₂, were determined by fitting circles to the droplet profiles. Consequently, for a sessile droplet with a contact angle less than π/2, the center of the circular profile is depicted below the droplet base, whereas for a droplet with a contact angle above π/2, the center is located above the base.

The pillars and connecting channel were fabricated from Perspex using a micromiller. Prior to experimentation, the fabricated device was soaked in deionized water with 10% of a surface active cleaning agent (Decon 90) for 10 hours, and then rinsed in deionized water and dried with compressed nitrogen. This cleaning procedure prevented surface contamination – a result of the micromilling process – from affecting the surface tension of the droplets.²⁴ Before droplets were deposited on the pillars, the connecting channel was filled manually using a syringe. The droplets were then carefully deposited upon their respective pillars using a micropipette (Acura 826 XS, Socorex) with range 0.5–10 μL and accuracy 0.02–1.8%. The contact line of each droplet remained pinned at the pillar edge throughout the experiment,²⁵ such that droplet contact areas are constant. The diameter of pillars and initial droplet volumes were varied in experiments but the dimensions of the connecting channel remained the same; a horizontal section of length 10 mm and cross-section of 0.5 mm × 0.5 mm was connected at its ends to the base of each pillar via vertical cylindrical channels (see SI for exact channel geometry). The length of the connecting channel maintained the two droplets sufficiently far away from one another (i.e., greater than 10 droplet radii apart) that collective evaporation effects are considered to be small; the strength of vapour-mediated interactions is predicted to scale with 1/r¹⁹ in quiescent environments however in the laboratory, air currents produced by e.g., air-conditioning, likely diminish the interaction range considerably. Droplets were illuminated with uniform, white back-lighting and side-view profiles were imaged using a CCD camera (Genie, Teledyne Dalsa) with a spatial resolution of 100 pixels per mm and a frame rate of 40 images per second.

Throughout experiments, the Bond number Bo = ρgH²/γ, where g is the acceleration of gravity and H is the droplet height, is less than 0.86, which implies that hydrostatic contributions to the pumping flow between droplets were small yet non-negligible relative to Laplace pressure contributions. However, droplet shape remained well-approximated by a spherical cap, so we acquired measurements of the radius, r(t), height, H(t) and volume V(t) of each droplet during evaporation by fitting circles to droplet profiles and assuming axisymmetry; fits and measured profiles agreed to within 2.4% throughout.

The pressure difference between channel ends were so small, however, that they were beyond the lower resolution limit of commercial pressure sensors. Hence, rather than measuring pressure differences, we assessed the pumping flow in the connecting channel using flow visualization. To measure flow fields and flow velocities, the motion of tracer particles in the interconnecting channel was imaged using a microscope (AD4113ZTL, Dino-Lite) and we performed particle image velocimetry (PIV) on images acquired with a time interval of 1 s, using the PIVlab toolkit in Matlab. To measure the in-plane flow field in the connecting channel, the camera was placed underneath the channel, focused on the centre-line of the channel, and imaged from below with 130.5× magnification and a spatial resolution of 1.10 μm per pixel. For modelling purposes, we are primarily concerned with the flow rate inside the connecting channel, where u is the axial velocity and S the cross-sectional area. We approximate Q ≈ ūS where ū is the depth-averaged velocity that, for a channel of square cross-section, is related to the maximum speed along the center line by u_max(z = 0, y = 0) = 2.115ū.²⁶ To visualize the flow inside the droplets, the mid-plane of the droplets was illuminated using a light sheet of width 200 μm. To minimize the influence of heat on the evaporation rate of the droplet, we generated the light sheet using a cold light source (KL 2500 LED, Schott) in combination with a cylindrical lens.

3 Theory

The pressure at the channel ends, points O_i in Fig. 1, comprises hydrostatic and Laplace pressure contributions,where P_atm is atmospheric pressure, γ is the surface tension of water, H is the height of droplet and r the radius of curvature of the free surface of a drop, g is gravitational acceleration and i = 1,2 denotes the droplet number. (We take the convention that i = 1 represents the droplet on the left-hand side of side-view images and i = 2 represents the droplet on the right).

Fluid is exchanged between interconnected droplets via a pumping flow in the connecting microchannel. We assume the liquid to be both incompressible and Newtonian, and model the flow in the microchannel as fully-developed laminar flow with no-slip boundary conditions, such that the flow rate Q is related to the pressure drop between the two droplets by the Hagen–Poiseuille equation:

Q = K⁻¹ΔP, (2)

where ΔP = P₁ − P₂ is the pressure difference at the two ends of the connecting channel and K = 8.73 × 10⁹ kg m⁻⁴ s⁻¹ is the resistance coefficient of the connecting channel (see SI for further details).

Droplet evaporation is controlled by the diffusion of vapour from the free surface to the ambient environment. The rate of volume loss of an individual sessile droplet with constant contact radius is given by²¹

where

The volume of each droplet, V_evap,i, evolves in time because of evaporation, which is driven by a pressure difference Δp = p₀ − p_∞ between the saturation vapour pressure of water, p₀ = 2.065 kPa, and ambient vapour pressure far away from the surface of the droplets p_∞; here, we determined the ambient vapour pressure from the measured relative humidity RH = p_∞/p₀ = 0.47, giving Δp = p₀(1 − RH) = 1.094 kPa. Since the contact radius of the droplets, a_i, is fixed by the pillar geometry, the apparent contact angle of the droplets, θ_i, evolves in time. Here, D = 2.3 × 10⁻⁵ m² s⁻¹ is the diffusion coefficient of water vapour molecules in air, M = 18 g mol⁻¹ is the molar mass of water, and R = 8.3 J mol⁻¹ K⁻¹ and T are the gas constant and ambient temperature, respectively, and f(θ) is an empirical function^27,28 given by

f(θ) = 0.00008957 + 0.6333θ + 0.116θ² − 0.08878θ³ + 0.01033θ⁴, (5)

for 0.175 ≤ θ ≤ π, which matches the range of angles we can resolve experimentally.

For a pair of interconnected droplets, the volume of each individual droplet is therefore influenced both by evaporation and the pumping flow within the connecting channel. As such, the evolution of individual droplet volumes is coupled, and can be written in terms of volumetric change associated with both evaporation and pumping flow,

Assuming that the droplets take the form of a spherical cap we invoke the following geometrical relations between droplet volume V_i, height H_i, radius r_i and contact angle θ_i:

This permits us to rewrite eqn (6) in terms of droplet contact angles, θ_i, as

where

Since the evaporative terms are dissipative, eqn (10) has one stationary solution, θ₁ = θ₂ = 0, which physically implies that evaporation can drive a pair of interconnected droplets to evolve to have net zero volume, V₁ = V₂= 0, irrespective of initial conditions. The evolution of θ_i with time can be calculated numerically with the built-in initial value problem solver ODE45 in Matlab, for initial values of θ_i,t = 0.

4 Results

4.1 Experimental observations

Our experiments reveal that the evolution of an interconnected droplet pair occurs on two timescales, see Fig. 2, as fluid is exchanged through the connecting channel and evaporation occurs. Within t ≲ 1 s of being deposited on the pillars, a rapid redistribution of fluid between droplets causes sudden changes in the size of individuals droplets. This early-time behaviour demonstrates how the pumping flow – driven by a pressure difference between the two ends of the channel, ΔP, and moderated by channel resistance, as described by eqn (1) and (2), respectively – exchanges fluid between droplets in the absence of evaporation. In experiments, the time taken for this initial fast exchange flow to stop varies because of differences in initial droplet volumes (rather than differences in channel resistance). On a longer timescale, t ≫ 1 s, droplet volumes change slowly as evaporation decreases the total volume of the system. However, a pumping flow may persist during the slow evaporation stage, and act to replenish fluid lost through evaporation such that a single droplet may maintain its size as the volume of the other droplet decreases.

Fig. 2 Evolution of a pair of connected evaporating droplets. Images acquired at various instance in time (as indicated) demonstrate a rapid redistribution of fluid between droplets at early times (t ≲ 1 s) followed by a slow evolution of droplet shape at later times (t ≫ 1 s). (A) Two droplets of low initial volume on connected pillars of equal radii adopt a wetting configuration, such that both droplets have initial contact angles θ₁, θ₂ < π/2. They rapidly acquire the same volume (t ≲ 25 ms) and then remain approximately identical in size and shape as they evaporate (see Video S1 in SI). (B) Two droplets of markedly-different initial volume on connected pillars of equal radii adopt an asymmetric configuration with θ₁ > π/2 and θ₂ < π/2. As they evaporate, the volume of the larger droplet decreases while the volume of the smaller droplet increases until the two have the same volume, which indicates that there is a net flow in the connecting channel from the larger drop to the smaller drop (see Video S2 in SI). (C) Droplets on unequal pillars; larger droplet on the larger pillar. After the rapid transition period (t ≲ 500 ms), θ₁ > π/2 and θ₂ < π/2. With evaporation, the volume of the smaller droplet remains nearly constant, indicating that there is a net flow in the connecting channel from the larger droplet to the smaller droplet (see Video S3 in SI). (D) Droplets on unequal pillars; smaller droplet on the larger pillar. During the slow evaporation stage the droplet volumes switch suddenly; at t ∼ 400 s the position of the largest drop rapidly changes from pillar (2 → 1) (see Video S4 in SI).

First, we consider the case of droplets on connecting pillars of equal radii i.e., for a₁ = a₂= 0.5 mm, see Fig. 2(A) and (B). When both droplets have initial contact angles θ₁, θ₂ < π/2, they rapidly acquire the same volume, as shown in Fig. 2(A) for which V_1,t = 0 = 0.15 μL and V_2,t = 0 = 0.25 μL. Initially, within t ≲ 25 ms of deposition, both the hydrostatic and Laplace contributions to pressure combine to drive fluid from the larger droplet to the smaller droplet, see inset to Fig. 3(A). As they then slowly evaporate on a longer timescale, t ≫ 1 s, the two droplets remain approximately identical in size and shape, see Fig. 3(A). As such, the pressure differential across the channel is null, and the evaporative flux of each droplet is equivalent, so the mean flow velocity within the connecting channel is zero, see Fig. 3(B).

Fig. 3 Connected droplets on equal pillars. Evolution of droplet volumes (A), (C) and channel flow speed (B), (D) measured as a function of time t for pillars of equal radii a₁ = a₂ = 0.5 mm. A and B: V_1,t = 0 = 0.15 μL and V_2,t = 0 = 0.25 μL; the droplets are comparable in size during evaporation and the pumping flow between the droplets is measured to be approximately zero throughout. (C) and (D) V_1,t = 0 = 1.45 μL and V_2,t = 0 = 0.8 μL; the difference in droplet size generates a net pumping flow from drop 1 to drop 2 (i.e., from the larger drop to the smaller drop) during evaporation. Experimental data, from Fig. 2(A) and (B), represented by markers (see legend), numerical data represented by black lines. Insets: Zooming in on early times (t ≲ 1) highlights the rapid redistribution of fluid between two droplets that occurs immediately after they are connected.

If instead, the contact angle of at least one droplet is θ_i > π/2, then at early times the smaller droplet shrinks and the larger droplet grows so that an asymmetric configuration is rapidly established with notable differences in droplet volume and contact angle between the drops. For V_1,t = 0 = 1.45 μL and V_2,t = 0 = 0.8 μL this rapid redistribution of fluid occurs for t ≲ 100 ms, as shown in Fig. 2(B), and suggests that the pumping flow is initially dominated by a Laplace pressure differential. However, as the drops evaporate, the larger droplet shrinks while the smaller droplet grows slightly until the droplets match in size and evaporate out in unison, see Fig. 3(C). Although a flow driven by differences in evaporative flux might be expected to wick fluid towards the larger droplet, which has greater surface area and thus larger evaporative flux, our measured flow velocities in the channel instead provide evidence of a unidirectional flow from the larger droplet to the smaller droplet, see Fig. 3(D), implying that the pumping flow is instead governed by a pressure differential based on capillary statics i.e., Laplace and hydrostatic contributions. Our measurements show that the flow velocity diverges as the droplet volumes become approximately equal, after which it reduces to zero. This implies that pressure differences between the drops continue to drive the flow in the connecting channel as evaporation causes the droplets to change shape.

The distinction between these two aforementioned cases arises from the nonlinear pressure-contact angle relationship of a single sessile droplet, given by eqn (11), which has a local maximum for θ ≈ π/2 that corresponds to a hemispherical drop i.e., for V_i ≈ 2πa_i³/3. In Fig. 4, we show calculated pressures as a function of contact angle for individual droplets with equal and unequal contact radii. The pressure is the sum of Laplace pressure and hydrostatic pressure contributions, see eqn (11), and exhibits a local maximum for θ ≈ π/2 with an amplitude that depends on contact radius. (We highlight that the contact angle of maximum pressure is slightly greater than π/2 because of the additional contribution of hydrostatic pressure.) This nonlinear relation between pressure and contact angle is geometric in origin and can render a droplet bistable for a given pressure, although the corresponding droplet shapes have different volumes.

Fig. 4 Droplet pressure. Numerical values of the pressure calculated from eqn (11) as a function of contact angle for individual sessile droplets of (A) equal contact radii, a₁ = a₂ = 0.5 mm, and (B) unequal contact radii, a₁ = 0.75 mm a₂ = 0.5 mm. In all cases, the relation is non-monotonic with a local maximum for θ ≈ π/2. (In the absence of gravity, the maximum pressure occurs for θ = π/2.)

A pair of droplets, each with maximal pressure, therefore have a total volume equal to that of a sphere with an equatorial cross-section given by the contact area (i.e., pillar base) and this serves as a threshold value, V^crit_tot ≈ 4πa³/3 where a = a₁ = a₂. For a pair of connected droplets with equal contact area and θ₁, θ₂ < π/2, such that V_tot ≲ V^crit_tot, instantaneous pressure equalization can only be achieved with equal volume drops, since the pressure–volume relation of a single drop is monotonic for θ < π/2. Instead, for V_tot ≳ V^crit_tot instantaneous equilibrium can be attained with either an unequal pair (e.g., θ₁ < π/2 and θ₂ > π/2), since contact angle is a double-valued function of pressure in the range 0 < θ <π, or with an equal pair for which θ₁ = θ₂ < π/2. In Section 4.3, we demonstrate that the latter case is an unstable configuration, however, such that perturbations will drive the pair to instead take on different volumes and contact angles.

The interplay of droplet geometry and pressure contributions can be further disentangled by examining interconnected droplets on pillars of unequal radii i.e., for a₁ > a₂, see Fig. 2(C) and (D). First we consider the case in which both droplets have sufficient volume that their initial contact angles θ₁, θ₂ > π/2, and the larger droplet is positioned on the larger pillar – as shown in Fig. 2(C) for which V_1,t = 0 = 2 μL, V_2,t = 0 = 1.2 μL, a₁ = 0.75 mm and a₂ = 0.5 mm. Initially, fluid is rapidly redistributed from the smaller to the larger drop (see inset to Fig. 5(A)), suggesting that the Laplace pressure differential dominates the pumping flow at early times. After this rapid transition, an instantaneous pressure balance is reached with θ₁ > π/2 and θ₂ < π/2. As evaporation occurs, and the system evolves, fluid flows from the larger to the smaller drop such that the larger droplet shrinks whilst the smaller droplet maintains its size, see Fig. 5(A) and (B).

Fig. 5 Connected droplets on unequal pillars. Evolution of droplet volumes (A) and (C) and channel flow speed (B) and (D) measured as a function of time t for pillars of radii a₁ = 0.75 mm and a₂ = 0.5 mm, respectively. A and B: V_1,t = 0 = 2 μL and V_2,t = 0 = 1.2 μL; during evaporation, the difference in droplet size generates a net pumping flow from drop 1 to drop 2 (i.e., from the larger drop to the smaller drop). (C) and (D) V_1,t = 0 = 0.8 μL and V_2,t = 0 = 2.3 μL; initially the larger drop is positioned on the smaller pillar and the pumping flow is from drop 2 to drop 1 but, for t ∼ 400 s, this configuration inverts suddenly and flow reversal is observed in the connecting channel. Experimental data, Fig. 2(C) and (D), represented by markers (see legend), numerical data represented by black lines. Insets: Zooming in on early times (t ≲ 1) highlights the rapid redistribution of fluid between two droplets that occurs immediately after they are connected.

In all of the cases presented so far, the pumping flow established during evaporation has been unidirectional (from the larger to the smaller droplet) or null (for equal volume droplets on equal radii pillars). However, if the smaller droplet, θ₁ < π/2, is initially positioned on the larger pillar and the larger droplet, θ₂ > π/2, on the smaller pillar, we instead see flow reversal during the slow evaporation stage that is associated with a sudden switching of droplet shapes, see Fig. 2(D) and 5(C), (D). Following droplet-shape switching, the pumping flow in the connecting channel continues to feed the smaller droplet as they evaporate out. In Fig. 6, we show measured flow fields inside the two droplets and the connecting channel, and velocity profiles within the connecting channel, at various instants during droplet evaporation, for the same initial conditions as the experiment shown in Fig. 5(C) and (D). Flow visualization demonstrates how droplet-shape switching and flow reversal combine to maintain the pumping flow from larger droplet to smaller droplet. Visualization of flow fields inside the droplets show that before droplet-shape switching, for t ≲ 380 s, the smaller droplet (droplet 1 in Fig. 6(A)) is fed by the larger droplet (droplet 2 in Fig. 6(A)) and likewise after switching, for t ≳ 380 s, the larger droplet continues to feeds the smaller droplet. Flow velocities in the channel were acquired by spatial averaging in the x-direction of instantaneous flow fields, exhibit the parabolic profile typical of fully-developed Hagen–Poiseuille flow, and show how the flow in the connecting channel reverses direction for t ≳ 380 s, see Fig. 6(C). Using measured internal flow velocities, we estimate the scale of the pressure difference induced by the internal flow within the droplets to be ΔP_f ∼ μ(Δv/Δy) ∼ 10⁻⁶ Pa, where μ is the dynamic viscosity of fluids, and Δv is the difference in velocity in the y-direction, between the droplet apex and the channel outlet, such that the distance over which the pressure difference occurs scales with droplet height Δy ∼ H_i. Since this flow-induced contribution to the pressure is much smaller than either the Laplace or hydrostatic contributions i.e., on the order of compared to to , respectively, we neglect the pressure contributions due to internal droplet flows in our analysis. (Although we note that the inertia of this internal flow was responsible for flow reversal in Ju et al.²² and Ju et al.²³).

Fig. 6 Visualising flow reversal. (A) Flow fields, measured on the mid-plane, of two connected droplets at various instances in time following deposition, as indicated; the colorbar represents the magnitude of the velocity in the y-direction. The droplets’ initial total volume V_{tot,t = 0} = 3 μL and, as they evaporate, the droplets switch shape and the flows inside the droplets and in the connecting channel reverse direction. Droplet-shape switching and flow reversal occurs at t ∼ 400 s. (B) Instantaneous flow fields measured in the connecting channel before and after droplet switching demonstrate flow reversal; the colorbar represents the magnitude of the velocity in the x-direction. (C) Measured velocity profiles within the connecting channel at various instances in time before and after droplet-shape switching, as indicated. Data in B and C is from the same experiments as that shown in Fig. 5(C) and (D). Data in A is from a different experiment but with the same initial conditions.

4.2 Droplet dynamics

Numerical solutions to the contact angle evolution equation, eqn (10), enable us to construct a pseudo-phase portrait of an interconnected droplet pair; initial values of θ_i were chosen for portrait clarity. In Fig. 7(A) and (B), we present dynamical portraits of an interconnected droplet-pair in the phase space of the respective contact angles of the two droplets, (θ₁,θ₂) for pillars of equal and unequal radii, and for initial conditions θ_i,t = 0 < π. (We note that for finite contact radii droplets, the θ₁ = θ₂ = π solutions are degenerate, corresponding to infinite volume droplets with zero Laplace pressure.) Trajectories represent the time-evolution of an inter-connected droplet pair; as evaporation progresses and the total volume decreases, the contact angles of each droplet evolve until reaching the stationary point where θ₁ = θ₂ = 0. To highlight the two distinct timescales of the system, we plot solutions for t < 1 s as faint grey lines and solutions for t > 1 s as solid black lines, and include experimental data that is colour-coded by time after droplet deposition.

Fig. 7 Trajectories in pseudo-phase space. Phase portrait of the evolutionary trajectories of an evaporating, interconnected droplet pair, in the pseudo-phase space of the their respective contact angles, for (A) equal pillars of radii a₁ = a₂ = 0.5 mm, and (B) unequal pillars of radii a₁ = 0.75 mm and a₂ = 0.5 mm. Numerical data of the evolution of the contact angles of the two droplets is represented by lines (faint, grey lines represent early-time behaviour for t < 1 s and black lines represent late-time behaviour for t > 1 s), with arrows indicating the direction of droplet evolution. The red star (at θ₁ = θ₂ = 0) indicates the stationary solution of eqn (12), which is the terminal point of the evaporating droplets. Experimental data (for the same parameters as indicated in Fig. 3(C) and 4(C), respectively) is represented by markers where the colour indicates time t, as indicated in the legend colormap. The red arrows indicate the direction of the droplets evolution with time in the experiments. The black arrowheads are positioned at t = 10 ms along trajectories.

For pillars of equal radii, a₁ = a₂ = 0.5 mm, we see that droplet states rapidly converge onto slow evaporation pathways, see Fig. 7(A). If both droplets are initially in a wetting configuration i.e., for θ_1,t = 0 < π/2 and θ_2,t = 0 < π/2, the system quickly converges onto a slow evaporation pathway defined by θ₁ = θ₂ < π/2, and remain on that path as the total volume of the system decreases until θ₁ = θ₂ = 0. As the droplets have equal contact angles and equal contact areas, there exists zero pressure difference between channels ends in this scenario, and hence zero pumping flow along this path (droplets evaporate out together in unison). However, if either of the droplets has an initial contact angle π/2 < θ_i,t = 0 < π, and V_{tot,t = 0} ≳ 4πa_i³/3, then the system rapidly converges onto an evaporation path with one increasing contact angle and one decreasing contact angle, for which a pumping flow feeds the smaller droplet, before eventually reaching the low-volume evaporation pathway defined by θ₁ = θ₂ < π/2. Our numerical results also indicate trajectories rapidly diverge from the initial conditions π/2 < θ₁ = θ₂ < π which suggests this an unstable state; the stability of the system is discussed further in Section 4.3.

In the case of unequal pillars, with a₁ = 0.75 mm and a₂ = 0.5 mm, the symmetry of the system is broken, and this is reflected in the resulting pseudo-phase portrait, see Fig. 7(B). We find that initial conditions rapidly converge onto one of two evaporative pathways. On the θ₂ > θ₁ path, θ₂ decreases while θ₁ increases, as time progresses. On the θ₁ > θ₂ path, θ₁ decreases while θ₂ initially increases and then decreases, as time evolves. However, our experiments (see Fig. 2(D)) show that a droplet pair can instantaneously jump from one pathway to the other, which manifests in experiments as a droplet-shape switching event and flow reversal in the connecting channel.

4.3 Stability of an interconnected droplet pair

The typical timescales of fluid transport between droplets is much less than the evaporative lifetime of the droplets. Indeed, in the limit where evaporation is a much slower process than the pumping flow, we can consider the dynamics on the short timescale associated with fluid transportation between droplets, i.e., for a constant volume, and neglect the slow dynamics associated with evaporation. The instantaneous configuration of the droplet pair can thus be considered to be a quasi-steady state, corresponding to specific droplet volumes and a corresponding set of contact angles; as the droplets evaporate, the system steps through a sequence of quasi-steady states that are determined by the instantaneous total volume of the system.

To assess the stability of different droplet-pair configurations for a given total volume, V_tot = V₁ + V₂, we note that, on the timescale of the pumping flow, the equations that govern the dynamics of the contact angle can be expressed as:

Eqn (12) has two sets of stationary solutions: one solution set is given by P(θ₁,a₁,γ,ρ,g) − P(θ₂,a₂,γ,ρ,g) = 0, which corresponds to a solution governed by the exchange of fluid between droplets; the other set of stationary solutions, θ₁ = θ₂ = π, represents spherical droplets, which is not relevant to our experiments. To determine the stability of stationary points, we perturb the relevant set of solutions and numerically calculate eigenvalues of the corresponding Jacobian matrix; stability is inferred from the sign of the least negative eigenvalue.

We present the findings of our linear stability analysis of droplet-pair configurations in a parameter space based on the difference in contact angle Δθ = θ₁ − θ₂ and the total volume V_tot = V₁ + V₂ in Fig. 8. Stable solutions are represented by solid lines while the unstable solution is represented by the thick dashed line.

Fig. 8 Quasi-static stability of a pair of connected evaporating droplets. Bifurcation diagrams of the difference in contact angle Δθ = θ₁ − θ₂ measured as a function of the total volume V_tot = V₁ + V₂ for (A) equal pillars, a₁ = a₂ = 0.5 mm, manifests a pitchfork bifurcation, (B) unequal pillars, a₁ = 0.75 mm and a₂ = 0.5 mm, manifests an imperfect pitchfork bifurcation, and (C) varying ratio of pillar radii in the range 1<a₁/a₂<3.5 demonstrates the unfolding of the pitchfork bifurcation for increasing a₁/a₂. Results of a numerical stability analysis of the solutions of the governing equation on the timescale of pumping i.e., without evaporation (eqn (12)), are represented by black lines; solid lines indicate stable solutions, while dashed lines indicate unstable solutions. (D) Volume V^crit_tot, and difference in contact angle Δθ_crit (inset), at the limit-point calculated as a function of imperfection parameter a₁/a₂. In A and B, experimental data (for the same parameters as indicated in Fig. 3(C), (D) and 4(C), (D), respectively) is represented by markers where the colour represents time, as indicated in the legend colormap, and demonstrates how the system evolves through a series of stable, quasi-static configurations because of evaporation.

In the case of pillars of equal radii, a₁ = a₂, we uncover a supercritical pitchfork bifurcation, see Fig. 8(A), that breaks droplet-shape symmetry. For low total volumes, one stable fixed point exists and corresponds to a symmetric droplet-pair configuration with Δθ = 0. As the total volume increases, the system bifurcates from having one fixed point to three fixed points; the branching point, V^crit_tot ≈ 4πa³/3, is a consequence of both droplets having constant contact area and hence a nonlinear pressure–volume relation with maxima for θ ≈ π/2 (as discussed in Section 4.1). In the triple-valued region, asymmetric configurations are stable, and correspond to an unequal droplet pair with pumping flow from the larger to the smaller droplet, while the symmetric solution becomes unstable. This implies that, as time evolves and volume decreases because of evaporation, droplets regain symmetry (i.e., smoothly transition from an asymmetric to a symmetric configuration) once V_tot ≤ V^crit_tot, before evaporating out together.

Breaking the symmetry of the system, by using pillars of unequal radii, leads to an imperfect pitchfork bifurcation, see Fig. 8(B). The upper branch now transitions smoothly down to V_tot = 0 while the stable lower branch meets the unstable symmetric solution at a finite V^crit_tot. As such, a droplet pair that begins on the disconnected, stable asymmetric branch must jump catastrophically onto the connected asymmetric branch once evaporation reduces the system volume below this critical value. Indeed, the droplet-shape switching and accompanying flow reversal, shown in Fig. 2(D), 5(D) and 6, are manifestations of the discontinuous jump from one stable branch, that corresponds to the smaller droplet having larger contact area, to the other stable branch, that corresponds to the larger droplet having larger contact area, which occurs as evaporation reduces the total volume of the system.

Our stability analysis shows the unfolding of this pitchfork bifurcation is governed by the ratio of the two pillar radii, a₁/a₂, see Fig. 8(C). Increasing this imperfection parameter causes the size of the disconnection to increase. Physically, in the range 1 < a₁/a₂ ≤ π, the limit point corresponds to the total volume of the droplet pair when the droplet with larger contact area is hemispherical, with maximal local pressure for θ₁ ≈ π/2, while the droplet with smaller contact area has a non-wetting configuration with a volume that supplies the same pressure, such that V^crit_tot ≈ 2πa₁/3 + V₂(P₂ = P₁(θ₁ ≈ π/2), θ > π/2), as shown in Fig. 8(D). For a₁/a₂ > π, only one stable branch exists and that corresponds to a single droplet residing on the larger pillar.

5 Conclusions & discussion

We have shown that the dynamics of an interconnected droplet pair is governed by two timescales; a short timescale associated with the pumping flow that exchanges fluid between the two drops, and is determined by the pressure differential across the connecting channel, and a long timescale associated with the loss of fluid through evaporation. Initially, the two droplets rapidly adjust their volumes (on a timescale t ≲ 1 s) via a pressure-driven pumping flow and attain a stable configuration. Evaporation then causes the droplet pair to slowly step-through a sequence of quasi-stationary states determined by the instantaneous volume of the system. If the two droplets have the same contact area, a₁ = a₂ = a, and total volume V_tot ≳ 4πa³/3, then the droplets take on an asymmetric configuration with the larger droplet feeding the smaller droplet, so the larger droplet shrinks due to evaporative losses while the smaller droplet maintains approximately the same size and shape. Once V_tot ≲ 4πa³/3, the two droplets take on the same shape, the pumping flow between them reduces to zero, and they evaporate out in unison. The evaporation-driven time-evolution of the droplet pair, from an asymmetric to a symmetric configuration, is underpinned by a perfect pitchfork bifurcation (of quasi-static nature) in the case of droplets of equal contact area. Introducing pillars of unequal radii, and hence contact area, breaks the symmetry of the system and leads to an unfolding of the bifurcation with the distance between the connected and disconnected branches increasing with imperfection parameter a₁/a₂. In experiments, jumping from the disconnected to the connected branch manifests as a sudden droplet-shape switching which changes the sign of the pressure difference between channel ends and causes the flow in the connecting channel to reverse.

Our findings permit engineering control of flow reversal in capillary micropumps via platform design, which could be harnessed for bidirectional transport in point-of-care applications,¹¹ and a framework for avoiding flow reversal in systems where a unidirectional flow is instead required throughout the entire evaporation process. Future research directions include examination of interactions between multiple evaporating droplets connected via a network of channels, and the influence of the vapour phase on the interaction of nearby droplets connected via a short channel.

We finish by drawing a parallel between the behaviour of interconnected droplets and that of two connected rubber balloons²⁹. In this canonical example of elasticity and stability,^30–32 two neo-Hookian rubber balloons are inflated to the same size and connected to the two ends of a single pipe with a closed valve. Once the valve is opened the two balloons can exchange air and, if filled to an intermediate size, the initial configuration is unstable and the symmetry is broken instantly as air flows from one balloon to the other until a stable equilibrium is reached with an asymmetric configuration comprising one smaller and one larger balloon. The two-balloon demo is analogous to our droplet experiments performed on equal pillars since, in both cases, multiplicity arises from nonlinear pressure–volume relations and equilibrium states can be achieved for both unequal and equal shapes of finite (i.e., non-zero) radii. However, in the balloon experiment, multiplicity emerges from the nonlinear pressure–volume relation that is typically a consequence of both geometric and constitutive (i.e., material) nonlinearities, while droplets dynamics are solely governed by a geometric nonlinearity between volume and contact angle. We postulate that a pair of connected deflating balloons would traverse an analogous bifurcation pathway as is exhibited by interconnected droplets as they evaporate; balloon shapes would be regulated by air flow between balloons, while the slow volume loss of the entire system would instead be controlled by the deflation rate. Similarly, our findings on unequal pillars suggest that the underlying symmetry of the two balloons could be broken by using balloons with different intrinsic radii of curvature, for example.

Recently, the inflation of connected balloons has been proposed as a means of generating interactions between individual hysteretic elements (or ‘hysterons’) that can switch between states in a manner that depends on the path history of their states.³³ Interacting mechanical hysterons can store and process information and therefore be programmed to function as state machines.³⁴ Pushing further the analogy between our interconnected droplets and connected balloons, we propose that interconnected droplets are an example of interacting fluidic hysterons. Hence multiple interconnected droplets have the potential to exhibit coupled switching orders with evaporation-rate dependent pathways and memory effects, and could even behave as a (short-lived) state machine. Indeed, our results demonstrate that the transition pathways of these interacting fluidic hysterons can be programmed, in a similar vein to microfluidic bubble logic,³⁵via careful engineering of system properties; in particular, through manipulation of droplet contact area, the relative humidity of the environment and the volatility of the droplets.

Conflicts of interest

The authors declare no conflict of interest and the funders had no role in the design of the study; in the collection, analyses, or interpretation of data; in the writing of the manuscript, or in the decision to publish the results.

Data availability

The reported experimental data are available in the supplementary information (SI). Supplementary information is available. See DOI: https://doi.org/10.1039/d5sm00738k.

Acknowledgements

CR acknowledges support from the China Scholarship Council (grant no. 202106230077). FB acknowledges support from the Royal Society (URF/R1/211730 & RF/ERE/210192).

Notes and references

1. V. Narayanamurthy, Z. Jeroish, K. Bhuvaneshwari, P. Bayat, R. Premkumar, F. Samsuri and M. M. Yusoff, RSC Adv., 2020, 10, 11652–11680.
2. E. Walsh, A. Feuerborn, J. Wheeler, A. Tan, W. Durham, K. Foster and P. Cook, Nat. Commun., 2017, 8, 816.
3. F. Nebuloni, P. R. Cook and E. J. Walsh, J. Fluid Mech., 2023, 969, A28.
4. T. De Groot, K. Veserat, E. Berthier, D. Beebe and A. Theberge, Lab Chip, 2016, 16, 334–344.
5. S. Berry, J. Lee, J. Berthier, E. Berthier and A. Theberge, Anal. Methods, 2019, 11, 4528–4536.
6. E. Berthier, A. M. Dostie, U. N. Lee, J. Berthier and A. B. Theberge, Anal. Chem., 2019, 91, 8739–8750.
7. A. L. McPherson and G. M. Walker, Microfluid. Nanofluid., 2010, 9, 897–904.
8. G. M. Walker and D. J. Beebe, Lab Chip, 2002, 2, 131–134.
9. E. Berthier and D. J. Beebe, Lab Chip, 2007, 7, 1475–1478.
10. I.-J. Chen, E. C. Eckstein and E. Lindner, Lab Chip, 2009, 9, 107–114.
11. S. Xing, R. S. Harake and T. Pan, Lab Chip, 2011, 11, 3642–3648.
12. S. Moradi Mehr, M. A. Charsooghi, L. Businaro, M. Habibi and A.-R. Moradi, AIChE J., 2023, 69, e17847.
13. T. Moragues, D. Arguijo, T. Beneyton, C. Modavi, K. Simutis, A. R. Abate, J.-C. Baret, A. J. deMello, D. Densmore and A. D. Griffiths, Nat. Rev. Methods Primers, 2023, 3, 32.
14. N. S. Lynn and D. S. Dandy, Lab Chip, 2009, 9, 3422–3429.
15. C. Nie, A. Frijns, R. Mandamparambil and J. den Toonder, Biomed. Microdevices, 2015, 17, 47.
16. X.-M. Chen, Y.-J. Li, D. Han, H.-C. Zhu, C.-D. Xue, H.-C. Chui, T. Cao and K.-R. Qin, Micromachines, 2019, 10, 457.
17. G. Castanet, L. Perrin, O. Caballina and F. Lemoine, Int. J. Heat Mass Transfer, 2016, 93, 788–802.
18. A. J. D. Shaikeea and S. Basu, Langmuir, 2016, 32, 1309–1318.
19. A. Wray, B. Duffy and S. Wilson, J. Fluid Mech., 2020, 884, A45.
20. H. Masoud, P. D. Howell and H. A. Stone, J. Fluid Mech., 2021, 927, R4.
21. S. K. Wilson and H.-M. D’Ambrosio, Annu. Rev. Fluid Mech., 2023, 55, 481–509.
22. J. Ju, J. Y. Park, K. C. Kim, H. Kim, E. Berthier, D. J. Beebe and S.-H. Lee, J. Micromech. Microeng., 2008, 18, 087002.
23. K. Javadi, H. Moezzi-Rafie, V. Goodarzi-Ardakani, A. Javadi and R. Miller, Colloids Surf., A, 2017, 518, 30–45.
24. C. Ren, PhD thesis, University of Manchester, 2025.
25. M. Arogeti and A. Shapiro, Chem. Eng. Sci., 2024, 292, 119999.
26. B. R. Munson, D. F. Young and T. H. Okiishi, Fundamentals of fluid mechanics, Wiley, 1990.
27. R. Picknett and R. Bexon, J. Colloid Interface Sci., 1977, 61, 336–350.
28. F. Schönfeld, K.-H. Graf, S. Hardt and H.-J. Butt, Int. J. Heat Mass Transfer, 2008, 51, 3696–3699.
29. I. Müller and P. Strehlow, Rubber and Rubber Balloons, Paradigms of Thermodynamics, Springer-Verlag, Berlin Heidelberg, 2004.
30. F. Weinhaus and W. Barker, Am. J. Phys., 1978, 46, 978–982.
31. D. Merritt and F. Weinhaus, Am. J. Phys., 1978, 46, 976–977.
32. W. Dreyer, I. Müller and P. Strehlow, Q. J. Mech. Appl. Math., 1982, 35, 419–440.
33. G. Muhaxheri and C. D. Santangelo, Phys. Rev. E, 2024, 110, 024209.
34. M. van Hecke, Phys. Rev. E, 2021, 104, 054608.
35. M. Prakash and N. Gershenfeld, Science, 2007, 315, 832–835.

---