Dynamical Theory of the Inverted Cheerios Effect

TL;DR

This paper develops a dynamical theory for the inverted Cheerios effect, describing how liquid drops interact on deformable substrates by considering elasticity, capillarity, and viscoelasticity, and validates it against experiments.

Contribution

It introduces a comprehensive dynamical model for droplet interactions on deformable substrates, incorporating elasticity, capillarity, and viscoelastic effects, advancing understanding of this phenomenon.

Findings

The theory accurately predicts droplet attraction and repulsion velocities.

Both the full and simplified models agree at larger droplet separations.

Experimental data confirms the models' validity across different conditions.

Abstract

Recent experiments have shown that liquid drops on highly deformable substrates exhibit mutual interactions. This is similar to the Cheerios effect, the capillary interaction of solid particles at a liquid interface, but now the roles of solid and liquid are reversed. Here we present a dynamical theory for this inverted Cheerios effect, taking into account elasticity, capillarity and the viscoelastic rheology of the substrate. We compute the velocity at which droplets attract, or repel, as a function of their separation. The theory is compared to a simplified model in which the viscoelastic dissipation is treated as a localized force at the contact line. It is found that the two models differ only at small separation between the droplets, and both of them accurately describe experimental observations.

Equations

$$\mathbf{\sigma }(\mathbf{x},t)={\int }_{-\infty }^t\mathrm{d}{t}^{\prime }\Psi (t-t^{\prime })\dot{\mathbf{ϵ}}(\mathbf{x},t^{\prime })-p(\mathbf{x},t)\mathbf{I}.$$

$$\nabla .\mathbf{\sigma }=0,$$

$${{\sigma }_{yy}|}_{y=h_0}=T(x,t)+{\gamma }_s\frac{{\partial }^{2}h}{\partial x^2},$$

$${{\sigma }_{yx}|}_{y=h_0}=0,$$

$$\mathbf{u}(x,0,t)=0.$$

$$\hat{\tilde{h}}(q,\omega )=\hat{\tilde{T}}(q,\omega ){\left[{\gamma }_s{q}^2+\frac{\mu (\omega )}{\tilde{\mathcal{K}}(q)}\right]}^{-1}.$$

$$\mu (\omega )=\mathrm{i}\omega {\int }_0^{\infty }\Psi (t){\mathrm{e}}^{-\mathrm{i}\omega t}\mathrm{d}t=G^{\prime }(\omega )+\mathrm{i}{G}^{\prime \prime }(\omega ).$$

$$\tilde{\mathcal{K}}(q)=\left[\frac{\sinh (2q{h}_0)-2q{h}_0}{\cosh (2q{h}_0)+2(q{h}_0{)}^2+1}\right]\frac{1}{2q}.$$

$$\mu (\omega )=G\left[1+(\mathrm{i}\omega \tau {)}^n\right].$$

$$T(x,t)=\gamma \delta \left(x+\frac{d(t)}{2}\right)+\gamma \delta \left(x-\frac{d(t)}{2}\right).$$

$$h(x,t)=h_1(x,t)+h_2(x,t).$$

$${\tilde{H}}_{-v}(q)=\gamma {\left[{\gamma }_s{q}^2+\frac{\mu (vq)}{\tilde{\mathcal{K}}(q)}\right]}^{-1}.$$

$$h(x,t)=H_{-v}\left(x+\frac{d}{2}\right)+H_v\left(x-\frac{d}{2}\right).$$

$$H_{-v}^{\prime }(d)={\varphi }_v,$$

$$H_{-v}^{\prime }(d)={\int }_{-\infty }^{\infty }\mathrm{i}q{\tilde{H}}_{-v}(q){\mathrm{e}}^{\mathrm{i}qd}\frac{\mathrm{d}q}{2\pi },$$

$$\begin{array}{rl}{\varphi }_v& =-\frac{1}{2}\left[H_v^{\prime }\left(0^{+}\right)+H_v^{\prime }\left(0^{-}\right)\right]\\ & =-{\int }_{-\infty }^{\infty }\Re \left[\mathrm{i}q{\tilde{H}}_v(q)\right]\frac{\mathrm{d}q}{2\pi }.\end{array}$$

$$\begin{array}{rl}& \overline{q}=q{h}_0,\overline{x}=\frac{x}{h_0},\overline{H}=\frac{H}{h_0},\bar{d}=\frac{d}{h_0},\bar{v}=\frac{v\tau }{h_0},\\ & \overline{\omega }=\omega \tau ,\overline{\tilde{\mathcal{K}}}=\frac{\tilde{\mathcal{K}}}{h_0},\overline{\tilde{H}}=\frac{\tilde{H}}{h_0^2},\overline{\mu }=\frac{\mu }{G}.\end{array}$$

$${\tilde{H}}_{-v}(q)=\alpha {\left[{\alpha }_s{q}^2+\frac{\mu (vq)}{\tilde{\mathcal{K}}(q)}\right]}^{-1},$$

$$\alpha =\frac{\gamma }{G{h}_0},{\alpha }_s=\frac{{\gamma }_s}{G{h}_0},$$

$$H_{-v}^{\prime }\simeq \alpha {\int }_{-\infty }^{\infty }\mathrm{i}q{\left[{\alpha }_s{q}^2+\frac{1}{\tilde{\mathcal{K}}(q)}\right]}^{-1}{\mathrm{e}}^{\mathrm{i}qd}\frac{\mathrm{d}q}{2\pi }+\mathcal{O}(v).$$

$$H_{-v}^{\prime }\simeq \alpha C({\alpha }_s){\mathrm{e}}^{-\tilde{\lambda }({\alpha }_s)d},$$

$${\varphi }_v\simeq \alpha \frac{n{2}^{n-1}}{{\alpha }_s^{n+1}\cos \frac{n\pi }{2}}{v}^n.$$

$$v={\left[\frac{{\alpha }_s^{n+1}C({\alpha }_s)}{A(n)}\right]}^{\frac{1}{n}}{\mathrm{e}}^{-\frac{\tilde{\lambda }({\alpha }_s)}{n}d},$$

$$\overline{f}=f/\gamma =H_{v=0}^{\prime }(d).$$

$$\begin{array}{rl}& T_d(x)=\gamma \sin {\theta }_{eq}\left[\delta (x+R)+\delta (x-R)\right]\\ & -\gamma \frac{\sin {\theta }_{eq}}{R}\Theta (R-|x|),\end{array}$$

$${\tilde{T}}_d(q)=2\gamma \sin {\theta }_{eq}\left[\cos (qR)-\frac{\sin (qR)}{qR}\right].$$

$$\tilde{T}(q)={\tilde{T}}_d(q){\mathrm{e}}^{\mathrm{i}q\ell /2}+{\tilde{T}}_d(q){\mathrm{e}}^{-\mathrm{i}q\ell /2},$$

$$\begin{array}{rl}{E}_{e\ell }& =\frac{1}{2}{\int }_{-\infty }^{\infty }T(x)h(x)\mathrm{d}x\\ & =\frac{1}{2G}{\int }_{-\infty }^{\infty }\tilde{T}(q)\tilde{T}(-q)\tilde{K}(q)\frac{\mathrm{d}q}{2\pi }.\end{array}$$

$$E_{int}=\frac{1}{G}{\int }_{-\infty }^{\infty }{\tilde{T}}_d(q){\tilde{T}}_d(-q)\tilde{K}(q)\left({\mathrm{e}}^{\mathrm{i}q\ell }+{\mathrm{e}}^{-\mathrm{i}q\ell }\right)\frac{\mathrm{d}q}{4\pi }.$$

$${\tilde{T}}_d(q)\simeq \frac{(-\mathrm{i}q{)}^2}{2}{Q}_2,$$

$$\begin{array}{rl}{E}_{int}& \simeq \frac{1}{G}{\int }_{-\infty }^{\infty }\frac{(\mathrm{i}q{)}^4}{2^2}{Q}_2^2\tilde{K}(q)\left({\mathrm{e}}^{\mathrm{i}q\ell }+{\mathrm{e}}^{-\mathrm{i}q\ell }\right)\frac{\mathrm{d}q}{4\pi }\\ & =\frac{Q_2^2}{8G}\left(K^{\prime \prime \prime \prime }(\ell )+K^{\prime \prime \prime \prime }(-\ell )\right).\end{array}$$

$$E_{int}=\frac{4}{3\pi }\frac{(\gamma \sin {\theta }_{eq}{)}^2}{G}\frac{R^4}{{\ell }^4}.$$

Full text

I Introduction

The clustering of floating objects at the liquid interface is popularly known as the Cheerios effect Vella and Mahadevan (2005). In the simplest scenario, the weight of a floating particle deforms the liquid interface and the liquid-vapor surface tension prevents it from sinking Vella (2015). A neighboring particle can reduce its gravitational energy by sliding down the interface deformed by the first particle, leading to an attractive interaction. The surface properties of particles can be tuned to change the nature of interaction, but two identical spherical particles always attract Kralchevsky and Nagayama (2000). Anisotropy in shape of the particles or curvature of the liquid interface adds further richness to this everyday phenomenon Loudet et al. (2005); Botto et al. (2012); Ershov et al. (2013). Self-assembly of elongated mosquito eggs on the water surface provides an example of this capillary interaction in nature Loudet and Pouligny (2011), while scientists have exploited the effect to control self-assembly and patterning at the microscale Bowden et al. (1997); Lewandowski et al. (2010); Furst (2011); Cavallaro et al. (2013); Stamou et al. (2000).

The concept of deformation-mediated interactions can be extended from liquid interfaces to highly deformable solid surfaces. Similar to the Cheerios effect, the weight of solid particles on a soft gel create a depression of the substrate, leading to an attractive interaction Chakrabarti et al. (2015); Chakrabarti and Chaudhury (2014, 2013). Recently, it was shown that the roles of solid and liquid can even be completely reversed: liquid drops on soft gels were found to exhibit a long-ranged interaction, a phenomenon called the inverted Cheerios effect Karpitschka et al. (2016). An example of such interacting drops is shown in fig. 1(a). In this case, the droplets slide downwards along a thin, deformable substrate (much thinner as compared to drop size) under the influence of gravity, but their trajectories are clearly deflected due to a repulsive interaction between the drops. Here capillary traction of the liquid drops instead of their weight deforms the underlying substrate (cf fig. 1(b)). The scale of the deformation is given by the ratio of liquid surface tension to solid shear modulus ($\gamma/G$), usually called the elasto-capillary length. Surprisingly, drops on a thick polydimethylsiloxane (PDMS) substrate (much thicker than the drop sizes) were found to always attract and coalesce, whereas for a thin substrate (much thinner than the drop sizes) the drop-drop interaction was found to be repulsive Karpitschka et al. (2016). This interaction has been interpreted as resulting from the local slope of the deformation created at a distance by one drop, which can indeed be tuned upon varying the substrate thickness.

Here we wish to focus on the dynamical aspects of the substrate-mediated interaction. Both the Cheerios effect and the inverted Cheerios effect are commonly quantified by an effective potential, or equivalently by a relation between the interaction force and the particle separation distance Vella and Mahadevan (2005); Karpitschka et al. (2016). However, the most direct manifestation of the interaction is the motion of the particles, moving towards or away from each other. In the example of fig. 1(a), the dark gray lines represent drop trajectories and the arrows represent their instantaneous velocities. When drops are far away, the motion is purely vertical due to gravity with a steady velocity $\vec{v}_{g}$. As the two drops start to interact, a velocity component along the line joining the drop centers develops ($\vec{v}_{i}-\vec{v}_{ig}$, $i=1,2$) that move the drops either away or towards each other. In fig. 1(c), we quantify the interaction by $\Delta(\vec{v}-\vec{v}_{g})=(\vec{v}_{1}-\vec{v}_{1g})-(\vec{v}_{2}-\vec{v}_{2g})$, representing the relative interaction velocity as a function of inter drop distance ($d$) for several interaction events. Even though the dominant behavior is repulsive, we find an attractive regime at small distances. Similar to the drag force on a moving particle at a liquid interface Dörr et al. (2016), viscoelastic properties of the solid resist the motion of liquid drops on it Carre et al. (1996). Experiments performed in the over damped regime allow one to extract the interaction force from the trajectories, after proper calibration of the relation for the drag force as a function of velocity Karpitschka et al. (2016).

In this paper, we present a dynamical theory of elastocapillary interaction of liquid drops on a soft substrate. For large drops, the problem is effectively two dimensional (cf. fig. 1(b)) and boils down to an interaction between two adjacent contact lines. By solving the dynamical deformation of the substrate, we directly compute the velocity-distance curve based on substrate rheology, quantifying at the same time the interaction mechanism and the induced dynamics. In sec. II we set up the formulation of the problem and propose a framework to introduce the viscoelastic properties of the solid. In sec. III we present our main findings in the form of velocity-distance plots for the case of very large drops. These results are compared to the previous theory for force-distance curves, and we investigate the detailed structure of the dynamical three phase contact line. Subsequently, our results are generalised in sec. IV to the case of finite-sized droplets, and a quantitative comparison with experiments is given.

II Dynamical Elastocapillary interaction

II.1 Qualitative description of the interaction mechanism

Before developing a detailed dynamical theory of the inverted cheerios effect, we start by a qualitative discussion of the interaction mechanism. The problem of drop-drop interaction on soft substrates can be characterised by four length scales: the radius of the drop ($R$), the thickness of the substrate ($h_{0}$), the elastocapillary length ($\gamma/G$), and of course the gap separating the drops ($d$). In the case of drop size much larger than all other length scales: $R\gg\textup{max}(h_{0},d,\gamma/G)$, one can neglect the curvature of the contact line and the wetting ridge becomes quasi two-dimensional as in fig. 1(b). Furthermore the drop-drop interaction is limited to an interaction between the neighboring contact lines. This provides an important simplification: the liquid contact angles at the outer contact lines of each drop are essentially unaffected by the presence of the second drop, and the slow relaxation of the viscoelastic solid ensures circular liquid - vapor interface. As a result, the liquid contact angles at the interacting contact lines also remain constant to satisfy the circular cap shape of the liquid-vapor interface. In the example of fig. 1(b) this liquid contact angle with respect to horizontal ($\theta_{\ell}$) is $\pi/2$ by symmetry, since we assume $\gamma_{sv}=\gamma_{s\ell}=\gamma_{s}$.

To understand the interaction mechanism, we magnify the deformation around neighbouring contact lines in fig. 2(a). It shows a sketch of two repelling contact lines with velocity $v$. Our goal is to find the value of $v$ for a given distance $d$. Deformation of the substrate at the three phase contact line is fully defined by the balance of three respective surface tensions, known as the Neumann condition Marchand et al. (2012); Lubbers et al. (2014); Park et al. (2014). The assumption of $\gamma_{sv}=\gamma_{s\ell}=\gamma_{s}$ leads to the symmetric wetting ridges of fig. 2(a). As demonstrated in the next section, in first approximation the solid profile is obtained by the superposition of the deformations due to the two moving contact lines shown in fig. 2b. However, motion on a viscoelastic medium causes a rotation of the wetting ridge by an angle $\phi_{v}$ that depends on the velocity: the larger the velocity, the larger the rotation angle Karpitschka et al. (2015). At the same time, superposition results in a rotation of a wetting ridge due to the local slope of the other deformation profile ($H^{\prime}_{-v}(d)$) at that location. Hence, the force balance with the capillary traction from an unrotated liquid interface can be maintained when both rotations exactly cancel each other. Therefore, by equating $H^{\prime}_{-v}(d)$ to $\phi_{v}$, we obtain a relation between velocity and distance.

II.2 Formulation

In this section, we calculate the dynamic deformation of a linear, viscoelastic substrate below a moving contact line. This step will enable us to quantify the interaction velocity, according to the procedure described above.

Due to the time-dependent relaxation behaviour of viscoelastic materials, the response of the substrate requires a finite time to adapt to any change in the imposed traction at the free surface. Under the assumption of linear response, the stress-strain relation of a viscoelastic solid is given by Ferry (1961)

[TABLE]

Here $\bm{\sigma}$ is the stress tensor, $\bm{\epsilon}$ is the strain tensor, $\Psi$ is the relaxation function, $p$ is the isotropic part of the stress tensor, $\bm{I}$ is the identity tensor, and $\bm{x}=\{x,y\}$. Throughout we will assume plain strain conditions, for which all strain components normal to the page vanish. We are interested in the surface profile of the deformed substrate, given by the vertical ($y$)-component of displacement vector $\bm{u}$. For small deformations this imply $u_{y}(x,h_{0},t)=h(x,t)$. The soft substrate is attached to a rigid support at its bottom surface and subjected to traction $T(x,t)$ on the top surface due to the pulling of the contact line. In addition, the surface tension of the solid $\gamma_{s}$ provides a solid Laplace pressure that acts as an additional traction on the substrate. For inertia-free dynamics, the Cauchy equation along with the boundary conditions for the plane strain problem are given by

[TABLE]

Following previous works Long et al. (1996); Karpitschka et al. (2015), we solve eqns. (1) and (2) using a Green’s function approach. Applying Fourier transforms in both space and time ($x\rightarrow q$ noted by ‘~’, $t\rightarrow\omega$ noted by ‘^’), the deformation can be written as,

[TABLE]

Forward and backward Fourier transforms (with respect to space) of any function is defined as, $\widetilde{f}(q)=\int_{-\infty}^{\infty}f(x)\mathrm{e}^{-\mathrm{i}qx}\mathrm{d}x$ and $f(x)=\int_{-\infty}^{\infty}\widetilde{f}(q)\mathrm{e}^{\mathrm{i}qx}\mathrm{d}q/2\pi$. Transformation rules are equivalent for $t\leftrightarrow\omega$. Here $\mu(\omega)$ is the complex shear modulus of the material and consists of the storage modulus $G^{\prime}(\omega)=\Re[\mu(\omega)]$, and loss modulus $G^{\prime\prime}(\omega)=\Im[\mu(\omega)]$:

[TABLE]

The spatial dependence is governed by the kernel $\widetilde{\mathcal{K}}(q)$, which for an incompressible, elastic layer of thickness $h_{0}$ reads

[TABLE]

Two backward transforms of (3) give the surface deformation $h(x,t)$.

In what follows, we perform the calculations assuming a power-law rheology of the soft substrate, characterised by a complex modulus,

[TABLE]

This rheology provides an excellent description of the reticulated polymers like PDMS, as are used in experiments on soft wetting Jerison et al. (2011); Style et al. (2013a); Karpitschka et al. (2015); Style et al. (2013b). In particular, the loss modulus has a simple power-law form, ${G}^{\prime\prime}\sim G(\omega\tau)^{n}$, where $G$ is the static shear modulus and $\tau$ is the relaxation timescale. The exponent $n$ depends on the stoichiometry of the gel, and varies between 0 and 1 Winter and Chambon (1986); Scanlan and Winter (1991); Chambon and Winter (1987). The experimental data presented in fig. 1 is for a PDMS gel with a measured exponent of $n=0.61$ and timescale $\tau=0.68s$.

II.3 Two steadily moving contact lines

Now we turn to two moving contact lines, as depicted in fig. 2(a). Laplace pressure of the liquid drops scales with $R^{-1}$, and vanish in this limit of large drops. Hence, the applied traction is described by two delta functions (noted $\delta$), separated by a distance $d(t)$:

[TABLE]

Owing to the linearity of the governing equations, the deformations form a simple superposition

[TABLE]

Here $h_{1}(x,t)$ and $h_{2}(x,t)$ respectively are the solutions obtained for contact lines moving to the left and right.

Before considering the superposition, let us first discuss the deformation due to a contact line centered at the origin and moving steadily with velocity $-v$. The traction and deformation are then of the form $\gamma\delta(x+vt)$, and $H_{-v}(x+vt)$. When computing the shape in a co-moving frame, it can be evaluated explicitly from (3) as

[TABLE]

The effect of velocity is encoded in the argument of the complex modulus $\mu$, and hence couples to the dissipation inside the viscoelastic solid. Setting $v=0$, this expression (9) gives a perfectly left-right symmetric deformation of the contact line. As argued, for $v\neq 0$, the motion breaks the left-right symmetry and induces a rotation of the wetting ridge.

Now, for two contact lines located at $\pm d(t)/2$ moving quasi-steadily with opposite velocities $v=\pm\frac{1}{2}\dot{d}$, the deformation in eq. (8) can be written as

[TABLE]

Under this assumption of quasi-steady motion, there is no explicit time-dependence; the shape is fully known once the velocity $v$ and the distance $d$ are specified. However, as discussed in sec. IIA, the total deformation of eq. (10) should be consistent with a liquid angle of $\pi/2$. This is satisfied when the ridge rotation due to motion, $\phi_{v}$, perfectly balances the slope induced by the second drop. Summarizing this condition as

[TABLE]

we can compute the velocity $v$ for given separation $d$ of the two contact lines. Explicit expressions can now be found from equation (9), as

[TABLE]

and (See ‘Methods’ section of Karpitschka et al. (2015) for details)

[TABLE]

In what follows we numerically solve eq. (11), and obtain the interaction velocity as a function of $d$.

II.4 Dimensionless form

The results below will be presented in dimensionless form. We use the substrate thickness $h_{0}$ as the characteristic lengthscale, relaxation timescale of the substrate $\tau$ as the characteristic timescale, and introduce the following dimensionless variables:

[TABLE]

In the remainder of the paper we drop the overbars, unless specified otherwise. Eqns (11), (12), (13) remain unaltered in dimensionless form, but the scaled deformation reads

[TABLE]

This form suggests that the two main dimensionless parameters of the problem are,

[TABLE]

which can be seen as two independent elastocapillary numbers. $\alpha$ compares the typical substrate deformation to its thickness, and as discussed in the next section $\alpha_{s}$ measures the decay of substrate deformation compared to $h_{0}$. The assumption of linear (visco)elasticity demands that the typical slope of the deformation is small, which requires $\alpha/\alpha_{s}=\gamma/\gamma_{s}\ll 1$. In the final section we will also consider drops of finite size $R$. In that case, the ratio $R/h_{0}$ appears as a third dimensionless parameter.

III Results

III.1 Relation between velocity and distance

Here we discuss our main result, namely the contact line velocity $v$ of the inverted Cheerios effect, as a function of distance $d$. Figure 1(d) represents a typical velocity-distance plot for $\alpha_{s}=1$ and $n=0.7$ obtained by numerically solving eq. (11). The theory shows qualitative agreement with the experiments: the interaction velocity is repulsive at large distance, but displays an attraction at small distances. The change in the nature of interaction can be understood from the local slope of deformation. As can be seen from fig. 2(b), the slope $H^{\prime}_{-v}$ is negative at small distances from the contact line. This induces a negative velocity (attraction) to the contact line located at $d/2$. However, incompressibility of the material causes a dimple around the wetting ridge, and at $d\sim 1$ the slope changes its sign and so does the interaction velocity. A fully quantitative comparison between theory and experiment, taking into account finite drop size effects, will be presented at the end of the paper in Sec. V.

We now study the velocity-distance relation in more detail. In fig. 3(a), we show the interaction velocity on a semi-logarithmic scale, for different values of the rheological parameter $n$. In all cases the interaction is attractive at small distance, and repulsive at large distance. The point where the interaction changes sign is identical for all these cases. The straight line at large $d$ reveals that the interaction speed decays exponentially, $v\sim\mathrm{e}^{-\lambda d}$. The exponent $\lambda$ is not universal, but varies as $\lambda\sim 1/n$. Indeed by rescaling the velocity as $|v|^{n}$, we observed a collapse of the curves (cf. fig. 3(b)) in the ranges where the interaction velocities are small. The inset of fig. 3(b), however, shows that the collapse does not hold at very small distances where the velocities are larger.

To understand these features, we extract the far-field asymptotics of the velocity-distance curve. As $v$ decreases for large $d$, we first expand eq. (12) for small $v$ to get

[TABLE]

This effectively corresponds to the slope of a stationary contact line and the viscoelastic effects drop out. Hence, at large distances away from the contact line the deformed shape is essentially static. The substrate elastocapillary number ($\alpha_{s}$) governs the deformed shape in far field. In the Appendix we evaluate this static slope, which yields

[TABLE]

indeed providing an exponential decay. This relation confirms that $\alpha_{s}$ governs the decay of substrate deformation. Similarly, we expand $\phi_{v}$ in eq. (13) for small $v$ as (see Karpitschka et al. (2015)),

[TABLE]

Hence, the rotation of the ridge inherits the exponent $n$ from the rheology. Comparison of eqs (18) and (19) gives the desired asymptotic relation

[TABLE]

with $A(n)=n2^{n-1}/\cos(n\pi/2)$. For $\alpha_{s}=1$, $\widetilde{\lambda}=0.783$ and $C=0.069$. The red dashed line in fig. 3(a) confirms this asymptotic relation. The quantity $A(n)v^{n}$ is independent of the substrate rheology and collapses to a universal curve for large distances. Only at small separation, a weak dependency on $n$ appears (see the inset of fig. 3(b)).

III.2 Force-distance

We now interpret these velocity-distance results in terms of an interaction force, mediated by the substrate deformation. The underlying physical picture is that the interaction force between the two contact lines induces a motion to the drops. In a regime of overdamped dynamics, dissipation within the liquid drop can be neglected, and the interaction force must be balanced by a dissipative force due to the total viscoelastic drag inside the substrate. In fact, eq. (11) can be interpreted as a force balance: the slope $H^{\prime}_{-v}(d)$ induced by the other drop is the interaction force, while the dissipative force is given by the rotation $\phi_{v}$ due to drop motion. Viscoelastic dissipation given by the relation $G^{\prime\prime}(\omega)\sim\omega^{n}$ manifests itself in the ridge rotation through, $\phi_{v}\sim v^{n}$. The blue line in fig. 4, shows the corresponding dissipative force as a function of distance for $n=0.7$.

In Ref. Karpitschka et al. (2016) we presented a quasi-static approach to compute the interaction force. There, the main hypothesis was that the dissipative force can be modelled as a point force $f$ that is perfectly localized at the contact line. We are now in the position to test this hypothesis by comparing its predictions to those of the dynamical theory developed above, where the viscoelastic dissipation inside the solid is treated properly. In the localized-dissipation model, the deformation profiles are treated as perfectly static and thus involve $H^{\prime}_{v=0}(d)$. However, in this approximation the change in deformation profile now leads to a local violation of the Neumann balance (cf. inset figure 5(a)), since the static model has no ridge rotation, $\phi_{v=0}=0$. This imbalance in the surface tensions must be compensated by the localized dissipative force $f$, which can subsequently be equated to the interaction force. Hence, for the case of large drops, the localized-force model predicts

[TABLE]

Like before, we drop the overbar and denote the dimensionless force by $f$. The result of the localized-force model is shown as the red data points in fig. 4, showing $f$ versus $d$. At large distance this gives a perfect agreement with the dynamical theory. For large $d$, $H^{\prime}_{-v}(d)\simeq H^{\prime}_{v=0}(d)$, and the dissipative force $\phi_{v}$ perfectly agrees with the interaction force $f$. In this limit, slope of the deformation decays exponentially with distance, and so is the interaction force. However, as $d\rightarrow 0$, the static profiles differ from the dynamic profiles, inducing a small error in the calculation of the interaction force.

The differences between the localized-dissipation model and the full viscoelastic dynamic calculation can be interpreted in terms of true versus apparent contact angles. Figure 5(a) shows the profile in the approximation of the localized-force model. The inset shows the violation of the Neumann balance: the red dashed line is the dynamic profile obtained from our viscoelastic theory. A more detailed view on the contact angles is given in fig. 5(b,c). Here we compare the dynamical slope (red dashed line) to that of the localized-force model (gray solid line), on both sides of the contact line at $x=d/2$. The two profiles have the same slope beyond $x\sim\gamma_{s}/Gh_{0}$. However, differences appear at small distance near the contact line where viscoelastic dissipation plays a role. While the dynamic profiles are left-right symmetric, the localized-force model gives a symmetry breaking: the difference in slope between the two approach at the contact line gives the magnitude of $f$.

In summary, the inverted Cheerios effect can be characterised by an interaction force, which at separations ($d$) larger than $\gamma_{s}/Gh_{0}$ is independent of substrate rheology. In this regime the interaction force can be computed assuming that all dissipation is localized at the contact line – this gives rise to an apparent violation of the Neumann balance, due to a dissipative force $f$. In reality the Neumann balance is restored at small distance, but the total dissipative force is indeed equal to $f$. At separations $d<\gamma_{s}/Gh_{0}$, however, the spatial distribution of dissipation has significant effects and the localized-force model gives rise to some (minor) quantitative errors.

IV Drops of finite size

We now show that interaction is dramatically changed when considering drops of finite size $R$, compared to the substrate thickness $h_{0}$. A fully consistent dynamical theory is challenging for finite-sized drops, since the two contact lines of a single drop are not expected to move at the same velocity. However, the previous section has shown that the localized-force model provides an excellent description of the interaction force, which becomes exact at separations $d>\gamma_{s}/Gh_{0}$. We therefore follow this approximation to quantify the inverted Cheerios effect for finite sized drops.

IV.1 Effect of substrate thickness

We assume that during interaction, radius of both the drops does not change and the liquid cap remains circular (i.e. no significant dissipation occurs inside the liquid). The traction applied by a single drop on the substrate comprised of two localized forces pulling up at the contact lines and a uniform Laplace pressure pushing down inside the drop. The traction applied by a single drop of size $R$, centered at the origin, is thus given by

[TABLE]

expressed in dimensional form. Here $\Theta$ is the Heaviside step function. The corresponding Fourier transform reads

[TABLE]

Considering two drops centered at $x=-\ell/2$ and $x=\ell/2$ respectively, with corresponding tractions of $T_{d}(x+\ell/2)$ and $T_{d}(x-\ell/2)$, the total traction induced by the two drops reads, in Fourier space

[TABLE]

where the inter-drop separation $d=\ell-2R$. Using the static Green’s function approach, we compute the corresponding deformation profiles. A typical profile obtained from this traction is shown in fig. 6.

The problem is closed by determining the liquid contact angles at both the contact lines ($\theta_{i}$, $\theta_{o}$) self-consistently with the elastic deformation Lubbers et al. (2014). The interaction is found from reestablishing a violated Neumann balance by a dissipative force (for details we refer to Karpitschka et al. (2016)). As such we determine the force-distance relation for varying drop sizes. The results are shown in fig. 6(b), for various drop sizes $h_{0}/R$. Clearly, upon decreasing the drop size, or equivalently increasing the layer thickness, the interaction force becomes predominantly attractive. For a thin substrate, incompressibility leads to a sign change of the interaction force. As the substrate thickness is increased the profile minimum is moved away from the contact line. For an elastic half space, i.e. $h_{0}/R=\infty$, the deformation monotonically decays away from the contact line giving rise to a purely attractive interaction. The inset of fig. 6(b) shows the interaction force for an elastic half space.

IV.2 Far-field asymptotics on a thick substrate

When two drops are far away from each other, $\theta_{i}=\theta_{o}\simeq\theta_{eq}$ and the liquid cap resembles its equilibrium shape. In this limit, the drop-drop interaction is mediated purely by the elastic deformation. Here we calculate the large-$d$ asymptotics of the interaction force. We follow a procedure similar to that proposed for elasticity-mediated interactions between cells Schwarz and Safran (2002); Bischofs et al. (2004), where the cells are treated as external tractions. The total energy of the system $E_{tot}$ then consists of the elastic strain energy $E_{e\ell}$ plus a term $E_{w}$ originating from the work done by the external traction. At fixed separation distance, the mechanical equilibrium of the free surface turns out to give $E_{e\ell}=-\frac{1}{2}E_{w}$ (the factor 1/2 being a consequence of the energy in linear elasticity). So, we can write $E_{tot}=E_{e\ell}+E_{w}=-E_{e\ell}$ Bischofs et al. (2004). Subsequently, the interaction force between the droplets is computed as $f=-\frac{\mathrm{d}E_{tot}}{\mathrm{d}d}=\frac{\mathrm{d}E_{e\ell}}{\mathrm{d}d}$. The task will thus be to evaluate $E_{e\ell}$ from the traction induced by two drops.

In the absence of body forces, the elastic energy (per unit length) is written as

[TABLE]

Plugging in the traction applied by two drops given in eq. (24) into eq. (25) we find that there are two parts of the elastic energy: one is due to the work done by individual tractions (the self-energy of individual drops), and an interaction term ($E_{int}$) given by

[TABLE]

To obtain the far-field interaction energy ($\ell\gg R$), we now perform a multipole expansion of the traction $T_{d}(x)$, or equivalently a long-wave expansion of its Fourier transform. The quadrupolar traction of eq. (23) gives, after expansion

[TABLE]

were $Q_{2}$ is the second moment of the traction and is given by, $Q_{2}=\int_{-\infty}^{\infty}T_{d}(x)x^{2}\mathrm{d}x={4\over 3}\gamma\sin\theta_{eq}R^{2}$. So, the interaction energy in the far-field reads

[TABLE]

For drops on a finite thickness, the Green’s function $K(x)$ decays exponentially and again gives an exponential cut-off of the interaction at distances beyond $h_{0}$.

A more interesting result is obtained for an infinitely thick substrate, formally corresponding to $h_{0}\gg\ell$. The kernel for a half space is $K(x)=-{\log|x|\over 2\pi}$, which can e.g. be inferred from the large thickness limit of (5). Consequently, $K^{\prime\prime\prime\prime}(\ell)=K^{\prime\prime\prime\prime}(\ell)={3\over\pi\ell^{4}}$, so that the interaction energy becomes

[TABLE]

This is identical to the interaction energy between capillary quadrupoles Stamou et al. (2000); Danov et al. (2005); Danov and Kralchevsky (2010), except for the prefactor that depends on elastocapillary length in the present case. We calculate the force of interaction by $f=\frac{\mathrm{d}E_{int}}{\mathrm{d}\ell}$, giving

[TABLE]

where the result is collected in dimensionless groups. Here we replaced $\ell\simeq d$, which is valid when the drop separation is large compared to $R$. The red dashed line in the inset of fig. 6(b) confirms the validity of the above asymptotic result, including the prefactor. Hence, in the limit of large thickness the attractive interaction force decays algebraically, as $1/d^{5}$, which is in stark contrast to the exponential decay observed for thin substrates.

It is instructive to compare eq. (29) with the interaction energy for capillary quadrupoles formed by rigid anisotropic particles embedded in a fluid-fluid interface Stamou et al. (2000); Botto et al. (2012). In that case an undulated triple line is formed owing to contact line pinning Stamou et al. (2000) or to the change in curvature of the solid surface along the triple line Loudet et al. (2006); Lewandowski et al. (2010). For capillary quadrupoles the far-field interaction energy is proportional to $\gamma H_{p}^{2}(R/\ell)^{4}$ , where $H_{p}$ is the amplitude of the quadrupolar contact line distortion. Despite an identical power-law dependence, the prefactors in eq. 29 and in the corresponding expression for capillary quadrupoles are different, reflecting a different physical origin. A capillary quadrupole is formed by four lobes located at an average distance R from the quadrupole center: two positive meniscus displacements of amplitude $H_{p}$ along one symmetry axis of the particle, and two negative displacements along the perpendicular axis. The presence of a second particle changes the slope of equal sign meniscii along the center-to-center line, resulting in a net capillary attraction. The corresponding interaction force per unit length of contact line is proportional to $\gamma(H_{p}/R)^{2}(R/d)^{5}$, to be compared with eq. 30. The non-dimensional ratio $(H_{p}/R)$ is a geometric feature of the particle that depends primarily on the particle shape and contact angle Botto et al. (2012). This ratio is independent of the particle size. On the contrary, the non-dimensional ratio $\frac{\gamma}{GR}$ does depend on the particle size, as well as on the material-dependent scale $\gamma/G$. A differentiating feature is also that in the elastocapillary quadrupole the upward and downward distortions of the elastic surface are located along the line of interaction, thus forming a linear quadrupole. The identical power-law dependence for capillary and elasto-capillary quadrupoles is incidental. It is a consequence of the fact that the $\log$ function is the Green’s function for both the small-slope Young-Laplace equation of capillarity and the equation governing 2D linear elasticity.

V Discussion

Finally, let us turn to experiments of the inverted Cheerios effect. The results in fig. 1 already show a qualitative agreement between theory and experiment. A quantitative comparison in terms of “velocity vs distance” is challenging, since the experimental motion follows is due to a combination of drop-drop interaction and gravity – the latter was not taken into account in the dynamical theory. However, gravity is easily substracted when representing the result as “force vs distance”, as was also done in Karpitschka et al. (2016). The experimental data are reproduced in fig 7: On thick substrates the interaction is purely attractive (panel b), while a change in sign of the interaction force is observed on thin substrates (panel a). The two-dimensional model presented in this paper predicts a force per unit length $f$. Hence, to compare to experiment one needs to multiply the result with a proper length over which this interaction is dominant. The blue lines are obtained by choosing a length of $\pi R/6$, i.e. a small fraction of the perimeter, giving a good agreement with experiment. In Karpitschka et al. (2016) we provided numerically a three-dimensional result, based on the same localized-force approach. The result is shown as the red lines in fig. 7 and is similar to the two-dimensional case. This good agreement with experiment shows that the elasto-capillary theory presented here provides an accurate description of the physics. The similarity between the two- and three-dimensional models shows that the interaction is dominated by the regions where the contact lines are closest.

In conclusion, we present a theoretical model to capture the dynamical interaction between two liquid drops on a soft, solid surface. For large drops the interaction is limited to neighboring contact lines, and we find the interaction velocity-distance relation in this limit. For large distances, the velocity decays exponentially with an exponent that depends both on substrate rheology and surface tension. The over damped dynamics of drop motion allow us to relate the velocity to an effective viscoelastic drag force at the contact line. We also connect our theory to a quasi static description of the interaction, quantified in terms of force-distance curves, and find that the two approaches agree well for large inter drop distances. Motivated by this, we extend the quasi-static force calculation to finite sized drops to show how the ratio of substrate thickness to drop radius governs the nature of the interaction. In the other physically relevant limit of small drops on a very thick substrate, the interaction force is found to decay algebraically for large distances. We anticipate that the dissipation mediated interaction mechanism studied here to play a key role in collective dynamics of drops on soft surfaces Style et al. (2013a); Bueno et al. (2017), and dynamics of soft, adhesive contacts.

Appendix A

We evaluate the complex counterpart of the integrand in eq. (17) on a contour $C$, given by

[TABLE]

Here $\mathbf{q}$ is the complex wave vector. The contour $C$ consists of a semicircle of radius $r$ centered at the origin and a straight line connecting the ends of the semicircle. We evaluate the above integral using Residue theorem. As $r\rightarrow\infty$, the integrand vanishes on the semicircle and the contour integral simplifies to eq. (17). For any $\alpha_{s}$ the integrand of eq. 31 has an infinite number of poles on the complex plane. The leading order behavior of the integration is governed by the pole with smallest imaginary part ($\mathbf{q}_{pole}$) as other poles lead to much faster decay.

The substrate elastocapillary number $\alpha_{s}$ dictates the value of $\mathbf{q}_{pole}$. Figure 8 shows that below a critical $\alpha_{s}$, $\mathbf{q}_{pole}$ is purely imaginary suggesting an purely exponential decay of the deformation. Beyond the critical $\alpha_{s}$, $\mathbf{q}_{pole}$ has a nonzero real and imaginary part giving rise to oscillatory behavior of slope as predicted by Long et. al Long et al. (1996). Hence, the large distance approximation of slope is given by,

[TABLE]

This expression simplifies to eq. (18) in the main text.

Acknowledgements.

L.B. acknowlwdges financial support from European Union Grant CIG 618335. S.K. acknowledges financial support from NWO through VIDI Grant No. 11304. A.P. and J.H.S. acknowledge financial support from ERC (the European Research Council) Consolidator Grant No. 616918.