<figure><img src="image_1.png" /></figure>

HAL
open science

# Drops Can Perpetually Bounce over a Vibrating Wettable Solid

Lebo Molefe, Tomas Fullana, François Gallaire, John Kolinski

► To cite this version:

Lebo Molefe, Tomas Fullana, François Gallaire, John Kolinski. Drops Can Perpetually Bounce over a Vibrating Wettable Solid. Physical Review Letters, 2025, 135 (14), pp.144001. 10.1103/w3qq-cnj3. hal-05314880

HAL Id: hal-05314880

https://hal.science/hal-05314880v1

Submitted on 14 Oct 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
Drops Can Perpetually Bounce over a Vibrating Wettable Solid

Lebo Molefe<sup>1,2,*</sup> Tomas Fullana<sup>1,2,*†</sup> François Gallaire<sup>1,2</sup> and John M. Kolinski<sup>1</sup>

<sup>1</sup>Engineering Mechanics of Soft Interfaces (EMSI) Laboratory, EPFL, Lausanne, Switzerland

<sup>2</sup>Laboratory of Fluid Mechanics and Instabilities (LFMI), EPFL, Lausanne, Switzerland

<figure><img src="image_3.png" /></figure>

(Received 7 May 2025; accepted 12 August 2025; published 2 October 2025)

On the surface of a vibrating liquid bath, instead of coalescing, a drop will continually bounce on a thin film of air between the drop and the free surface, giving rise to rich chaotic dynamics and quantum analog behavior. However, perpetual bouncing is yet to be demonstrated on a vibrating rigid solid, where the control of the drop's motion is not limited by the bath dynamics. Here we show that vibration of an atomically smooth mica surface prolongs a drop's hovering state by several orders of magnitude, increasing the time to wet from less than a second to several minutes. The excitation of the second spherical harmonic mode determines a transition between a bouncing state with high-amplitude rebounds, and a bound state, where the drop's motion is locked onto the vibrating solid above a thin air layer. We further develop a coupled linear spring model, accounting for the drop's deformation, that predicts bouncing trajectories without fitting parameters. Our results extend the scope of bouncing drop studies beyond the liquid bath and highlight the potential for precision manipulation of small liquid quantities in air.

DOI: 10.1103/w3qq-cnj3

*Introduction*—Reducing liquid-solid contact facilitates the incredible performance of solid surfaces with antifouling [1,2], anti-icing [3,4], antibacterial [5], and self-cleaning [6,7] properties. The Leidenfrost effect [8,9], where a drop levitates on its own vapor film above a hot surface for an extended period, eliminates wetting completely, but heat could damage sensitive chemistry in the drop, and leads to evaporation of the entire liquid mass. A stable hovering state for a liquid drop on a solid surface at room temperature is yet to be demonstrated.

In this Letter, we experimentally show that surface vibration can extend a hovering state for a drop on a smooth solid [10–13] by several orders of magnitude, demonstrating a purely kinetic analog of the Leidenfrost effect. Drops hover near the solid in a no-contact state, with their bouncing amplitude controlled by the surface vibration amplitude and frequency. Paradoxical gravity-defying states driven by vibration [14–16] have been observed for drops on vibrating liquid baths, where the bath's free surface deformation drives the drop's rich dynamics [17–27]. Here, because the solid is rigid, the drop's resonant behavior is governed by its deformation; we find that a simple reduced-order model that couples the drop's center

of mass and its second harmonic mode deformation is sufficient to predict drop trajectories without fitting parameters.

For a drop impacting a stationary solid surface, the dynamics are governed by the contributions of inertia, gravity, and viscosity relative to surface tension. The three dimensionless groups are the Weber number $We = \rho r_0 V^2 / \gamma$, the Bond number $Bo = \rho g r_0^2 / \gamma$, and the Ohnesorge number $Oh = \eta (\gamma \rho r_0)^{-1/2}$ ($\rho$: liquid density, $\gamma$: surface tension, $\eta$: dynamic viscosity, $r_0$: drop radius, $V$: impact velocity, and $g$: gravitational acceleration). In our experiments, surface tension is a significant restoring force [10,28], which allows bouncing without contact on initial impact due to the presence of a thin compressed air film between the drop and the surface [11–13,29–31]. In the case of a vibrating surface with position given by $z_s(t) = \frac{1}{2} A_s \sin(\omega_s t + \phi_s)$ [Fig. 1(a)], the vibration frequency $f_s = \omega_s / (2\pi)$ and amplitude $A_s$ control the drop's rebound dynamics via two additional parameters: the vibration number $\Omega = \omega_s / \omega_d$ and the dimensionless acceleration $\Gamma = \frac{1}{2} A_s \omega_s^2 / g$, where $\omega_d = (\gamma / \rho r_0^3)^{1/2}$ is the drop's Rayleigh frequency [32].

*Generation of bouncing drops*—We use a high-speed camera to measure the dynamics of bouncing silicone oil drops with $\rho = 941 \pm 9$ kg/m$^3$, $\eta = 20.17 \pm 0.07$ mPa.s, $\gamma = 21.2 \pm 0.2$ mN/m, and $r_0 = 0.78 \pm 0.05$ mm. A smooth mica sheet, with roughness on the order of a few nanometers (see Supplemental Material [33]), is attached to the vibrating stage and the drop is released from a syringe located between $z_0 = 2–4$ mm above the surface depicted in Fig. 1(a), off of which it will bounce

<footnote>*These authors contributed equally to this work.</footnote>

<footnote>†Contact author: tomas.fullana@epfl.ch</footnote>

<footnote>Published by the American Physical Society under the terms of the Creative Commons Attribution 4.0 International license. Further distribution of this work must maintain attribution to the author(s) and the published article's title, journal citation, and DOI.</footnote>
<figure><img src="image_4.png" /><figcaption>FIG. 1. Drop bouncing on air above a vibrating smooth solid. (a) Schematic of the experimental setup. (b) Still experimental images with dashed lines indicating the surface position and arrows indicating its direction of motion. (c) Corresponding snapshots from the DNS. (d),(e) Kymographs from a cut through the drop center with the stage position indicated by a solid line; here $r_0 = 0.76$ and $A_s = 0.21$ mm. (d) Bouncing case with $f_s = 40$ Hz ($\Omega \approx 1.2$ and $\Gamma \approx 0.7$). (e) Bound case with $f_s = 180$ Hz ($\Omega \approx 5.2$ and $\Gamma \approx 13.7$). (f),(g) Drop center-of-mass position $z(t)$, stage position $z_s(t)$, and air layer thickness for both cases. Insets show image sequences of a cycle, with $\Delta t = 0.014$ (bouncing case) and $\Delta t = 0.0026$ s (bound case). (h),(i) Variation of drop's width $\Delta w$ and height $\Delta h$ as a percentage of the drop diameter. (f)–(i) DNS matching the drops' trajectories and deformations in dash-dotted lines. (j),(k) Drop energies $\mathcal{E}$ normalized by the initial total energy $T_0$ in DNS for the bouncing and bound cases.</figcaption></figure>

[Fig. 1(b)]. The corresponding dimensionless numbers are $We \approx 1.25-2.95$ with $V \approx 0.20-0.28$ m/s, $Bo \approx 0.26$, and $Oh \approx 0.16$, which indicate that the drop will initially bounce off the stage even in the absence of surface vibration [28]. The stage position $z_s(t)$, driven with a function generator and a mechanical shaker, has amplitude $A_s = 0.15-0.25$ mm, frequency $f_s = 20-180$ Hz, and an arbitrary phase shift at release $\phi_s$.

On the vibrating surface, the drop's time to wet $t_w$—the time from initial drop impact to contact initiation—is extended by over 2 orders of magnitude. Compared to a typical time to wet of $t_w \sim 1$ s for a stationary surface [10,12,13], we observe bouncing for several minutes ($t_w \sim 300$ s, movie S1 [33]). In experiments, this prolongation is only limited by the drop's lateral motion, as the drop may encounter occasional defects while roaming the mica surface [30]. Otherwise, our direct numerical simulations (DNS), carried out with the free software Basilisk [48,49] (see Supplemental Material [33]), suggest that a drop retains enough kinetic energy to bounce for an

extended period, and possibly indefinitely, as the energy injected through the air layer balances the viscous dissipation.

**Dynamical regimes**—Typical trajectories depicting the motion of the drop's center-of-mass and the air gap show distinct dynamical states for the same stage amplitude, but with different driving frequencies (movies S2 to S5 [33]). At intermediate frequencies ($f_s = 40$ Hz, $A_s = 0.21$ mm, $\Omega \approx 1.2$, $\Gamma \approx 0.7$), the drop adopts a “bouncing” state with high-amplitude rebounds, leaving a visible air gap during each bounce, as shown in the kymograph in Fig. 1(d), and images in Fig. 1(f), inset. By contrast, at higher frequencies ($f_s = 180$ Hz, $A_s = 0.21$ mm, $\Omega \approx 5.2$, $\Gamma \approx 13.7$), the center-of-mass is hardly displaced, indicating that only a thin air film remains and the drop appears dynamically “bound” to the stage [Figs. 1(e) and 1(g), inset]. In this bound state, the drop's trajectory and deformation [Figs. 1(g) and 1(i)] are locked onto the forcing frequency, with the interaction mediated by a thin air film. In both cases, we observe large drop deformations; in particular,
significant flattening in the bouncing case [Fig. 1(f), inset] and alternating triangular shapes in the bound case [Fig. 1(g), inset] that are accurately captured by the DNS (dash-dotted lines) in Figs. 1(f)–1(i). Figures 1(j) and 1(k) show a normalized energy budget for both cases, with the drop’s total energy given by $T = P + K + S$ ($P$: potential energy, $K$: kinetic energy, and $S$: surface energy). At touchdown, the kinetic energy gained during the initial free fall is transferred into surface energy as the drop spreads, and dissipated through the integrated dissipation energy in the liquid $D$, which must be counterbalanced by the injected stage energy $W$ for a drop to bounce indefinitely (see Supplemental Material [33] for details on energy calculations). We note, in particular, that the dissipation is higher in the bound case than in the bouncing case, suggesting that the drop retains insufficient energy to fully rebound [Figs. 1(j)–1(k) and 2(c)].

**Resonant behavior**—In light of the two very different dynamical states, we explore a wide range of frequencies to better understand the nonlinear response of the system. For a drop radius $r_0 = 0.78 \pm 0.05$ mm with natural frequency $f_d = 35$ Hz and stage amplitude $A_s = 0.21$ mm ($A_s/r_0 \approx 0.3$), we vary the forcing frequency over a range from 20 to 180 Hz corresponding to $0.6 \lesssim \Omega \lesssim 5.2$ and $0.2 \lesssim G \lesssim 13.7$. The drop’s root-mean-square bouncing amplitude $A_d$ (Supplemental Material [33] Sec. S4) is shown as a function of $f_s$ for measurements (black points) and numerical simulations (solid lines) in Fig. 2(a). $A_d$ has a sensitive dependence on the initial impact phase [50] [Fig. 2(b)]; thus, the bouncing amplitude in the DNS is averaged over phase shifts $\phi_s = n\pi/8$ with $n = \{0, 1, \dots, 15\}$. Furthermore, to probe the experimental spread in drop radius (see Fig. S1 [33]), simulation results are shown for three radii $r_0 = \{0.73, 0.78, 0.83\}$ mm (gray, orange, and black lines). The first peak ($f_s \approx f_d$, $\Omega \approx 1$) corresponds to an harmonic response of the system, and, in the second peak ($f_s \approx 2f_d$, $\Omega \approx 2$), the experimental spread for a given $f_s$ comes from the phase shift $\phi_s$ at initial impact [Fig. 2(a)] and lies within one standard deviation of the DNS (orange band). The dependence of $A_d$ on $\phi_s$ persists over longer times (Fig. S3 [33]), suggesting that the drop retains a memory of its initial impact. These two coexisting states around the second peak correspond to sub-harmonic and harmonic responses of the system for high and low bouncing amplitudes, respectively (Figs. S4 and S5 [33]).

At higher frequencies, we observe a sharp bouncing-to-bound transition at a critical frequency $f_c \approx 100$ Hz [Figs. 2(a) and 2(b), orange dash-dotted line for $r_0 = 0.78$ mm], where we have defined the bound state by a low drop amplitude $A_d/r_0 < 0.1$. In addition, the simulations show that $f_c$ depends on drop size, and decreases as $r_0$ increases [Fig. 2(a), dash-dotted lines]. This behavior suggests that the transition is related to the Rayleigh frequency $\omega_d$, as confirmed by the dimensionless

<figure><img src="image_5.png" /><figcaption>FIG. 2. Resonance curve and the bouncing-to-bound transition. (a) Mean drop amplitude $A_d$ for a fixed stage amplitude $A_s = 0.21$ mm. Experimental results with $r_0 = 0.78 \pm 0.05$ mm (black points) compared to phase-averaged DNS results with $r_0 = \{0.73, 0.78, 0.83\}$ mm (gray, orange, and black lines), for a numerical time span of 0.5 s. Proposed bouncing-to-bound cutoff frequencies $f_c \approx \{89, 98, 108\}$ Hz (dash-dotted colored lines) for each radius in the DNS collapse to a single cutoff at $\Omega_c = \sqrt{8}$ in the dimensionless plot (black dash-dotted line in inset). Experimental spread, shown by dashed vertical lines, falls within 1 standard deviation of the numerical average (orange shaded band), and is attributable to variation in $\phi_s$. (b) Influence of $\phi_s$ on $A_d$ in the DNS for the same $A_s$ and varied $f_s$ with $r_0 = 0.78$ mm. (c) Normalized integrated dissipation $D/T_0$ for the same DNS cases as in (a).</figcaption></figure>

plot [Fig. 2(a), inset]. Moreover, the total dissipation sharply increases for forcing frequencies above the critical ones $f_s > f_c$ [Fig. 2(c)], hinting at the excitation of higher order deformation modes in the bound cases.

**Spherical harmonic decomposition**—We decompose the drop’s axisymmetric profile using spherical harmonics, such that

$$r(\theta, t) = r_0 + \sum_{\ell=2}^{\infty} a_{\ell}(t) P_{\ell}(\cos \theta). \quad (1)$$

Here, $\theta$ is the polar angle, $a_{\ell}(t)$ is the time-varying mode amplitude and $P_{\ell}(\cdot)$ the Legendre polynomial of the $\ell$th
<figure><img src="image_6.png" /><figcaption>FIG. 3. Spherical harmonic decomposition of the drop deformation. (a) Illustration of axisymmetric spherical harmonic modes $l = \{2, 3, 4\}$ with amplitudes $a_l/r_0 = \pm 0.3$. (b) Spherical harmonic mode amplitudes of the bouncing case [Fig. 1(d)] with each line color corresponding to one mode illustrated in (a). (c) Detailed view of a single cycle denoted by the shaded region in (b) where the peak negative values of $a_{l=2,3}$ are marked by a point (minimal values over a rebound period). (d), (e) Identical description of the bound case [Fig. 1(e)]. (f) Dimensionless peak amplitudes $-a_{l=2,3}/r_0$ as a function of $\Omega$ with the cutoff frequency $\Omega_c$ as a black dash-dotted line. DNS results are plotted as dash-dotted lines in (b)–(e) and as solid lines in (f) with $r_0 = \{0.73, 0.78, 0.83\}$ mm.</figcaption></figure>

axisymmetric spherical harmonic mode. The sum starts with $l = 2$ because the 0th mode violates mass conservation, while the 1st mode is a translation (not considered in the center-of-mass frame). In practice, the experimental and numerical data are fitted up to $l = 10$, and we find that only the first few modes contribute significantly to the deformation [13] [Fig. 3(a)]. We neglect the air viscosity [51,52], and thus the natural angular frequency $\omega_l$ associated with each mode is

$$ \omega_l = \sqrt{l(l-1)(l+2)}\omega_d. \qquad (2) $$

In the bouncing case [Figs. 3(b) and 3(c)], the drop spreads horizontally (negative $a_2$) before elongating vertically (positive $a_2$); this leads to an upward translation of the center of mass, and to the lift-off of the drop. By contrast, bound drops [Figs. 3(d) and 3(e)] deform into alternating triangular shapes with nearly equal amplitudes for modes

$l = 2$ and $l = 3$, resulting in minimal center-of-mass motion; this prevents escape from the stage. These kinematics reveal that adequate coupling to the $l = 2$ mode is necessary for a drop to completely rebound from the surface; $f_2$ defines a critical transition frequency, corresponding to $\Omega_c = \omega_2/\omega_d = \sqrt{8}$. The predicted cutoff value is in remarkable agreement with our experimental and numerical observations, where the bouncing-to-bound transition frequency collapses to the value $\Omega_c$ for all drop sizes [Fig. 2(a), inset]. Beyond $\Omega_c$, $a_2$ decays whereas $a_3$ increases, corresponding to a limited drop rebound height [Fig. 3(f)].

**Coupled linear spring model**—The pivotal role played by the second spherical harmonic mode motivates the development of a simple coupled linear spring model to predict the drop motion. Whereas most prior models on a liquid bath [20,21,53] neglect the drop deformation, here we couple the dimensionless drop’s center of mass $\tilde{z}(t)$ to its shape, characterized by its dimensionless second deformation mode $\tilde{a}_2(t)$

$$ \left\{ \begin{aligned} \frac{d^2 \tilde{z}}{d\tilde{t}^2} - H(-\tilde{z}^*) [\max\{-\tilde{b}_2 \frac{d\tilde{z}}{d\tilde{t}} - \tilde{z}^*, 0\}] &= -Bo^*, \\ \frac{d^2 \tilde{a}_2}{d\tilde{t}^2} + 2\tilde{b}_2 \frac{d\tilde{a}_2}{d\tilde{t}} + \tilde{\omega}_2^2 \tilde{a}_2 &= H(-\tilde{z}^*) [5(\tilde{\alpha} - \tilde{a}^2/3)], \end{aligned} \right. \qquad (3) $$

with $\tilde{b}_{l=2} = (l-1)(2l+1)\text{Oh} = 5$ Oh the dimensionless dissipation coefficient [51], $\tilde{\omega}_{l=2} = \sqrt{l(l-1)(l+2)} = \sqrt{8}$ the dimensionless natural frequency associated to the second mode, $Bo^* = \text{Bo}[1 + \Gamma \sin(\Omega \tilde{t} + \phi_s)]$ the effective Bond number and $\tilde{\alpha} = d^2\tilde{z}/d\tilde{t}^2 + Bo^*$ the dimensionless instantaneous center-of-mass acceleration. The Heaviside step function $H(-\tilde{z}^*)$ triggers the coupling when $\tilde{z}^* = \tilde{z} - \tilde{a}_2 < 0$, where $\tilde{z}^*$ is the displaced center-of-mass that accounts for the drop deformation (see Appendix for the complete model derivation).

We compare the drop’s trajectory and deformation from the DNS to the model without any fitting parameters. In the model, the dimensionless width $\tilde{w} = 2 - \tilde{a}_2$ and height $\tilde{h} = 2 + 2\tilde{a}_2$ are extracted using Eq. (1). The coupled linear spring can capture the rich drop dynamics in both the bouncing and bound states [Figs. 4(b) and 4(c)]. Over a wide range of driving accelerations and frequencies, the drop amplitude $A_d$ predicted from the model agrees nearly exactly with that of the DNS; furthermore superimposed experimental measurements show excellent agreement with both predictions [Figs. 4(d) and 4(e)]. Bouncing states exist for accelerations as low as $\Gamma = 0.3$, in agreement with similar phase diagrams for drops bouncing on a liquid bath [20]. These results confirm the robustness of the critical frequency $\Omega_c = \omega_2/\omega_d = \sqrt{8}$, above which the drop is bound even for high $\Gamma$ [dash-dotted line in Figs. 4(e) and 4(f)].

**Control and manipulation of drops**—Our results demonstrate a novel hovering state for a liquid drop on a
<figure><img src="image_7.png" /><figcaption>FIG. 4. Predictive model captures bouncing and bound dynamics. (a) Schematic of coupled linear spring model defined in Eq. (3). (b), (c) Comparison of the DNS drop center-of-mass position $z(t)$ (light green line) with the model (dark green line) in a bouncing case with $f_s = 80$ Hz (b) and in a bound case with $f_s = 140$ Hz (c). Insets show the comparison of DNS and model for the drop's width and height variation as a percentage of the original drop diameter. (d) DNS $\Omega-\Gamma$ phase diagram showing the dimensionless drop mean amplitude $A_d/r_0$ with $r_0 = 0.78$ mm, averaged over phase shifts $\phi_s = n\pi/2$ with $n = \{0, 1, 2, 3\}$. Experimental data points are overlaid as circles, with colors indicating measured values of $A_d/r_0$. The cutoff frequency $\Omega_c$ is indicated as a black dash-dotted line. (e) Identical phase diagram computed with the model Eq. (3), showing strong agreement with both numerical simulations and experimental data.</figcaption></figure>

vibrating smooth solid. We identify two distinct regimes depending on $f_s$ and $A_s$—a bouncing and a bound state—with a transition governed by excitation of the drop’s second spherical harmonic mode. From this insight, we develop a coupled linear spring model that predicts drop trajectories without fitting parameters.

Prospective drop manipulation in this no-contact state is facilitated by much more rapid transport than above a liquid layer (Fig. S6 [33]), which has been similarly observed for larger puddles levitated by airflow [54] and Leidenfrost drops above a heated solid [8,9]. Unlike the Leidenfrost effect, the vibration method does not require high temperatures, allowing for handling of sensitive chemistries. As a proof-of-concept, we demonstrate a pneumatic control method [Fig. S7], opening avenues for precision micro-manipulation of liquid in air.

**Acknowledgments**—T. F. acknowledges the Swiss Government Excellence Scholarship ESKAS-Nr: 2023.0357.

J. M. K. conceived the study. L. M. conducted the experiments. T. F. performed the numerical simulations. T. F. and F. G. derived the reduced order model. All authors analyzed the data. L. M. and T. F. drafted the manuscript. All authors contributed to the final version of the manuscript.

*Data availability*—The data that support the findings of this article are openly available [55].

[1] S. Amini, S. Kolle, L. Petrone, O. Ahanotu, S. Sunny, C. N. Sutanto, S. Hoon, L. Cohen, J. C. Weaver, J. Aizenberg, N. Vogel, and A. Miserez, Science **357**, 668 (2017).

[2] S. A. McBride, J. R. Lake, and K. K. Varanasi, J. Chem. Phys. **158**, 134721 (2023).

[3] C. Antonini, M. Innocenti, T. Horn, M. Marengo, and A. Amirfazli, Cold Reg. Sci. Technol. **67**, 58 (2011).

[4] J. Chen, R. Dou, D. Cui, Q. Zhang, Y. Zhang, F. Xu, X. Zhou, J. Wang, Y. Song, and L. Jiang, ACS Appl. Mater. Interfaces **5**, 4026 (2013).

[5] R. Y. Siddiquie, A. Gaddam, A. Agrawal, S. S. Dimov, and S. S. Joshi, Langmuir **36**, 5349 (2020).

[6] W. Barthlott and C. Neinhuis, Planta **202**, 1 (1997).

[7] M. Zhang, S. Feng, L. Wang, and Y. Zheng, Biotribology **5**, 31 (2016).

[8] A.-L. Biance, C. Clanet, and D. Quéré, Phys. Fluids **15**, 1632 (2003).

[9] D. Quéré, Annu. Rev. Fluid Mech. **45**, 197 (2013).

[10] D. Richard and D. Quéré, Europhys. Lett. **50**, 769 (2000).

[11] J. Kolinski, L. Mahadevan, and S. Rubinstein, Europhys. Lett. **108**, 24001 (2014).

[12] J. de Ruiter, R. Lagraauw, D. van den Ende, and F. Mugele, Nat. Phys. **11**, 48 (2015).
[13] J. de Ruiter, R. Lagraauw, F. Mugele, and D. van den Ende, J. Fluid Mech. 776, 531 (2015).

[14] I. Blekhman, L. Blekhman, L. Vaisberg, V. Vasil'kov, and K. Yakimova, Dokl. Phys. 53, 520 (2008).

[15] B. Apffel, F. Novkoski, A. Eddi, and E. Fort, Nature (London) 585, 48 (2020).

[16] B. Apffel, S. Hidalgo-Caballero, A. Eddi, and E. Fort, Proc. Natl. Acad. Sci. U.S.A. 118, e2111214118 (2021).

[17] Y. Couder, E. Fort, C.-H. Gautier, and A. Boudaoud, Phys. Rev. Lett. 94, 177801 (2005).

[18] Y. Couder, S. Protière, E. Fort, and A. Boudaoud, Nature (London) 437, 208 (2005).

[19] S. Protière, A. Boudaoud, and Y. Couder, J. Fluid Mech. 554, 85 (2006).

[20] J. Moláček and J. W. Bush, J. Fluid Mech. 727, 582 (2013).

[21] D. Terwagne, F. Ludewig, N. Vandewalle, and S. Dorbolo, Phys. Fluids 25, 122101 (2013).

[22] Øistein Wind-Willassen, J. Moláček, D. M. Harris, and J. W. M. Bush, Phys. Fluids 25, 082002 (2013).

[23] Y. Couder and E. Fort, Phys. Rev. Lett. 97, 154101 (2006).

[24] A. Eddi, E. Fort, F. Moisy, and Y. Couder, Phys. Rev. Lett. 102, 240401 (2009).

[25] D. M. Harris, J. Moukhtar, E. Fort, Y. Couder, and J. W. M. Bush, Phys. Rev. E 88, 011001(R) (2013).

[26] P. J. Sáenz, T. Cristea-Platon, and J. W. M. Bush, Sci. Adv. 6, eaay9234 (2020).

[27] P. J. Sáenz, G. Pucci, S. E. Turton, A. Goujon, R. R. Rosales, J. Dunkel, and J. W. M. Bush, Nature (London) 596, 58 (2021).

[28] V. Sanjay, P. Chantelot, and D. Lohse, J. Fluid Mech. 958, A26 (2023).

[29] J. M. Kolinski, S. M. Rubinstein, S. Mandre, M. P. Brenner, D. A. Weitz, and L. Mahadevan, Phys. Rev. Lett. 108, 074503 (2012).

[30] R. Kaviani and J. M. Kolinski, Phys. Rev. Fluids 8, 103602 (2023).

[31] J. E. Sprittles, Annu. Rev. Fluid Mech. 56, 91 (2024).

[32] Lord Rayleigh, Proc. R. Soc. 29, 71 (1879).

[33] See Supplemental Material at http://link.aps.org/ supplemental/10.1103/w3qq-cnj3 for details on methods and additional results, which includes Refs. [34–47].

[34] Bluestar Silicones, BLUESIL TM: FLD 47 V 3 bis 47 V 20 (2011), https://my.elkem.com/product2/bluesil-fld-47v20- a2/01tb1000000UD1RQAW.

[35] L. Molefe, G. A. Zampogna, J. M. Kolinski, and F. Gallaire, J. Fluid Mech. 1001, A59 (2024).

[36] A. Daerr and A. Mogne, J. Open Res. Software 4, 3 (2016).

[37] T. J. Senden and W. A. Ducker, Langmuir 8, 733 (1992).

[38] S. Popinet, J. Comput. Phys. 190, 572 (2003).

[39] S. Afkhami and M. Bussmann, Int. J. Numer. Methods Fluids 57, 453 (2008).

[40] S. Afkhami and M. Bussmann, Int. J. Numer. Methods Fluids 61, 827 (2009).

[41] R. Scardovelli and S. Zaleski, Annu. Rev. Fluid Mech. 31, 567 (1999).

[42] S. Popinet and S. Zaleski, Int. J. Numer. Methods Fluids 30, 775 (1999).

[43] G. Tryggvason, R. Scardovelli, and S. Zaleski, Direct Numerical Simulations of Gas-Liquid Multiphase Flows (Cambridge University Press, Cambridge, England, 2011).

[44] J. U. Brackbill, D. B. Kothe, and C. Zemach, J. Comput. Phys. 100, 335 (1992).

[45] J. B. Bell, P. Colella, and H. M. Glaz, J. Comput. Phys. 85, 257 (1989).

[46] B. Zhang, V. Sanjay, S. Shi, Y. Zhao, C. Lv, X.-Q. Feng, and D. Lohse, Phys. Rev. Lett. 129, 104501 (2022).

[47] P. K. Sharma and H. N. Dixit, Phys. Fluids 32, 112107 (2020).

[48] S. Popinet, J. Comput. Phys. 228, 5838 (2009).

[49] S. Popinet, J. Comput. Phys. 302, 336 (2015).

[50] A. Potnis and A. Saha, J. Fluid Mech. 988, A15 (2024).

[51] C. Miller and L. Scriven, J. Fluid Mech. 32, 417 (1968).

[52] A. Prosperetti, J. Fluid Mech. 100, 333 (1980).

[53] K. Okumura, F. Chevy, D. Richard, D. Quéré, and C. Clanet, Europhys. Lett. 62, 237 (2003).

[54] W. Bouwhuis, K. G. Winkels, I. R. Peters, P. Brunet, D. van der Meer, and J. H. Snoeijer, Phys. Rev. E 88, 023017 (2013).

[55] L. Molefe, T. Fullana, F. Gallaire, and J. M. Kolinski, Dataset for “Drops Can Perpetually Bounce over a Vibrating Wettable Solid,” 10.5281/zenodo.16875198 (2025).

[56] J. Moláček and J. W. M. Bush, Phys. Fluids 24, 127103 (2012).

End Matter

Appendix: Model derivation—We now include the complete derivation of the coupled linear spring model. We nondimensionalize the lengths by the drop’s radius $r_0$ and time by the characteristic period $\omega_d^{-1}$ and place ourselves in the frame of reference of the vibrating surface with the stage position fixed at $\tilde{z} = -1$ such that an undeformed drop would be acted upon by the stage when its center-of-mass location is negative ($\tilde{z} < 0$) [left panel of Fig. 4(a)].

Away from the stage, the equation of the center-of-mass position is

$$ \frac{d^2 \tilde{z}}{d \tilde{t}^2} = -Bo^*, \qquad (A1) $$

and the equation of second mode amplitude is

$$ \frac{d^2 \tilde{a}_2}{d \tilde{t}^2} + 2 \tilde{b}_2 \frac{d \tilde{a}_2}{d \tilde{t}} + \tilde{\omega}_2^2 \tilde{a}_2 = 0. \qquad (A2) $$

When the stage acts on the drop, we model the center-of-mass equation as a damped linear spring with a reaction force proportional to the displaced drop center-of-mass
position $(\tilde{z} - \tilde{a}_2)$ [20] and a dissipation term proportional to the drop velocity $d\tilde{z}/d\tilde{t}$, with a dissipation coefficient $\tilde{b}_2$, as $2d\tilde{z}/d\tilde{t} \sim d\tilde{a}_2/d\tilde{t}$ from geometrical considerations. Enforcing that the reaction force is always positive [20], the equation of motion reads

$$ \frac{d^2\tilde{z}}{d\tilde{t}^2} - \max\left\{-\tilde{b}_2 \frac{d\tilde{z}}{d\tilde{t}} - (\tilde{z} - \tilde{a}_2), 0\right\} = -Bo^*. \quad (A3) $$

Considering a quasistatic approximation of drop shapes in contact with a flat solid as derived by Moláček and Bush [56] [see Eqs. (13) and (42)], the second mode amplitude is forced by $f = \tilde{a}(2\ell + 1)/(\ell + 1)P'_2(1 - \tilde{a}/3)$, with $P'_2(\cdot)$ the derivative of the second Legendre polynomial $P_2(x) = \frac{1}{2}(3x^2 - 1)$. Specifying the expression with $\ell = 2$, we obtain that $f = 5(\tilde{a} - \tilde{a}^2/3)$, where

$\tilde{a} = d^2\tilde{z}/d\tilde{t}^2 + Bo^*$ is the dimensionless instantaneous center-of-mass acceleration. The equation of motion for the second mode of deformation when forced by the stage reads

$$ \frac{d^2\tilde{a}_2}{d\tilde{t}^2} + 2\tilde{b}_2 \frac{d\tilde{a}_2}{d\tilde{t}} + \tilde{\omega}_2^2 \tilde{a}_2 = 5(\tilde{a} - \tilde{a}^2/3). \quad (A4) $$

The coupling between the center-of-mass and the second mode of deformation is activated when $\tilde{z}^* = (\tilde{z} - \tilde{a}_2) < 0$. Contrary to the experiments and the DNS, where the drop never contacts the stage, the model allows the deformed drop to intersect the boundary at $\tilde{z} = -1$ [right panel of Fig. 4(a)]. Putting together Eqs. (A1)–(A4), we recover the coupled linear spring model defined in Eq. (3).