Significance The Cheerios effect is the attraction of solid particles floating on liquids, mediated by surface tension forces. We demonstrate experimentally that a similar interaction can also occur for the inverse case, liquid particles on the surface of solids, provided that the solid is sufficiently soft. Remarkably, depending on the thickness of the solid layer, the interaction can be either purely attractive or become repulsive. A theoretical model, in excellent agreement with the experimental data, shows that the interaction requires both elasticity and capillarity. Interactions between objects on soft substrates could play an important role in phenomena of cell–cell interaction and cell adhesion to biological tissues, and be exploited to engineer soft smart surfaces for controlled drop coalescence and colloidal assembly.

Abstract Solid particles floating at a liquid interface exhibit a long-ranged attraction mediated by surface tension. In the absence of bulk elasticity, this is the dominant lateral interaction of mechanical origin. Here, we show that an analogous long-range interaction occurs between adjacent droplets on solid substrates, which crucially relies on a combination of capillarity and bulk elasticity. We experimentally observe the interaction between droplets on soft gels and provide a theoretical framework that quantitatively predicts the interaction force between the droplets. Remarkably, we find that, although on thick substrates the interaction is purely attractive and leads to drop–drop coalescence, for relatively thin substrates a short-range repulsion occurs, which prevents the two drops from coming into direct contact. This versatile interaction is the liquid-on-solid analog of the “Cheerios effect.” The effect will strongly influence the condensation and coarsening of drops on soft polymer films, and has potential implications for colloidal assembly and mechanobiology.

The long-ranged interaction between particles trapped at a fluid interface is exploited for the fabrication of microstructured materials via self-assembly and self-patterning (1⇓⇓⇓–5) and occurs widely in the natural environment when living organisms or fine particles float on the surface of water (6, 7). In a certain class of capillary interactions, the particles deform the interface because of their shape or chemical heterogeneity (8⇓–10). In this case, the change in interfacial area upon particle–particle approach causes an attractive capillary interaction between the particles. In the so-called Cheerios effect, the interaction between floating objects is mainly due to the change in gravitational potential energy associated to the weight of the particles, which deform the interface while being supported by surface tension (11), and the same principle applies when the interface is elastic (12⇓–14). The name “Cheerios effect” is reminiscent of breakfast cereals floating on milk and sticking to each other or to the walls of the breakfast bowl.

Here, we consider a situation opposite to that of the Cheerios effect, liquid drops deposited on a solid. The solid is sufficiently soft to be deformed by the surface tension of the drops, resulting in a lateral interaction. Recent studies have provided a detailed view of statics of single-drop wetting on deformable surfaces (15⇓⇓⇓–19). The length scale over which the substrate is deformed is set by the ratio of the droplet surface tension γ and the substrate shear modulus G. The deformation can be seen as an elastocapillary meniscus, or “wetting ridge,” around the drop (Fig. 1 A and B). Interestingly, the contact angles at the edge of the drop are governed by Neumann’s law, just as for oil drops floating on water. In contrast to the statics of soft wetting, its dynamics has only been explored recently. New effects such as stick–slip motion induced by substrate viscoelasticity (20, 21) and droplet migration due to stiffness gradients (22) have been revealed. The possibility that elastocapillarity induces an interaction between neighboring drops is of major importance for applications such as drop condensation on polymer films (23) and self-cleaning surfaces (24⇓⇓–27). The interaction between drops on soft surfaces might also provide insights into the mechanics of cell locomotion (28⇓–30) and cell–cell interaction (31).

Fig. 1. The inverse Cheerios effect for droplets on soft solids. Two liquid drops sliding down a soft gel exhibit a mutual interaction mediated by the elastic deformation of the substrate. (A) Drops sliding down a thick elastic layer attract each other, providing a new mechanism for coalescence. (B) Drops sliding down a thin elastic layer (thickness h 0 ) repel each other. (C) Measurement of the horizontal relative velocities Δ v x of droplet pairs, as a function of separation distance d. This measurement quantifies the interaction strength. (D) Sliding velocity of isolated droplets on the thick layer as a function of their volume. These data are used to calibrate the relation between force (gravity) and sliding velocity.

Here, we show experimentally that long-ranged elastic deformations lead to an interaction between neighboring liquid drops on a layer of cross-linked polydimethylsiloxane (PDMS). The layer is sufficiently soft for significant surface tension-induced deformations to occur (Fig. 1). The interaction we observe can be thought of as the inverse Cheerios effect, because the roles of the solid and liquid phases are exchanged. Remarkably, the interaction can be either attractive or repulsive, depending on the geometry of the gel. We propose a theoretical derivation of the interaction force from a free-energy calculation that self-consistently accounts for the deformability of both the liquid drop and the elastic solid.

Experiment: Attraction Versus Repulsion Here, the inverted Cheerios effect is observed with submillimeter drops of ethylene glycol on a PDMS gel. The gel is a reticulated polymer network formed by crosslinking small multifunctional prepolymers—contrary to hydrogels, there is no liquid phase trapped inside the network. The low shear modulus of the PDMS gel gives an elastocapillary length ℓ = γ / G = 0.17 mm sufficiently large to be measurable in the optical domain. The interaction between two neighboring liquid drops is quantified by tracking their positions while they are sliding under the effect of gravity along the surface of a soft solid held vertically. The interaction can be either attractive (Fig. 1A) or repulsive (Fig. 1B): drops on relatively thick gel layers attract each other, whereas drops on relatively thin layers experience a repulsion. The drop–drop interaction induces a lateral motion that can be quantified by measuring the horizontal component of the relative droplet velocity, Δ v x ( Δ v x > 0 implies repulsion). In Fig. 1C, we report Δ v x as a function of the separation d, defined as the shortest distance between the surfaces of the drops. The drops ( R ≃ 0.5 − 0.8 mm ) exhibit attraction when sliding down a thick layer ( h 0 = 8 mm , black curve), whereas they repel on a thin layer ( h 0 = 0.04 mm , red curve). The value of Δ v x is larger at close proximity, signaling an increase in the interaction force. Spontaneous merging occurs where drops come into direct contact. Importantly, these interactions provide a new mechanism for droplet coarsening (or ordering) by coalescence (or its suppression) that has no counterpart on rigid surfaces. The interaction force F can be inferred from the relative velocities between the drops, by using an effective “drag law”, where the drag is due to sliding on the gel. We first calibrate this drag law by considering drops that are sufficiently separated, so that they do not experience any mutual interaction. The motion is purely downward and driven only by a gravitational force F g = M g (inertia is negligible). Fig. 1D shows that the droplet velocity v y approximately scales as F g 2 . This force–velocity calibration curve is in good agreement with viscoelastic dissipation in the gel, based on which one expects the scaling law (21): v ∼ ℓ τ ( F 2 π R γ ) 1 / n . [1]Here, n is the rheological exponent that emerges from the scale invariance of the gel network (32⇓–34), and τ is a characteristic timescale. The parameter values n ≃ 0.61 and τ ≃ 0.68 s are calibrated in a rheometer (SI Materials and Methods). Eq. 1 is valid for v below the characteristic rheological speed, ℓ / τ . Our approach is justified here because ℓ / τ ∼ 0.25 mm / s for the silicone gel, whereas the reported speeds reach at most ∼ 100 nm / s . The large viscoelastic dissipation in the gel exceeds the dissipation within the drop by orders of magnitude, and explains these extremely slow drop velocities observed experimentally (21, 35). In this case, it was also shown that all of the dissipation occurs in a very narrow region around the wetting ridge (21). Therefore, the dynamic substrate deformation approaches the corresponding static deformation beyond a distance v τ ≲ 60 nm from the contact line. The force–distance relation for the inverse Cheerios effect can now be measured directly using the independently calibrated force–velocity relation (Fig. 1D and Eq. 1). By monitoring how the trajectories are deflected with respect to the downward motion of the drops, we obtain F (see Materials and Methods for additional details). Despite the different origins of calibration and interaction forces, both are balanced by the same dissipative mechanism because the dissipative viscoelastic force is nearly perfectly localized at the contact line (21), which corroborates the validity of our calibration routine. The key result is shown in Fig. 2, where we report the interaction force F as a function of distance d. Fig. 2A shows experimental data for the attractive force ( F < 0 ) between drops on thick layers (black dots), together with the theoretical prediction outlined below. Movie S1 shows an example of attractive drop–drop interaction. The attractive force is of the order of micronewtons, which is comparable to both the capillary force scale γ R and the elastic force scale G R 2 . The force decreases for larger distance and its measurable influence was up to d ∼ R . Fig. 2. Measured interaction force F (symbols) as a function of their separation d, compared with the 3D theory (red lines, no adjustable parameters). (A) Attraction on a thick elastic layer ( h 0 ≈ 8 mm ≫ R ≫ ℓ ). (B) Repulsion and attraction on a thin layer ( R ≫ ℓ ≳ h 0 ≈ 40 μ m ). Each data point represents an average over ∼10 realizations, with the error bars giving the SD. Measurements are based on pairs of ethylene glycol drops whose radii are in the range R ∼ 0.7 ± 0.1 mm . The elastic substrate has a static shear modulus of 0.28 kPa . Fig. 2B shows the interaction force between drops on thin layers. The dominant interaction is now repulsive ( d ≳ h 0 ) (Movie S2). Intriguingly, we find that the interaction is not purely repulsive, but displays an attractive range at very small distance. It is possible to access this range experimentally in case the motion of the individual drops are sufficiently closely aligned (Movie S3). The “neutral” distance for which the interaction force changes sign appears when the separation is comparable to the substrate thickness h 0 , suggesting that the key parameter governing whether the drops attract or repel is the thickness of the gel.

Mechanism of Interaction: Rotation of Elastic Meniscus We explain the attraction versus repulsion of neighboring drops by computing the total free energy ℰ of drops on gel layers of different thicknesses. The interaction force between the drops is equal to the energy gradient with respect to the separation, − ∂ ℰ / ∂ d , which in the experiment is balanced by the forces due to viscoelastic dissipation in the vicinity of the contact line. In contrast to the normal Cheerios effect, which involves two rigid particles, both the droplet and the elastic substrate are deformable, and their shapes will change upon varying the distance d. Hence, the interaction force involves both the elastic and the surface tension contributions to the free energy. The free energy emerges from self-consistently computed shapes of the drops and elastic deformations. To reveal the mechanism of interaction, we first consider 2D drops, for which the free energy can be written as follows: ℰ [ h ] = ℰ e ℓ [ h ] + ∫ dry d x γ S V 1 + h ′ 2 + ∫ wet d x [ γ 1 + ℋ ′ 2 + γ S L 1 + h ′ 2 ] . [2]The geometry is sketched in Fig. 3 B and C, and further details are given in Supporting Information and Fig. S1. The elastic energy ℰ e ℓ is a functional of the profile h ( x ) describing the shape of the elastic solid: the functional explicitly depends on the layer thickness and is ultimately responsible for the change from attraction to repulsion. The function ℋ ( x ) represents the shape of the liquid–vapor interface. The integrals in Eq. 2 represent the interfacial energies; they depend on the surface tensions γ, γ S L , γ S V associated with the liquid–vapor, solid–liquid, and solid–vapor interfaces, respectively. For the sake of simplicity, we ignore here the possibility of a dependence of surface energy on the elastic strain. In absence of the Shuttleworth effect (36), surface stress and surface energy are equal. Fig. 3. Mechanism of interaction between two liquid drops on a soft solid. (A) Deformation h ( x ) induced by a single droplet on a thick substrate. The zoom near the contact line illustrates that the contact angles satisfy the Neumann condition. (B) A second drop placed on a thick substrate experiences a background profile due to the deformation induced by the neighboring drop on the right. This background profile is shown in red. As a consequence, the solid angles near the elastic meniscus rotate by an angle φ (see zoom). This rotation perturbs the Neumann balance, yielding an attractive force f. In the experiment, this force is balanced by the dissipative force due to the viscoelastic deformation of the wetting ridge. (C) On a thin substrate the single-drop profile yields a nonmonotonic elastic deformation. The zoom illustrates a rotation φ of the Neumann triangle in the opposite direction, leading to a repulsive interaction. Fig. S1. Problem setup. Two liquid droplets on an elastic half space. Energy is to be minimized with respect to ℋ ( x ) , h s v ( x ) , h s l ( x ) , x i , and x o . The liquid angles θ i and θ o follow from the minimization. The distance d between the adjacent contact lines is enforced with a Lagrange multiplier f acting on ( x i − d / 2 ) . At equilibrium, the droplet shapes can be found by analyzing the changes in the free energy upon variations of the functions h ( x ) and ℋ ( x ) . The variation of the contact line positions provide the relevant boundary conditions (19). However, when two drops are separated by a finite distance d, the drops are not at equilibrium: the gradient in free energy results into an overdamped motion in which changes in free energy are dissipated in the solid. To compute the interaction force f (per unit length in the 2D model), one therefore needs to consider the work done by the dissipative force −f that we can assume to be localized near the inner contact line. This allows one to determine the interaction force f = − ∂ ℰ / ∂ d , with the convention that attractive forces correspond to f < 0 (see Supporting Information and Fig. S2 for details). Fig. S2. Force–distance curve for 2D droplets on elastic half space (log–log scale), for γ / ( G R ) = 1 / 10 , γ / γ s = 10 . The Inset shows solid profiles at the indicated distances. The energy minimization reveals the mechanism of drop–drop interaction: the interaction force f appears in the boundary condition for the contact angles, f = γ ⁡ cos ⁡ θ + γ S L ⁡ cos ⁡ θ S L − γ S V ⁡ cos ⁡ θ S V , [3]where the angles are defined in Fig. 3. Eq. 3 can be thought of as an “imbalance” of the static Neumann boundary condition. The resulting interaction force due to the elastocapillary energy gradient is balanced by the dissipation due to the viscoelastic nature of the substrate. For a single droplet, the contact angles satisfy Neumann’s law, which is Eq. 3 with f = 0 (Fig. 3A). On a thick elastic layer, the overall shape of the wetting ridge is of the following form (18, 19): h ( x ) ∼ γ G Ψ ( x γ s / G ) , [4]where the horizontal scale is set by the elastocapillary length based on the solid surface tension γ s . The origin of f can be understood from the principle of superposition. Due to the substrate deformation of a single drop, a second drop approaching the first one will see a surface that is locally rotated by an angle φ ∼ h ′ ∼ γ / γ s . The elastic meniscus near the inner contact line of this approaching drop will correspondingly be rotated by an angle φ (Fig. 3B). Importantly, changes in the liquid angle θ scale as ∼ h / R ∼ γ / ( G R ) , which for large drops can be ignored. As a consequence, this meniscus rotation induces a net resultant surface tension forces according to Eq. 3, which is balanced by the dissipative force f due to the viscoelastic nature of the substrate (21). For small rotations, one obtains f ≃ γ φ , where φ follows from the single-drop deformation (Eq. 4). There is no resultant interaction force from the stress below the drop, which, due to deformability of the drop, results only in a uniform pressure on the solid-liquid interface. The inverted Cheerios effect is substantially different from the Cheerios effect between two particles floating at the surface of a liquid. Apart from the drop being deformable, we note that the energy driving the interaction is different in the two cases: whereas the liquid interface shape is determined by the balance between gravity and surface tension in the Cheerios effect, the solid shape is determined by elastocapillarity in the inverted Cheerios effect. Another difference is the mechanism by which the interaction is mediated. The Cheerios effect is primarily driven by a change in gravitational potential energy, which implies a vertical displacement of particles: a heavy particle slides downward, like a bead on a string, along the deformation created by a neighboring particle (11). A similar interaction was recently discussed for rigid cylinders that deform an elastic surface due to gravity (14). In contrast, the inverted Cheerios effect discussed here does not involve gravity and can be totally ascribed to elastocapillary tilting of the solid interfaces, as shown in Fig. 3. The rotation of contact angles explains why the drop–drop interaction can be either attractive or repulsive. On a thick substrate, the second drop experiences solid contact angles that are rotated counterclockwise, inducing an attractive force (Fig. 3B). In contrast, on a thin substrate, the elastic deformation induced by the second drop has a nonmonotonic profile h ( x ) . This is due to volume conservation of the substrate: lifting the gel near the contact line creates a depression at larger distances (Fig. 3C). The rotation of the contact angles thus changes sign, and, accordingly, the interaction force changes from attractive to repulsive. The relevant length scale for this phenomenon is set by the layer thickness h 0 .

Three-Dimensional Theory The extension of the theory to three dimensions is straightforward and allows for a quantitative comparison with the experiments. For the 3D case, we compute the shape of the solid numerically, by first solving for the deformation field induced by a single drop using an axisymmetric elastic Green's function (18). Adding a second drop on this deformed surface gives an intricate deformation that is shown in Fig. 4 and Fig. S3. The imbalance of the Neumann law applies everywhere around the contact line: the background deformation induces a rotation of the solid contact angles around the drop. According to Eq. 3, these rotations result into a distribution of force per unit length contact line f → = f ( β ) e → r , where e → r is the radial unit vector and β is the azimuthal angle (Fig. 4). The resultant interaction force F → is obtained by integration along the contact line, F → = R ∫ d β f → ( β ) (see Supporting Information and Fig. S4 for details). By symmetry, this force is oriented along the line connecting the two drops. Fig. 4. Three-dimensional calculation of interface deformation for a pair of axisymmetric drops. The elastocapillary meniscus between the two drops is clearly visible, giving a rotation of the contact angle around the drop. The total interaction force F → is obtained by integration of the horizontal force f → (indicated in red) and is related to the free-energy gradient associated with a change in separation between the drops. Parameter values are ℓ / R = 0.1 , γ / γ s = 1 . Fig. S3. Problem setup for axisymmetric drops on an elastic half space. The solid shape is shown for a distance d = 0.3 R between the droplets and γ / γ s = 10 . x and y are in units of R, and z is in units of γ / G . Fig. S4. Force–distance curve for axisymmetric droplets on a thin elastic layer. (Inset) f r ( β ) at the contact line of a droplet on a thin layer in the presence of a second droplet at d / h 0 = 2 . h 0 / R = 1 / 20 , γ / ( G R ) = 1 / 10 , γ / γ s = 10 . The interaction force obtained by the 3D theory is indicated by the red curves in Fig. 2 A and B. The theory gives an excellent description of the experimental data without adjustable parameters. The quantitative agreement indicates that the interaction mechanism is indeed caused by the rotation of the elastic meniscus.

Discussion In summary, we have shown that liquid drops can exhibit a mutual interaction when deposited on soft surfaces. The interaction is mediated by substrate deformations, and its direction (repulsive versus attractive) can be tuned by varying the thickness of the layer. The measured force/distance relation is in quantitative agreement with the proposed elastocapillary theory. The current study reveals that multiple “pinchings” of an elastic layer by localized tractions γ lead to an interaction having a range comparable to γ / G . The key insight is that interaction emerges from the rotation of the elastic surface, providing a generic mechanism that should be applicable to a wide range of objects interacting on soft media. Our model provides general concepts that are applicable to a wide range of experimental settings, whenever objects exert dipolar or quadrupolar forces on their substrate [the integral force must vanish in this case, however, as is the case, e.g., for cells (37)]. The length scale of interaction is governed by the ratio of two quantities: the force (per unit length) γ, and the substrate shear modulus G. This can range from nanometers for small forces or stiff substrates, to hundreds of microns for strong forces or soft substrates. In biological settings, elastocapillary interactions may play a role in cell–cell interactions, which are known to be sensitive to substrate stiffness (31). One example would be stem cell aggregates that interact with their extracellular matrix (38). In addition, the elastic interaction could also play a role in cell–extracellular matrix interactions, as a purely passive force promoting aggregation between anchor points on the surface of adhered biological cells. For example, it has been demonstrated that a characteristic distance of about 70 nm between topographical features enables the clustering of integrins. These transmembrane proteins are responsible for cell adhesion to the surrounding matrix, mediating the formation of strong anchor points when cells adhere to substrates (39, 40). Assuming that the topographical features “pinch” the cell with a force likely comparable to the cell’s cortical tension, which takes values in the range 0.1–1 mN/m (41⇓⇓–44), and an elastic modulus of 103–104 Pa in the physiological range of biological tissues (45), one predicts a range of interaction consistent with observations. More generally, substrate-mediated interactions could be dynamically programmed using the responsiveness of many gels to external stimuli (pH, temperature, electric fields). Possible applications range from fog harvesting and cooling to self-cleaning or anti-fouling surfaces, which rely on controlling drop migration and coalescence. The physical mechanisms revealed here, in combination with the fully quantitative elastocapillary theory we propose, pave the way for new design strategies for smart soft surfaces.

Materials and Methods Supporting Information provides further technical information, the derivation underlying Eq. 3, and the numerical scheme used for the calculations of Figs. 2 and 4. Movies S1–S3 show typical experiments of drop–drop interactions. Substrate Preparation. The two prepolymer components (Dow Corning; CY52-276 A and B) were mixed in a ratio of 1.3:1 (A:B). Thick elastic layers (∼8 mm) were prepared in Petri dishes (diameter, ∼90 mm). Thin layers (∼40 μm) were prepared by spin-coating the gel onto silicon wafers. The thickness was determined by color interferometry. See SI Materials and Methods for details on substrate curing and rheology (Fig. S5). Fig. S5. Spectra of storage ( G ′ ) and loss modulus ( G ″ ) of the substrates. Curves for −10 °C and 80 °C were shifted to obtain a room temperature master curve. Determining the Interaction Between Drops. Droplets of ethylene glycol (V ∼ 0.3–0.8 μL) were pipetted onto a small region near the center of the cured substrate. The sample was then mounted vertically so that gravity acts along the surface ( − y direction; compare Fig. 1 A and B). The droplets were observed in transmission (thick layers) or reflection (thin layers) with collimated illumination, using a telecentric lens (JenMetar 1×) and a digital camera (pco 1200). Images were taken every 10 s. The contours of the droplets were determined by a standard correlation technique. At large separation, droplets move downward due to gravity. The gravitational force on each droplet is proportional to its volume. The relation between force and velocity follows the same power law as the rheology, as was explained recently (21). Individual droplets have different volumes and move with different velocities. Thus their distances change with time. Whenever two droplets approach each other, their trajectories change due to their interaction. Drops on thick substrates (Fig. 1C, black) attract and eventually merge. On a rigid surface, these droplets would not have merged. The opposite holds for droplets on thin layers (red): the droplets repel each other, which prevents coalescence. To determine the interaction forces, we first evaluate the velocity vector of each individual droplet. The droplets move in a quasi-stationary manner, and the total force vector acting on each droplet is aligned with its velocity vector. The magnitude of the total force is obtained through calibration from the data shown in Fig. 1D. The interaction force is obtained by subtracting the gravitational force from the total force. Fig. 2A shows data from nine individual droplet pairs, corresponding to different times and different locations on the substrate. Fig. 2B shows data from 18 different droplet pairs. The raw data have been averaged over distance bins, taking the SD as error bar.

Supporting Information describes technical details of the theory by which we compute the drop–drop interaction as well as experimental information. We first discuss the 2D variational problem and numerical solution, and the scheme is then extended to axisymmetric 3D droplets. SI Materials and Methods describe details of the experiment, including the rheological calibration.

Interaction of 2D Droplets Free Energy. The problem setup is shown in Fig. S1. Two liquid drops of radius R sit on a soft, elastic substrate with shear modulus G. The distance between liquid drops is d. The soft substrate has surface free energies γ s v (dry part) and γ s l (wet part). γ is the liquid surface tension. The proximal contact lines of the two drops are at x i and − x i , whereas the distal contact lines are at x o and − x o . x i = d / 2 and x o = d / 2 + 2 R . θ i and θ o are the inner and outer liquid contact angles and are unknown at this point. The deformed profile of the substrate is defined by h s l ( x ) and h s v ( x ) , for wet and dry parts, respectively. Due to the extremely slow motion of the droplets, all dissipative effects are localized at the contact line (21). We use a quasistationary approximation and incorporate the total dissipative force into a Lagrange multiplier f, acting on the contact line distance between the drops. Owing to the left–right symmetry, we use the energy of the right half of the system, so that the total free energy of the system is 2 ℰ : ℰ = ℰ e l [ h s v , h s l ] + γ ∫ x i x o 1 + ℋ ′ 2 d x + γ s v ∫ 0 x i ( 1 + h s v ′ 2 ) d x + γ s l ∫ x i x o ( 1 + h s l ′ 2 ) d x + γ s v ∫ x o ∞ ( 1 + h s v ′ 2 ) d x + Λ 1 ( ℋ ( x i ) − h s v ( x i ) ) + Λ 2 ( ℋ ( x i ) − h s l ( x i ) ) + Λ 3 ( ℋ ( x o ) − h s v ( x o ) ) + Λ 4 ( ℋ ( x o ) − h s l ( x o ) ) + f ( x i − d / 2 ) + P ( V − ∫ x i x o ( ℋ − h s l ) d x ) . [S1]The first term, ℰ e l , describes the elastic energy. Its explicit form will be discussed later. Terms with prefactors γ, γ s v , and γ s l are the surface energies of the drop and the dry and the wet part of the solid, respectively. This form is presented as Eq. 2 in the main text. The remaining terms are Lagrange multipliers to enforce the various boundary conditions on the three phase contact lines: • Λ 1 to Λ 4 enforce contact between the three interfaces ℋ , h s v , and h s l at x i and x o , and thereby also the continuity of the substrate surface.

• f enforces the distance d between the inner contact lines of the two droplets. Hence, f quantifies the interaction force. Positive value of f corresponds to a repulsive interaction.

• P is the Laplace pressure inside the drop, which enforces a drop volume V. We solve for stationary ℰ at fixed d. The variation of d gives δ ℰ / δ d = − f / 2 , which implies that f is indeed the interaction force (the factor 1/2 drops out because ℰ is only one-half of total free energy of the system). Variation of ℋ ( x ) . δ ℰ = γ ∫ x i x o ℋ ′ 1 + ℋ ′ 2 δ ℋ ′ d x + Λ 1 δ ℋ ( x i ) + Λ 2 δ ℋ ( x i ) + Λ 3 δ ℋ ( x o ) + Λ 4 δ ℋ ( x o ) − P ∫ x i x o δ ℋ ( x ) d x = ( Λ 1 + Λ 2 − γ ℋ ′ 1 + ℋ ′ 2 ) δ ℋ ( x i ) + ( Λ 3 + Λ 4 + γ ℋ ′ 1 + ℋ ′ 2 ) δ ℋ ( x o ) − ∫ x i x o ( P + γ ℋ ″ ( 1 + ℋ ′ 2 ) 3 / 2 ) δ ℋ ( x ) d x . [S2]Because δ ℋ is arbitrary, and δ ℰ = 0 at equilibrium, we obtain the following relations for the Lagrange multipliers Λ 1 to Λ 4 : Λ 1 + Λ 2 = γ ℋ ′ ( x i ) 1 + ℋ ′ ( x i ) 2 = γ ⁡ sin ⁡ θ i , Λ 3 + Λ 4 = γ − ℋ ′ ( x o ) 1 + ℋ ′ ( x o ) 2 = γ ⁡ sin ⁡ θ o . [S3]In addition, the integrand in Eq. S2 should vanish in the domain ( x o , x i ) , which gives the differential equation for the drop shape: γ ℋ ″ ( 1 + ℋ ′ 2 ) 3 / 2 = − P . [S4] Variation of h s v ( x ) and h s l ( x ) . The Lagrange multipliers Λ i enforce the continuity of the solid profile at x i and x o . With h ( x ) = { h s v ( x ) for x < x i h s l ( x ) for x i < x < x o , h s v ( x ) for x > x i [S5]we can write for the elastic free energy: ℰ e l = ∫ 0 ∞ σ n ( x ) h ( x ) d x , [S6]where σ n ( x ) is the normal stress at the surface of the elastic medium. The variation gives the following: δ ℰ = ∫ 0 ∞ σ n ( x ) δ h ( x ) d x + P ∫ x i x o δ h s l d x + γ s v ∫ 0 x i h s v ′ 1 + h s v ′ 2 δ h s v ′ d x + γ s v ∫ x o ∞ h s v ′ 1 + h s v ′ 2 δ h s v ′ d x + γ s l ∫ x i x o h s l ′ 1 + h s l ′ 2 δ h s l ′ d x − Λ 1 δ h s v ( x i ) − Λ 2 δ h s l ( x i ) − Λ 3 δ h s v ( x o ) − Λ 4 δ h s l ( x o ) = ∫ 0 ∞ σ n ( x ) δ h ( x ) d x + P ∫ x i x o δ h s l d x + γ s v [ h s v ′ 1 + h s v ′ 2 δ h s v | 0 x i − ∫ 0 x i h s v ″ ( 1 + h s v ′ 2 ) 3 / 2 δ h s v d x ] + γ s l [ h s l ′ 1 + h s l ′ 2 δ h s l | x i x o − ∫ x i x o h s l ″ ( 1 + h s l ′ 2 ) 3 / 2 δ h s l d x ] + γ s v [ h s v ′ 1 + h s v ′ 2 δ h s v | x o ∞ − ∫ x o ∞ h s v ″ ( 1 + h s v ′ 2 ) 3 / 2 δ h s v d x ] − Λ 1 δ h s v ( x i ) − Λ 2 δ h s l ( x i ) − Λ 3 δ h s v ( x o ) − Λ 4 δ h s l ( x o ) . [S7]Due to symmetry, h s v ′ ( 0 ) = 0 . In addition, we impose that the deformation h s v (and thus its variation δ h s v ) shall vanish for x → ∞ . Hence, the elastic problem is defined by the normal stress as follows: σ n ( x ) = { γ s v h s v ″ ( 1 + h s v ′ 2 ) 3 / 2 for x < x i − P + γ s l h s l ″ ( 1 + h s l ′ 2 ) 3 / 2 for x i < x < x o . γ s v h s v ″ ( 1 + h s v ′ 2 ) 3 / 2 for x > x o [S8]In analogy to ref. 19, the elastic stress consists of a solid Laplace pressure term proportional to the curvature of the deformed substrate; inside the drop, there is in addition the capillary pressure P. As before, the boundary terms give the following: Λ 1 = γ s v h s v ′ ( x i ) 1 + h s v ′ ( x i ) 2 = γ s v ⁡ sin ⁡ θ s v , Λ 2 = γ s l − h s l ′ ( x i ) 1 + h s l ′ ( x i ) 2 = γ s l ⁡ sin ⁡ θ s l , Λ 3 = γ s v − h s v ′ ( x o ) 1 + h s v ′ ( x o ) 2 = γ s v ⁡ sin ⁡ θ s v , Λ 4 = γ s l h s l ′ ( x o ) 1 + h s l ′ ( x o ) 2 = γ s l ⁡ sin ⁡ θ s l . [S9]Combined with Eq. S3, this gives the vertical boundary conditions at the two contact lines, γ ⁡ sin ⁡ θ i , o = γ s v ⁡ sin ⁡ θ s v + γ s l ⁡ sin ⁡ θ s l . [S10] Variation of x i and x o . The key equation that provides the interaction force, Eq. 3 of the main text, is obtained from the variation of x i : δ ℰ δ x i = − γ 1 + ℋ ′ ( x i ) 2 + γ s v 1 + h s v ′ ( x i ) 2 − γ s l 1 + h s l ′ ( x i ) 2 + Λ 1 ( ℋ ′ ( x i ) − h s v ′ ( x i ) ) + Λ 2 ( ℋ ′ ( x i ) − h s l ′ ( x i ) ) + f + P ( ℋ ( x i ) − h s l ( x i ) ) . [S11]The last term drops out as the interfaces meet at the contact lines. Now using Eq. S9 to express Λ 1 and Λ 2 , the equilibrium condition δ ℰ / δ x i = 0 gives the following: f = γ 1 + ℋ ′ ( x i ) 2 − γ s v 1 + h s v ′ ( x i ) 2 + γ s l 1 + h s l ′ ( x i ) 2 , f = γ ⁡ cos ⁡ θ i + γ s l ⁡ cos ⁡ θ s l − γ s v ⁡ cos ⁡ θ s v . [S12]This is Eq. 3 of the main text, which demonstrates that f emerges as an imbalance in the horizontal Neumann condition of the contact angles at x i . For the variation of x o , we obtain the following: δ ℰ δ x o = γ 1 + ℋ ′ ( x o ) 2 + γ s l 1 + h s l ′ ( x o ) 2 − γ s v 1 + h s v ′ ( x o ) 2 + Λ 3 ( ℋ ′ ( x o ) − h s v ′ ( x o ) ) + Λ 4 ( ℋ ′ ( x o ) − h s l ′ ( x o ) ) − P ( ℋ − h s l ) . [S13]Again, with Eq. S9 and δ ℰ / δ x o = 0 : γ 1 + ℋ ′ ( x o ) 2 + γ s l 1 + h s l ′ ( x o ) 2 − γ s v 1 + h s v ′ ( x o ) 2 = 0 , γ ⁡ cos ⁡ θ o + γ s l ⁡ cos ⁡ θ s l − γ s v ⁡ cos ⁡ θ s v = 0. [S14]This is the horizontal Neumann condition at the outer contact line at x o . Numerical Solution. This variations of ℋ ( x ) , h ( x ) and the contact line positions x i , x 0 define the problem of two interacting drops. As usual, the constrained minimization (i.e., imposing the distance between the inner contact lines) only admits solutions for a particular value of the Lagrange multiplier f. By determining numerically the solutions ℋ ( x ) , h ( x ) consistent with the boundary conditions, we thus find a unique interaction force f for a given distance. The spherical cap solution for ℋ ( x ) , together with elementary geometry, allows to express the drop pressure in terms of the contact angles as follows: P = γ sin ⁡ θ i + sin ⁡ θ o 2 R . [S15]This appears as a traction to the elastic layer. The substrate shape h ( x ) is solved using linear elasticity and a Green’s function approach, h ( x ) = ∫ − ∞ ∞ K ( x − x ′ ) σ n ( x ′ ) d x ′ , where K is the Green’s function for the elastic medium. The resulting equation is best solved in Fourier space (15). We use the following convention: f ˜ ( q ) = ∫ − ∞ ∞ f ( x ) e − i q x d x , [S16]to obtain in the limit of small slopes and assuming γ s v = γ s l = γ s : h ˜ ( q ) = σ ˜ n ( q ) γ s q 2 + K ˜ ( q ) − 1 . [S17]Here, K ˜ ( q ) is the Fourier transform of Green’s function: K ˜ ∞ ( q ) = 1 2 G q , [S18] K ˜ h 0 ( q ) = 1 2 G q sinh ( 2 q h 0 ) − 2 q h 0 cosh ( 2 q h 0 ) + 2 ( q h 0 ) 2 + 1 , [S19]where K ˜ ∞ corresponds to an elastic half-space, and K ˜ h 0 to elastic layers with finite thickness h 0 . σ ˜ n ( q ) is the traction that the two droplets exert onto the elastic medium given by Eq. S8. Following the same steps as for the single-drop problem (19), we find the explicit form σ ˜ n ( q ) = 2 γ [ ( sin θ i + sin θ o ) ( cos ⁡ q R − sin q R q R ) cos ( d + 2 R ) q 2 + ( sin θ i − sin θ o ) sin q R sin ( d + 2 R ) q 2 ] , [S20]where we made use of Eq. S15 to eliminate the pressure inside the drop. Note that, within this formulation, the vertical Neumann balance is automatically ensured. The solution for the elastic layer h ( x ) is now fully specified in terms of θ i and θ o . These need to satisfy the two horizontal Neumann conditions (Eqs. S12 and S14), but this introduces another unknown f. The problem is closed by the requirement that the heights match at the contact lines, h ( x i , o ) = ℋ ( x i , o ) , which selects a solution with unique θ i , o and f for each distance d. These solutions are determined numerically. Fig. S2 shows a force distance curve for typical parameters, together with solid profiles (Inset) at given distances.

Interaction of Axisymmetric Droplets According to ref. 18, the Green’s functions for axisymmetric loads on elastic media are as follows: K ˜ ∞ ( ξ ) = 1 2 G ξ , [S21] K ˜ h 0 ( ξ ) = 1 2 G ξ sinh ( 2 ξ h 0 ) − 2 ξ h 0 cosh ( 2 ξ h 0 ) + 2 ξ 2 h 0 2 + 1 , [S22]for half-space and thin layers, respectively. The displacements due to the traction exerted by a droplet at the origin reads in Fourier space (18): h ˜ ( ξ ) = γ sin ⁡ θ R J 0 ( ξ R ) − 2 ξ J 1 ( ξ R ) K ˜ ( ξ ) − 1 + γ s ξ 2 . [S23]Real-space profiles are obtained by numerical Hankel transformation according to the following: h ( r ) = ∫ 0 ∞ ξ h ˜ ( ξ ) J 0 ( ξ r ) d ξ . [S24] Two Droplets. In the presence of a second droplet, the slope of the wetting ridge of the second droplet will rotate the solid angles at the contact line of the first drop. Therefore, two droplets nearby will not be the equilibrium configuration, but attract or repel, depending on the relative strengths of the forces acting on the contact line. To keep two droplets with circular footprints of radius 1 at a distance d, we impose a constraint on the radial coordinate of the contact line. Associated with this is a continuous Lagrange multiplier f r ( β ) that depends on the azimuthal angle β. f r ( β ) can be interpreted as a virtual traction f → = f r e → r that fixes the contact line position. The horizontal force balance becomes the following: f r ( β ) = γ ⁡ cos ⁡ θ + γ s ( cos ⁡ θ s l ( β ) − cos ⁡ θ s v ( β ) ) . [S25] Let h 2 ( x , y ) be the solid shape for the case of two droplets a and b, where drop a shall be located at the origin of the coordinate system, and drop b shall be located to the right of drop a, centered at x = 2 + d (Fig. S3). The full shape is composed as the superposition of two individual droplets as follows: h 2 ( x , y ) = h ( x , y ) + h ( 2 R + d − x , y ) . [S26]In polar coordinates, this becomes the following: h 2 ( r , β ) = h ( r ) + h ( r 2 ( r , β ) ) , [S27]where r 2 ( r , β ) = ( 2 R + d − r ⁡ cos ⁡ β ) 2 + ( r ⁡ sin ⁡ β ) 2 . In small-slope approximation, the horizontal force balance for drop a becomes the following: f r ( β ) = γ ⁡ cos ⁡ θ + γ s [ ( ∂ r h 2 ( r , β ) | r = 1 + ) 2 − ( ∂ r h 2 ′ ( r , β ) | r = 1 − ) 2 ] , [S28]where the superscripts + and − indicate the slopes on the dry and the wet side of the contact line, respectively. Again, this can be decomposed into a slope discontinuity and a rotation. The discontinuity equals γ / γ s , and the rotational part can be obtained by evaluating the backward transform exactly at the position of the discontinuity (46), and the force balance becomes the following: f r ( β ) = γ ( cos ⁡ θ + h ′ ( R ) + ( 1 − ( 2 R + d ) cos ⁡ β ) h ′ ( r 2 ) r 2 ) . [S29]The interaction force F ( d ) between the droplets is given by the integral of the projection of f r ( β ) onto e → x : F = ∫ 0 ∞ f r ( β ) cos ⁡ β d β . [S30] The calculation for axisymmetric droplets on thin elastic layers is identical except for the convolution kernel. As in the 2D case, the “dimple” in the profile causes a sign change in the slope, which results in a repulsive regime for distances d ≳ h 0 . However, now the contact line of one droplet always crosses the dimple of the other droplet. Therefore, the rotation of the solid angle and the net radial force on the contact line may change sign when plotted against the azimuthal angle (compare Fig. S4, Inset). This leads to a different f ( d ) compared with the 2D model: the sign change in the 3D model appears at smaller distances because a significant fraction of the contact line remains in the repulsive range.

SI Materials and Methods Experiments were performed using ethylene glycol (Sigma; purity, ≥ 98 % ) drops on a silicone gel (Dow Corning; CY52-276). Typical droplet radii were R ∼ 500–800 μm. The two prepolymer components were mixed in a ratio of 1.3:1 (A:B), stirred thoroughly, degassed, and poured into Petri dishes (diameter, ∼90 mm) to make thick (∼8 mm) elastic layers. Thin layers (∼40 μm) were prepared by spin-coating the gel onto silicon wafers. Thicknesses were determined by color interferometry. To ascertain smooth and flat gel surfaces, substrates were cured in two steps: first, the sample was kept for 48 h at room temperature on a stabilized optical table in a dust-free environment. After that, the gels were slowly heated to 90 °C in an oven and cured for 2 h at that temperature. To determine the substrate rheology, part of the mixed prepolymers was cured simultaneously (with the same temperature profile) in a parallel plate geometry ( d = 50 mm ) of a rheometer (Anton Paar MCR502 with Peltier bottom plate and hood). After curing, frequency sweeps with small strains were performed at 20 °C. Measurements are accurately fitted by the following form: G ′ ( ω ) + i G ″ ( ω ) = G [ 1 + ( i τ ω ) n ] , [S31]with exponent n ≃ 0.61 , a timescale τ = 0.68 s, and a long-time shear modulus G = 280 Pa. Combined with γ = 0.048 N/m, this provides the magnitude for elastic deformations dictated by the elastocapillary length ℓ = 0.17 mm . Fig. S5 shows the measured rheology of the substrates used in the experiments. High- and low-temperature measurements are shifted to give a room temperature master curve.

Acknowledgments L.B. acknowledges useful discussion on mechanobiology with Dr. Nuria Gavara, and financial support from European Union Grant CIG 618335. S.K. acknowledges financial support from NWO through VIDI Grant 11304. A.P. and J.H.S. acknowledge financial support from European Research Council Consolidator Grant 616918.

Footnotes Author contributions: S.K., L.B., B.A., and J.H.S. designed research; S.K., A.P., L.A.L., J.H.W., L.B., S.D., and J.H.S. performed research; S.K. and A.P. analyzed data; and S.K., A.P., L.B., B.A., and J.H.S. wrote the paper.

The authors declare no conflict of interest.

This article is a PNAS Direct Submission. K.D.-V. is a guest editor invited by the Editorial Board.

See Commentary on page 7294.

This article contains supporting information online at www.pnas.org/lookup/suppl/doi:10.1073/pnas.1601411113/-/DCSupplemental.