# Similarity solutions in elastohydrodynamic bouncing

Vincent Bertin<EMAIL_ADDRESS>Physics of fluids group, University of
Twente, 7500 AE Enschede, The Netherlands

###### Abstract

We investigate theoretically and numerically the impact of an elastic sphere
on a rigid wall in a viscous fluid. Our focus is on the dynamics of the
contact, employing the soft lubrication model in which the sphere is separated
from the wall by a thin liquid film. In the limit of large sphere inertia, the
sphere bounces and the dynamics is close to the Hertz theory. Remarkably, the
film thickness separating the sphere from the wall exhibits non-trivial self-
similar properties that vary during the spreading and retraction phases.
Leveraging these self-similar properties, we establish the energy budget and
predict the coefficient of restitution for the sphere. The general framework
derived here opens many perspectives to study the lubrication film in impact
problems.

††preprint: APS/123-QED

##  I Introduction

The impact or collision of a spherical object on a surface is a problem that
has been of great interest for decades. Typical examples can be found at
various length scales ranging from asteroid impact [1], drop impact [2], ball
sports [3], collision in granular media [4] or suspensions, or more commonly
children playing with rubber bouncing balls. Of particular interest is the
restitution coefficient, defined as the ratio between the bouncing speed
$V_{\infty}$ and the impact speed $V_{0}$. An elastic collision corresponds to
the case of a bouncing speed equal to the impact speed, without any loss of
energy. In real impacts, various effects of the solid tend to decrease the
restitution coefficient at large velocity [5] such as viscoelasticity [6, 7],
plastic deformations [8, 9], or sphere vibrations [10, 11, 12]. On top of
these, viscous dissipation is predominant when the impact occurs in liquid
environments, for instance in particle-laden flows [13], which has a large
variety of industrial and natural applications. Most numerical models of such
flows use a point-particle approximation or immersed boundary and do not
resolve the flow in the lubrication layer upon collision. Instead, empirical
collision laws are implemented when particles start overlapping [14, 15]. In
this article, we aim at describing particle collisions, focusing on the
canonical problem of a soft elastic sphere impacting a rigid surface in a
viscous liquid.

Several experimental studies have been devoted to the individual rebound of a
particle on a surface, either immersed in a liquid [16, 17], or with a thin
oil layer lubricating the surfaces [18, 19, 20]. Head-on and oblique particle-
particle collisions have also been investigated and follow a similar
phenomenology [21, 22]. The restitution coefficient increases with increasing
velocity, as the contact time decreases and the viscous friction dissipates
less energy (see Fig. 1(b)). Remarkably, for a large variety of sphere sizes,
material properties and fluid viscosities, the restitution coefficients
collapse when plotted versus the Stokes number $\mathrm{St}=mV_{0}/(6\pi\eta
R^{2})$, where $m$, $V_{0}$, $\eta$ and $R$ are the sphere mass, impact
velocity, viscosity and sphere radius. The Stokes number needs to exceed a
critical number, which is of order of $10$, for bouncing to occur. Then, the
restitution coefficient increases and reaches asympotically unity when the
Stokes number is $\sim 1000$. An empirical expression
$V_{\infty}/V_{0}=\exp(-35/\mathrm{St})$ has been shown to describe fairly
well the experimental data over the entire range of St [23] (see green dashed
lines in Fig. 1(b)), and is motivated by an analogy with a damped spring-mass
model.

The bouncing dynamics is mainly modeled in two ways. First, reduced linear
damped mass-spring models have been introduced [23]. The non-linearities of
the contact elasticity [24] and gravitational effects [6] have also been taken
into account. The second kind of models described the morphology of the
lubrication layer upon contact. Davis et al. were the first ones to derive an
elastohydrodynamic lubrication model, which describes the intricate coupling
between the thin liquid film and elastic deformations, when particles move
near soft interfaces [25]. The lubrication forces prevent the direct contact
between the sphere and the rigid surface, such that a thin liquid film always
separates the two surfaces. These models allow for a prediction of the
critical bouncing Stokes number, that typically takes the form
$\ln(H/\delta)$, where $H$ is the initial separation distance and $\delta$ is
a cut-off length. The cut-off length is an elastohydrodynamic length for
smooth surfaces, or the typical roughness for rough surfaces [19, 26].
Piezoviscous and compressible effects, which play an important role in
lubricants, have also been discussed [27, 28, 29]. We also point out that the
elastohydodynamic coupling gives rise to a very rich phenomenology for
particle sliding and rotating near soft interfaces [30].

The morphology of the lubrication layer during the contact has rarely been
addressed [31]. The central film thickness has been identified to scale as
$\mathrm{St}^{-1/2}$ in the Ref. [29], suggesting the presence of self-similar
solutions. We present in this article a detailed numerical and asymptotic
analysis of the elastohydrodynamic bouncing. The Stokes number is set in the
range $[10,10^{4}]$ in the numerical simulations of this paper, which
corresponds to the relevant experimental range (see Fig. 1(b)). The structure
of the lubrication layer during the impact dynamics is derived through self-
similar solutions. Furthermore, the energy budget allows to find an asymptotic
expression of the restitution coefficient at large Stokes number.

Figure 1:  (a) Schematic of an elastohydrodynamic bouncing of a soft sphere on
a rigid surface. The underformed sphere is indicated with red dashed lines. A
thin film of thickness $h(r,t)$ prevents direct contact, where $r$ denotes the
axisymmetric radial coordinate. The grey rectangle indicates a zoom in the
contact region. (b) Bouncing velocity as a function of the Stokes number
measured in the experiments of Gondret et al. [17], in the numerical
simulations and predicted by the asymptotic theory of Eq. (27). The prediction
of Legendre et al. [23], using a linear damped mass-spring model is shown in
green dashed lines.

##  II Soft lubrication model

###  II.1 Formulation

We briefly recall here the soft lubrication model, already employed to model
the bouncing of sphere of radius $R$ and mass $m$ on a rigid planar surface
[25, 29, 32]. We assume an impact normal to the surface such that the problem
is axisymmetric, where the axis of symmetry is the sphere axis normal to the
surface, referred later as vertical direction ($z$ in Fig. 1(a)). We introduce
the polar coordinates where the radial position is denoted $r$. The center of
mass vertical position of the sphere, shifted by one radius, is denoted $D(t)$
(see Fig. 1 (a)). We suppose that the sphere is submerged in a viscous fluid
of viscosity $\eta$. Due to lubrication forces, a thin film of liquid always
separates the sphere from the surface and there is no direct contact. Using
Newton’s second law along the vertical direction, the sphere’s dynamical
equation is $m\dot{V}(t)=F(t)$, where $F(t)$ is the the vertical hydrodynamic
viscous force acting on the particle, and $V(t)=\dot{D}(t)$ its vertical
velocity. We focus here on the contact dynamics, such that the sphere is close
enough to the surface to use the lubrication approximation. Therefore, the
flow in the liquid gap is a parabolic Poiseuille flow and the liquid film
thickness $h(r,t)$ follows the thin-film equation

$\frac{\partial h(r,t)}{\partial t}=\frac{1}{12\eta r}\frac{\partial}{\partial
r}\left(rh^{3}(r,t)\frac{\partial p(r,t)}{\partial r}\right),$ (1)

where $p(r,t)$ is the pressure field. Within the lubrication approximation,
the contribution of the outer region of the contact to the hydrodynamic
viscous force is negligible, and the hydrodynamic viscous force mainly comes
from the contact region and follows the lubrication force as
$F(t)=\int_{0}^{\infty}p(r,t)\,2\pi r\mathrm{d}r$. Additionally, the sphere
being very close to the surface, the underformed spherical particle can be
modelled with a parabolic approximation such that the liquid gap can be
written as [33]

$h(r,t)=D(t)+\frac{r^{2}}{2R}-w(r,t),$ (2)

where $w(r,t)$ represents the elastic deformations (see Fig. 1 (a)). These
deformations are modeled by using the linear elasticity theory and are related
to the pressure field with a convolution integral involving the elastic
Green’s function. Integrating the convolution integral in the azimuthal
direction leads to the integral relation [25]

$w(r,t)=-\frac{4}{\pi
E^{*}}\int_{0}^{\infty}\mathcal{M}(r,x)p(x,t)\,\mathrm{d}x,\quad\quad\mathcal{M}(r,x)=\frac{x}{r+x}K\left(\frac{4rx}{(r+x)^{2}}\right),$
(3)

where $E^{*}=E/(1-\nu^{2})$ is the plane strain elastic modulus with $E$ and
$\nu$ the Young’s modulus and Poisson ratio respectively. The function $K$ is
the complete elliptic integral of the first kind.

We introduce the typical sphere elastic deformation length $D_{0}$ during the
contact,

$D_{0}=\left(\frac{mV_{0}^{2}}{E^{*}\sqrt{R}}\right)^{2/5},$ (4)

obtained by balancing the kinetic energy $mV_{0}^{2}$ with the Hertz elastic
energy $E^{*}\sqrt{R}D_{0}^{5/2}$ (see below). Equations (1)-(3) are made
dimensionless, by using the scales of the Hertz theory, respectively $D_{0}$,
$\sqrt{RD_{0}}$, $D_{0}/V_{0}$, $p_{0}=2E^{*}\sqrt{D_{0}/R}/\pi$ and
$F_{0}=p_{0}RD_{0}$ as vertical/radial length, time, pressure and force
scales, as detailed in Appendix A. We then solved numerically the resulting
equations using the finite-difference scheme proposed in Ref. [34], as
discussed in Appendix B. As an initial condition, we suppose that the sphere’s
position is at a height $D_{0}$ with a downward velocity $V_{0}$. A single
dimensionless number characterizes the bouncing dynamics, which is the Stokes
number

$\mathrm{St}=\frac{mV_{0}}{6\pi\eta R^{2}}.$ (5)

The Stokes number can be interpreted as the ratio between sphere inertia
($mV_{0}^{2}$) to viscous dissipation ($6\pi\eta R^{2}V_{0}$), or
alternatively the viscous dissipation time scale $mD_{0}/(6\pi\eta R^{2})$ to
the bouncing time scale $D_{0}/V_{0}$. For small Stokes number (large
viscosity), or dissipation time smaller than the bouncing time, the sphere
does not bounce as the entire initial kinetic energy is dissipated before
contact. Therefore, in the context of bouncing, the Stokes number offers a
quantification of the viscous dissipation during bouncing. The critical Stokes
number for bouncing to occur depends slightly on the initial altitude and is
typically of the order of $10$ [25] (see Fig. 1 (b)). In this article, we
focus on the intermediate to large Stokes number regime (low viscosity fluid),
when the sphere bounces with a non-zero speed. We neglect the buoyancy forces
here that do not influence on the contact dynamics as long as the work done by
buoyancy during contact $\Delta\rho R^{3}gD_{0}$ is small compare to the
kinetic energy $mV_{0}^{2}$ [6], where $\Delta\rho$ is the density difference
between the sphere and the fluid. Lastly, we suppose that the impact speed is
not too large, such that the typical sphere deformation $D_{0}$ is small
compared to the sphere radius, corresponding to the condition
$(mV_{0}^{2}/(E^{*}R^{3}))^{2/5}\ll 1$.

In the following sections, we will find various asymptotic similarity
solutions for the film thickness, each with a specific scale different from
the typical deformation $D_{0}$ (see Eq. (4)). We shall notice that in the
infinite Stokes limit, the elastohydrodynamic model converges toward the
elastic collision, where the bouncing velocity is equal to the impact
velocity, as discussed in the next subsection.

##  III Impact dynamics

###  III.1 Dry limit, $\mathrm{St}\to\infty$: Hertz theory

In the absence of any surrounding fluid, or for $\mathrm{St}\to\infty$, an
exact solution of the elastic bouncing dynamics has been introduced previously
[5, 35]. The vertical force acting on the sphere is zero as long as the sphere
is out of contact, i.e. $D(t)>0$. During the contact phase, the force follows
the Hertz theory $4E^{*}R^{1/2}\delta^{3/2}(t)/3$, where we introduce
$\delta(t)=-D(t)$ as the (positive) indentation. Substituting this expression
in the Newton’s second law for the sphere, the indentation follows the
ordinary differential equation

$m\ddot{\delta}(t)=-\frac{4E^{*}}{3}R^{1/2}\delta^{3/2}(t),\quad\text{for
}\delta(t)>0,$ (6)

which corresponds to a nonlinear mass-spring equation, where the spring
stiffness increases with the indentation as $\sqrt{\delta}$. An exact solution
of Eq. (6) can be expressed in the form of an implicit equation

$\int_{0}^{\delta/D_{0}}\frac{\mathrm{d}x}{\sqrt{1-\frac{16}{15}x^{5/2}}}=\,_{2}F_{1}\left(\frac{2}{5},\frac{1}{2};\frac{7}{5};\frac{16(\delta(t)/D_{0})^{5/2}}{15}\right)\frac{\delta(t)}{D_{0}}=\frac{V_{0}t}{D_{0}},$
(7)

where ${}_{2}F_{1}$ is the Hypergeometric function. We notice that the
equivalent dynamical law for a linear spring model is
$\arcsin(\delta(t)/D_{0})=V_{0}t/D_{0}$, solution of Eq. (7) if one replaces
$16/15x^{5/2}$ by $x^{2}$ in the integral. The shape of Eq. (7) is fairly
similar to the usual sinusoidal law for harmonic oscillators. Furthermore, the
Hertz theory predicts that the pressure and elastic deformations are
piecewise-defined functions that read

$\displaystyle p_{\mathrm{Hz}}(r,t)$
$\displaystyle=\left\\{\begin{array}[]{l}\frac{2E^{*}}{\pi
R}\sqrt{a^{2}(t)-r^{2}},\quad r<a(t),\\\ 0,\quad r>a(t),\\\
\end{array}\right.$ (8c) $\displaystyle w_{\mathrm{Hz}}(r,t)$
$\displaystyle=\left\\{\begin{array}[]{l}-\delta(t)+\frac{r^{2}}{2R},\quad
r<a(t),\\\ -\frac{1}{\pi
R}\left[(2a^{2}(t)-r^{2})\arcsin\left(\frac{a(t)}{r}\right)+a(t)\sqrt{r^{2}-a^{2}(t)}\right],\quad
r>a(t),\\\ \end{array}\right.$ (8f)

where $a(t)=\sqrt{\delta(t)R}$ is the contact radius. The equations (8) are
solutions of the elasticity equation (Eq. (3)), satisfying the boundary
conditions of the Hertz contact problem: (i) stress free condition outside the
contact area $r>a$, and (ii) imposed elastic deformations in the contact zone
$r<a$ with a parabolic shape. Injecting the Hertz elastic deformation in the
film thickness definition Eq. (2), we can express the film thickness predicted
the Hertz theory

$h_{\mathrm{Hz}}(r,t)=\left\\{\begin{array}[]{l}0,\quad r<a(t),\\\
-\delta(t)+\frac{r^{2}}{2R}+\frac{1}{\pi
R}\left[(2a^{2}(t)-r^{2})\arcsin\left(\frac{a(t)}{r}\right)+a(t)\sqrt{r^{2}-a^{2}(t)}\right],\quad
r>a(t),\end{array}\right.$ (9)

that is indeed zero in within the contact zone $r<a(t)$. The limiting behavior
of the Hertz pressure and film thickness in the vicinity of the contact radius
is given by

$h_{\mathrm{Hz}}(r\to a^{+})=\frac{8\sqrt{2}\sqrt{a(t)}}{3\pi
R}\left(r-a(t)\right)^{3/2},\quad\quad p_{\mathrm{Hz}}(r\to
a^{-})=\frac{2E^{*}}{\pi R}\sqrt{2a(t)}\left(a(t)-r\right)^{1/2}.$ (10)

All these results will turn out relevant as “outer solutions” for the impact
dynamics at finite Stokes number.

###  III.2 Finite Stokes number

Figure 2:  Bouncing dynamics: (a) Snapshots of the soft sphere interface
during bouncing, which illustrates the different phases of the dynamics. The
dimensionless times from left to right are
$V_{0}t/D_{0}=0.8,1.2,2.0,2.8,3.5,4.1,4.3$ and the Stokes number is $100$.
Variation of the dimensionless sphere center of mass position (b), velocity
(c) and lubrication force (d), versus the dimensionless time for two Stokes
numbers $\mathrm{St}=50$ and $200$. The black dashed lines indicate the Hertz
theory (7), corresponding to the absence of viscous dissipation
($\mathrm{St}=\infty$). The points illustrate the five characteristic times
separating the different phases for the case $\mathrm{St}=50$ (see definition
in the text). The inset in the panel (d) displays a zoom on the viscous
adhesion phase, where both the horizontal and the vertical axis are rescaled
with $\mathrm{St}^{-2/5}$ and $\mathrm{St}^{-3/5}$ respectively. We also point
out that the origin of time has been shifted in the inset by $t_{3}$ where the
force vanishes.

The bouncing dynamics at finite Stokes number is illustrated in Fig. 2 (a)
(see also Supp. Video 1 [36]) for $\mathrm{St}=100$, corresponding to the
numerical solution of Eqs. (1)-(3). We identify five different stages during
the bouncing, namely the approach, spreading, retraction, adhesion and
bouncing phases. These are illustrated in Figure 2(a), where we also introduce
times $t_{i}$, for $i\in[1,5]$, that separate the phases. Initially, the
sphere approaches the surface and the elastic deformation remains small, i.e.,
of order $\mathrm{St}^{-1}$. After the time of contact $t_{1}=D_{0}/V_{0}$,
where the sphere would have touched the surface in the absence of surrounding
fluid, the elastic deformation gradually increases as if the sphere spreads on
the surface. The sphere reaches its maximal deformation at time $t=t_{2}$ and
its vertical velocity then changes sign, such that the contact radius now
retracts. The time $t=t_{3}$ corresponds to the instant where the lubrication
force vanishes, and is followed by a phase where a negative force is exerted
on the sphere such that the sphere effectively adheres to the surface.
Finally, at $t_{4}$, the sphere takes off and enters the bouncing phase, where
the elastic deformations are small. We define empirically $t_{4}$ as the local
inflexion point of the total energy of the system versus time curve, that
turns out to be a good estimation of the take-off time. The fifth time $t_{5}$
is defined as the time at which the sphere reaches back its original altitude
and bouncing velocity $V_{\infty}=V(t_{5})$.

The Figures 2(b)-(d) quantify the global bouncing dynamics via the sphere
vertical position, velocity and lubrication force for the cases
$\mathrm{St}=50$ and $200$. The specific times $t_{i}$, for $i\in[1,5]$, are
indicated with circles only for the case $\mathrm{St}=50$. The sphere dynamics
is fairly close to the Hertz theory of Eq. (7) (see black dashed lines in Fig.
2 (b)-(d)). The larger the Stokes number, the closer is the dynamics to Eq.
(7), as the viscous effects are less pronounced. Nevertheless, an important
difference is that the velocity after the bounce is smaller than the impact
velocity, i.e. $V_{\infty}/V_{0}<1$, due to the presence of dissipation.
Figure 1 (b) displays the coefficient of restitution $V_{\infty}/V_{0}$ versus
the Stokes number. The numerical results of coefficient of restitution
compares fairly well with the experimental results of Gondret et al. [17]
obtained for a variety of spheres of different materials/sizes and in various
liquids.

A remarkable feature occurs at the end of the contact, where the lubrication
force becomes negative, indicating some effective adhesion. Indeed, the sphere
appears to briefly “stick” to the surface (see adhesion phase in Fig. 2 (a)).
We stress that there is no surface adhesion here, and the negative force
arises from the viscous liquid, specifically due to the viscous resistance
during the filling of the liquid gap. Hence, the adhesion force here cannot be
described by JKR-like theory, but rather by the viscous adhesion that is also
called Stefan adhesion [37]. The adhesion becomes less pronounced at large
Stokes number, highlighting the viscous nature of the force. Interestingly, an
empirical collapse of the adhesion force is obtained when rescaling the force
by $\mathrm{St}^{-3/5}$ and time by $\mathrm{St}^{-2/5}$ (see inset in Fig. 2
(d)), suggesting some universality of the viscous adhesion [37]. The re-
scaling leads to dimensional force and time scales respectively given by $\eta
RV_{0}\left[\eta V_{0}/(RE^{*})\right]^{-2/5}$ and $R/V_{0}\left[\eta
V_{0}/(RE^{*})\right]^{2/5}$, which are mass free and correspond to the usual
elastohydrodynamic lubrication scales (see discussion in Appendix C). A
precise description of the flow and the elastohydrodynamic coupling in the
adhesion phase is left for future work.

The Hertz model is dissipation free, and thus cannot predict the coefficient
of restitution. To understand the viscous dissipation during bouncing, we
extend the Hertz theory and analyse the lubrication film thickness during the
contact. We first focus on the spreading phase of the dynamics in section IV,
and the retraction phase will be discussed in Section V. The main assumption
is that the Stokes number is large, so that the dynamics is close to Hertz
theory.

##  IV Spreading phase

The film thickness and pressure profiles are shown in Fig. 3(a) and (c)
respectively for three different times during the spreading phase for the case
of St=1000. The numerical pressure profiles agree very well with the solution
of the Hertz theory Eq. (8) within the contact region $r<a$. Hertz theory
neglects the lubrication pressure. However, the pressure gradient
$\mathrm{d}p/\mathrm{d}r$ is singular at the edge of the contact radius $r=a$
(see Eq. (10)), so that at finite velocity, the film thickness must be small
in order to balance Hertz and lubrication pressures. This gives rise to small
zone near the contact radius, smoothing out the pressure gradient, as shown in
Fig. 3(c). Moreover, the liquid film profile presents a dimple shape, similar
to drop impact [38], or the profiles obtained during the thin-film drainage
between droplets or bubbles [39, 40]. We separate the theoretical analysis of
the film into two zones: namely the edge of the contact called “neck region”
and the central region called “dimple region” (see schematics in Fig. 3).

###  IV.1 Neck region: soft lubrication inlet analogy

Figure 3:  Spreading phase: neck solution at different times. Typical
dimensionless film thickness (a) and pressure (c) as a function of the
dimensionless radius for three different times (resp. $t=1.3,\,1.5$ and
$1.8\,D_{0}/V_{0}$) during the spreading phase. The Stokes number is set to
$1000$. The black dashed lines show the Hertz theory. In (b) and (d), the
profiles are rescaled by the typical length and pressure scales in the neck
region, corresponding to Eq. (12). The different lateral scales of the problem
are shown in the schematic on top. The soft slider solution of Snoeijer et al.
[41] is shown in pink dashed lines.

The impact dynamics of a soft sphere exhibit resemblances with the problem of
an elastic sphere sliding near a wall under an applied load [29]. Indeed,
using as a reference frame $r\rightarrow r-a(t)$, where $a(t)$ is the contact
radius of the Hertz theory, we can map the impact dynamics in the neck to the
sliding of a soft sphere near a rigid wall, with a time-dependent velocity
$u=\dot{a}$ and under a load $4E^{*}a^{3}/(3R)$. In Refs. [42, 43, 41], it has
been shown that the steady soft sliding has self similar solutions near the
edge of the contact radius, that can be demonstrated rigorously using
asymptotic matching methods. Here, we first recall the typical scales in the
region using the scaling arguments [44]. To do so, we introduce
$h_{n},\ell_{n}$ and $p_{n}$ as the typical film thickness, lateral length and
pressure scales respectively in the neck region (see schematic in Fig. 3).
First, the fluid momentum balance, i.e. Stokes equation
$\eta\nabla^{2}\mathbf{u}=\nabla p$, in the horizontal direction yields to
$\eta u/h^{2}_{n}=p_{n}/\ell_{n}$ in the lubrication limit. Then, following
the linear elasticity theory, the pressure is proportional to the typical
strain which reads $p_{n}=E^{*}h_{n}/\ell_{n}$. Finally, the similarity
solution has to match the singular behavior of the Hertz theory near the
contact radius (see Eq. (10)). This imposes the geometric condition
$h_{n}=a^{1/2}\ell_{n}^{3/2}/R$. Combining the aforementioned expressions, one
finds

$h_{n}=\frac{a^{2}}{R}\lambda^{3/5},\quad\quad\ell_{n}=a\lambda^{2/5},\quad\quad
p_{n}=E^{*}\frac{a}{R}\lambda^{1/5},\quad\quad\lambda=\frac{\eta
uR^{3}}{E^{*}a^{4}}.$ (11)

The dimensionless group $\lambda$ is the relevant elastohydrodynamic
dimensionless number, which needs to be small to enforce the hierarchy of
scales $h_{n}\ll\ell_{n}\ll a\ll R$ of the asymptotic expansion. Following
this approach, we introduce a self-similarity ansatz of the form

$h(r,t)=h_{n}\mathcal{H}(\xi),\quad\quad
p(r,t)=p_{n}\mathcal{P}(\xi),\quad\quad\xi=\frac{r-a}{\ell_{n}}.$ (12)

Substituting Eq. (12) into Eq. (1), and using $\lambda\ll 1$, one gets

$\left(\frac{\dot{h}_{n}\ell_{n}}{h_{n}\dot{a}}\mathcal{H}-\frac{\dot{\ell}_{n}}{\dot{a}}\xi\mathcal{H}^{\prime}\right)-\mathcal{H}^{\prime}=\frac{1}{12}\left(\mathcal{H}^{3}\mathcal{P}^{\prime}\right)^{\prime},$
(13)

where the term in bracket in the left hand side of Eq. (13) reflects the
unsteadyness of the problem, given that both the velocity and load are time
dependent. In the early times after contact, i.e. when $t-t_{1}$ is small, the
vertical velocity is roughly constant $V\approx-V_{0}$, such that the contact
radius dynamics is governed by $a(t)\approx\sqrt{RV_{0}(t-t_{1})}$ and
$u(t)\approx\sqrt{RV_{0}/[4(t-t_{1})]}$. Hence, the elastohydrodynamic
parameter, which scales as $\lambda\propto(u/a^{4})\propto(t-t_{1})^{-5/2}$,
diverges as $t\to t_{1}$. The steady similarity solution develops after a
transient regime, once $\lambda\ll 1$. The typical scale of the transient time
can be estimated as the one when $\lambda=1$, leading to
$(t-t_{1})\sim(D_{0}/V_{0})\mathrm{St}^{-2/5}$, which is much smaller than the
bouncing time scale $D_{0}/V_{0}$ in the large $\mathrm{St}$ limit. We recover
the elastohydrodynamic lubrication time scales, also governing the adhesion
phase (see Appendix C). Hence, after this very brief transient state, the
unsteady terms of Eq. (13) are negligible and we can integrate to obtain

$\mathcal{P}^{\prime}=12\frac{\mathcal{H}_{0}-\mathcal{H}}{\mathcal{H}^{3}},$
(14)

where $\mathcal{H}_{0}$ is an integration constant. In the limit where the
length of the neck is smaller than the contact radius, i.e. $\ell_{n}\ll a$,
the elastic kernel $\mathcal{M}$ of Eq. (3) reduces to the dimensionless line
force Green’s function
$\mathcal{M}(r,r^{\prime})\propto-\frac{1}{2}\ln\mid\xi-\xi^{\prime}\mid$ [5].
Taking (twice) the derivative of Eq. (2), and combining with the line force
elastic kernel, we find

$\mathcal{H}^{\prime\prime}(\xi)=-\frac{2}{\pi}\int_{\mathbb{R}}\frac{\mathcal{P}^{\prime}(\xi^{\prime})}{\xi-\xi^{\prime}}\,\mathrm{d}\xi^{\prime},$
(15)

where we assume $\ell_{n}^{2}/(h_{n}R)\ll 1$, which follows from Eq. (11). The
similarity solution has to match the asymptotic behavior of the Hertz theory
near the contact radius (see Eq. (10)). Writing the matching condition, we
find

$\mathcal{H}(\xi\to\infty)=\frac{8\sqrt{2}}{2\pi}\xi^{3/2},\quad\quad\mathcal{P}(\xi\to-\infty)=\frac{2\sqrt{2}}{\pi}(-\xi)^{1/2}.$
(16)

The equations (14)-(16) are equivalent to the ones of the steady EHD sliding
in the inlet zone of Snoeijer et al. [41], up to a trivial prefactor
rescaling. As shown in Fig. 3 (b) and (d), the rescaled profiles with the
similarity variables of Eq. (12) indeed collapse for different times.
Furthermore, the similarity solution derived in [41] is in perfect agreement
with the profile, demonstrating that the structure of the advancing neck in
the bouncing problem maps perfectly to the inlet of the soft slider, in a
quasi-steady fashion. It is not trivial that the advancing neck region in the
bouncing problem is equivalent to the inlet of the steady soft sliding for two
reasons. First, the bouncing problem is unsteady, and second it implies normal
motion of the sphere versus lateral motion for the steady soft sliding.

Figure 4:  Spreading phase: Stokes number scaling. (a) Typical dimensionless
film thickness as a function of the dimensionless radius at $t=1.5D_{0}/V_{0}$
during the spreading phase. The three colors indicate different Stokes number,
respectively $50,\,200,\,\text{and}\,1000$ and the black dashed lines the
Hertz theory. In (b) (resp. (c)), the thickness profiles are rescaled by the
typical length scales in the neck (resp. dimple) region. The inset shows the
thickness at the Hertz contact radius (b) and the central film thickness (c)
versus the Stokes number in log-log, highlighting the Stokes number scaling
with a fitted line.

The analysis above shows how at a fixed Stokes number, the spreading phase can
be understood in time, by the analogy with a sliding contact problem. It is
also of interest to investigate the scaling of the neck thickness upon varying
the Stokes number – in particular since we expect the neck to vanish in the
limit $\mathrm{St}\to\infty$. In figure 4(a), we thus plot profiles at
different $\mathrm{St}$, at a fixed time ($t=1.5D_{0}/V_{0}$) during the
spreading phase. As expected, the larger the Stokes number, the thinner is the
lubrication film and the closer it gets to the Hertz theory. The neck film
thickness and lateral scales (see Eq. (11)) can be rewritten to make explicit
the Stokes number, as

$\frac{h_{n}}{D_{0}}=(6\pi\mathrm{St})^{-3/5}\left(\frac{uD_{0}}{V_{0}\sqrt{RD_{0}}}\right)^{3/5}\left(\frac{a}{\sqrt{RD_{0}}}\right)^{-2/5},\quad\frac{\ell_{n}}{\sqrt{RD_{0}}}=(6\pi\mathrm{St})^{-2/5}\left(\frac{uD_{0}}{V_{0}\sqrt{RD_{0}}}\right)^{2/5}\left(\frac{a}{\sqrt{RD_{0}}}\right)^{-3/5},$
(17)

Hence, we find that the Stokes number scaling of the film thickness and
lateral scale of the neck during the spreading phase are $\mathrm{St}^{-3/5}$
and $\mathrm{St}^{-2/5}$ respectively. Rescaling the profiles of Fig. 4(a)
with the corresponding Stokes number scaling, a collapse of film thickness
profiles is indeed observed in Fig. 4(b) close to the neck region.
Furthermore, the inset shows the film thickness at the radial position of the
Hertz contact radius which is in perfect agreement with $\mathrm{St}^{-3/5}$
over the full $\mathrm{St}$ number range explored numerically.

###  IV.2 Central region: dimple height model

Figure 5:  Dimple height. (a) Rescaled central film thickness height
$h_{0}(t)=h(r=0,t)$ as a function of time in both the spreading and retraction
phases. The colors indicate different Stokes numbers (same as in Fig. 4),
respectively $50$, $200$ and $1000$ from light to dark blue. As in Fig. 2, the
dots indicates the times separating the different phases of the bouncing
dynamics. The panel (b) shows the same data, where the vertical axis is
rescaled by $\mathrm{St}^{-1/2}$. The prediction of the dimple height Eq. (19)
is shown in red dashed lines, and the small time-to-contact asymptotic Eq.
(20) is displayed in yellow green line.

Toward the dimple region, the neck similarity solutions deviates from the
steady EHD inlet (see Fig. 3 (b)), that has a uniform thickness. Instead, in
the bouncing problem, the film thickness has some spatial variations in the
central region, taking the form of a dimple. The central region is not
described by the same scaling laws as the neck region, and a different
analysis must be adopted. We focus on the central height of the film
$h_{0}(t)=h(r\to 0,t)$, that is referred as dimple height in what follows. The
Figure 5(a) reports the dimple height as a function of time for various Stokes
numbers. As discussed above, the larger the Stokes number and the smaller the
film thickness, as the lubrication pressure gets smaller. Interestingly, the
dimple heights for different Stokes numbers collapse when rescaled by
$\mathrm{St}^{-1/2}$ (see Fig. 5(b)). In Fig. 4(c), we show the radial
dependence of the film thickness profiles rescaled by $\mathrm{St}^{-1/2}$,
where the profiles indeed collapse near the center. Nevertheless, the profiles
deviate away from the symmetry axis. Therefore, the $D_{0}\mathrm{St}^{-1/2}$
thickness scale is only valid in the vicinity of the symmetry axis and does
not collapse the film thickness in the entire dimple.

We turn to an analysis of the dimple height. The pressure fields is well
described by the Hertz profile Eq. (8c) (see Fig. 3(c)). Injecting the
Hertzian pressure field in the thin film equation Eq. (1) and investigating
the limit $r\to 0$, we obtain an ordinary differential equation for $h_{0}(t)$
as

$\left.\frac{\partial h(r,t)}{\partial t}\right|_{r\to
0}=\frac{\mathrm{d}h_{0}(t)}{\mathrm{d}t}=\left.\frac{1}{12\eta
r}\frac{\partial}{\partial r}\left(rh^{3}(r,t)\frac{\partial
p_{\mathrm{Hz}}(r,t)}{\partial r}\right)\right|_{r\to
0}=-\frac{E^{*}}{3\pi\eta R}\,\frac{h_{0}^{3}(t)}{a(t)}.$ (18)

Solving this equation with the initial condition $h_{0}(t\to
t_{1}^{+})=\infty$, we get

$h_{0}(t)=\left(\frac{2E^{*}}{3\pi\eta
R}\int_{t_{1}}^{t}\frac{\mathrm{d}\hat{t}}{a(\hat{t})}\right)^{-1/2}=D_{0}\mathrm{St}^{-1/2}\left(4\int_{Vt_{1}/D_{0}}^{Vt/D_{0}}\frac{\mathrm{d}\hat{t}}{a(\hat{t})/\sqrt{RD_{0}}}\right)^{-1/2}.$
(19)

Not only we recover the $\mathrm{St}^{-1/2}$ scaling, but the dynamics of the
dimple height is quantitatively described by Eq. (19) as shown in Fig. 5(b).
At small time-to-contact, i.e. $t-t_{1}\ll D_{0}/V_{0}$, the Hertz contact
radius follows $a(t)=\sqrt{RV_{0}(t-t_{1})}$. In this limit, the dimple height
expression simplifies to

$h_{0}(t)\simeq\left(\frac{4E^{*}}{3\pi\eta
R^{3/2}V_{0}^{1/2}}(t-t_{1})^{1/2}\right)^{-1/2}\propto(t-t_{1})^{-1/4}.$ (20)

Interestingly, the small time-to-contact asymptotic expression provides an
excellent description of the full expression, as shown with the yellow dotted
lines in Fig. 5(b). We point out that the $-1/2$ scaling in Stokes number of
the dimple height has already been discussed by Venner et al. [29].

##  V Retraction phase

Figure 6:  Retraction phase. Typical dimensionless film thickness (a) and
pressure (c) as a function of the dimensionless radius for three different
times (resp. $t=2.9,\,3.2$ and $3.5\,D_{0}/V_{0}$) during the retraction
phase, using the same notations as in Fig. 3. In (b)-(d), the profiles are
rescaled by the typical length and pressure scales in the neck during the
retraction phase, corresponding to Eq. (21). In contrast to the spreading
phase, there is no universal behavior, although the collapse is fairly good.

We now turn to the retraction phase of the bouncing dynamics, when the sphere
vertical velocity is positive. Figure 6(a) displays the film thickness and
pressure profiles at various times during the retraction phase. The neck
appears to be translated with minor changes of its vertical or lateral scales,
which differs qualitatively from the spreading phase (see Fig. 3(a)). The
central region is fairly flat, as if the dimple has disappeared and the
central thickness decreases very slowly in time, following Eq. (19).

Given the striking agreement between the neck region in the spreading phase
with the soft slider inlet, we use the same approach to the retraction phase.
The contact radius is now receding, such that the analogy must involve the
outlet region of the soft slider. We stress that the central region of the
soft slider has a uniform film thickness $\mathcal{H}_{0}h_{n}$ that is
universal and selected by the inlet profile, where $\mathcal{H}_{0}=2.478...$
is a numerical prefactor [41]. For the soft slider, the very same scaling
arguments of section IV.1 still apply in the outlet and a similarity solution
can be found with the same scales as Eq. (11). For the impact, however, the
situation is different: rescaling the numerical profiles in the neck with the
scales Eq. (11) does not collapse the data during the retraction phase. For
instance, Fig. 7(a) shows film thickness profiles at a fixed time
$t=3.2D_{0}/V_{0}$ for various Stokes numbers. The data are rescaled in Fig.
7(b) in the same manner as in the spreading phase (see Fig. 4(b)), where we no
longer observe a perfect collapse. The difference is further illustrated in
Fig. 7(c) where the value of the film thickness at the Hertz contact radius is
displayed versus the Stokes number in log-log scale. In the range explored
numerically, the film thickness does not follow the $\mathrm{St}^{-3/5}$
scaling as in the spreading phase, but is closer to the $\mathrm{St}^{-1/2}$
scaling of the dimple.

The fundamental difference between the inlet and outlet solution of the soft
slider is that the inlet has a unique solution with a universal dimensionless
flux $\mathcal{H}_{0}$, while the outlet solution has a family of solutions
with varying fluxes, where the adopted self-similar shape is found by matching
to the central region thickness [41]. We note that such a qualitative feature
is also found in the Bretherton bubbles in a tube [45], where the film
thickness is selected from the universal solution in the front dynamical
meniscus, while the rear adopts the same self-similarity but where the
solution is non-unique and selected by matching the film thickness to the
central solution. Hence, we conjecture that the typical thickness scale of the
neck region in the retraction phase is selected by the matching to the central
film thickness $h_{0}$, and not by the elastohydrodynamic film thickness
$h_{n}$. As a result, the similarity solution is no longer universal and
depends on a dimensionless number which is the instantaneous ratio between
these two thicknesses $h_{0}/h_{n}$, that is actually large given the
$\mathrm{St}$ scaling. For $\mathrm{St}=10^{3}$, we observed numerically that
the central film thickness in the retraction phase is roughly constant in time
and scales as $h_{0}\sim D_{0}\mathrm{St}^{-1/2}$ (see Fig. 6(a)). Hence, we
suppose that the film thickness scale in the neck region during the retraction
is $h_{r}=D_{0}\mathrm{St}^{-1/2}$, where the subscript $r$ stands for
retraction. As in section IV.1, we invoke the matching condition of the
retracting neck solution to the singular behavior of Hertz theory (see Eq.
(10)) to determine the lateral length and pressure scales in the neck as
$h_{r}=\ell^{3/2}_{r}a^{1/2}/R$ and $p_{r}=E^{*}a^{1/2}\ell_{r}^{1/2}/R$.
Combining the latter expressions, we find the scales

$h_{r}=D_{0}\mathrm{St}^{-1/2},\quad\quad\ell_{r}=\sqrt{RD_{0}}\mathrm{St}^{-1/3}\left(\frac{a}{\sqrt{RD_{0}}}\right)^{-1/3},\quad\quad
p_{r}=p_{0}\mathrm{St}^{-1/6}\left(\frac{a}{\sqrt{RD_{0}}}\right)^{1/3}.$ (21)

Once rescaled with Eq. (21), the thickness and pressure fields in the neck
collapse fairly well, as shown in Fig. 6(b) and (d). Nevertheless, as
expected, there is no clear universal self-similar solution, and some details
of the profiles are not described by the rescaling. Further investigations
would be necessary to determine the exact asymptotic structure of the solution
at large Stokes number. Such an analysis would required Stokes number much
larger than $10^{5}$ in order to distinguish between the two exponents.

Figure 7:  Retraction phase: Stokes number scaling. (a) Typical dimensionless
film thickness as a function of the dimensionless radius at $t=3.2D_{0}/V_{0}$
during the spreading phase. The four colors indicate different Stokes number,
respectively $50,\,200,\,1000\,\text{and}\,10000$ and the black dashed lines
the Hertz theory. In (b), the thickness profiles are rescaled by the typical
length scales in the neck region during the spreading phase (see Fig. 4). The
panel (c) shows the thickness at the Hertz contact radius versus the Stokes
number in log-log. The two lines indicates power laws with exponents $-3/5$
and $-1/2$ respectively.

##  VI Energy dissipation and coefficient of restitution

Figure 8:  Energy budget. (a) Decay of the rescaled total energy of the sphere
as a function of the dimensionless time. The colors indicate different Stokes
numbers (same as in Fig. 4), respectively $50$, $200$ and $1000$ from light to
dark blue. The energy dissipated in each phase is denoted $\Delta E_{i}$, for
$i\in[1,5]$ and shown in (b)-(f) as a function of the Stokes number. The black
dashed lines correspond to the fit of the numerical data with the asymptotic
predictions in the large-Stokes limit.

Having discussed in details the evolution of the lubrication film during the
contact dynamics, we now investigate the bouncing restitution coefficient. To
predict the restitution coefficient, we analyse the energy budget during
bouncing. The only source of dissipation comes from the liquid viscosity, such
that the energy dissipation rate is given by

$\frac{\mathrm{d}E}{\mathrm{d}t}=-\int_{\mathbb{R}^{3}}\eta(\nabla
v)^{2}\,\mathrm{d}^{3}x=-\int_{\mathbb{R}^{2}}\frac{h^{3}}{12\eta}(\nabla
p)^{2}\,\mathrm{d}^{2}x,$ (22)

where the total viscous dissipation in Eq. (22) is simplified to keep the
dominant term within the lubrication approximation (see appendix C in Ref.
[46]). The total energy of the system, including kinetic energy and elastic
energy of the sphere is denoted $E$ and is shown in Fig. 8(a). Interestingly,
the amount of energy dissipated in each phase of the dynamics is of the same
order of magnitude, which motivates a model accounting for all the phases of
the bouncing dynamics. We define the energy dissipated in each phase of the
dynamics by $\Delta E_{i}=E(t_{i})-E(t_{i-1})$, for $i\in[1,5]$ and where
$t_{0}=0$, and aim at describing the asymptotic behaviour of each phase in the
large $\mathrm{St}$ limit.

First, during the approach phase, the elastic deformations are small (see Fig.
2(a)), such that the lubrication force can be approximated by the lubrication
force of rigid planar surface $F(t)=6\pi\eta R^{2}V_{0}/D(t)$. Here we suppose
that the velocity is asymptotically equal to the impact velocity $V(t)\simeq-
V_{0}$, such that the distance decreases linearly in time, as $D(t)\simeq
H-V_{0}t$. Here, we generalize the result to an arbitrary initial height
(denoted $H$), which is $D_{0}$ in the numerical simulations of Fig. 8.
Without elastic deformations, the energy dissipation rate can be expressed by
$F(t)V(t)=-6\pi\eta R^{2}V_{0}^{2}/D(t)$. Integrating the dissipation rate, we
obtain the amount of energy dissipated as

$\Delta E_{1}=6\pi\eta
R^{2}V_{0}\ln\left(\frac{H}{\delta}\right)=mV_{0}^{2}\mathrm{St}^{-1}\ln\left(\frac{H}{\delta}\right),$
(23)

where one needs to introduce a cut-off length $\delta$, a priori unknown, to
regularize the integral. Similar expressions have already been derived
previously [25]. As shown in Appendix C, the typical film thickness scale at
the transition between the approach and spreading phase is given by
$D_{0}\mathrm{St}^{-2/5}$, which is a good candidate of regularizing length.
Hence, we fit the numerical energy drop in the approach phase with an
asymptotic law
$mV_{0}^{2}\mathrm{St}^{-1}\ln\left(A_{1}\mathrm{St}^{2/5}\right)$, where
$A_{1}$ is a fitting constant. An excellent agreement is found with the
numerical data (see Fig. 8(b)), where the fitting parameter is $A_{1}=1.71$.
The same arguments can be employed for the bouncing phase, such that the
energy drop should also follow the asymptotic law Eq. (23), but with a
different cut-off length. Fitting the numerical value of $\Delta E_{5}$ with
$mV_{0}^{2}\mathrm{St}^{-1}\ln\left(A_{5}\mathrm{St}^{2/5}\right)$ gives an
excellent agreement (see Fig. 8(f)), where $A_{5}=0.67$.

In the spreading and retraction phase, the liquid flux is maximum in the neck
region, such that the majority of the viscous dissipation occurs in this
region. Estimating the viscous dissipation rate in Eq. (22) as
$h_{n,r}^{3}(p_{n,r}/\ell_{n,r})^{2}(a\ell_{n,r})/\eta$, and using the scales
identified in sections IV.1 and V, the energy drop is asymptotically

$\frac{\Delta
E_{2}}{mV_{0}^{2}}=A_{2}\,\mathrm{St}^{-4/5},\quad\quad\frac{\Delta
E_{3}}{mV_{0}^{2}}=A_{3}\,\mathrm{St}^{-1/2},$ (24)

where $A_{2,3}$ are numerical constants. In this case, the fitting of the
numerical energy drop in the range of Stokes number explored numerically (i.e.
$\mathrm{St}\in[10,10^{4}]$) does not provide a very good fit. To get an
approximate expression more accurate to describe the experimentally relevant
range (see Fig. 1(b)), we fit the numerical energy drop in the range
$\mathrm{St}\in[10^{2},10^{3}]$, leaving the exponent as a fitting parameter.
A very good agreement is found by using

$\frac{\Delta
E_{2}}{mV_{0}^{2}}=1.26\,\mathrm{St}^{-0.676},\quad\quad\frac{\Delta
E_{3}}{mV_{0}^{2}}=0.358\,\mathrm{St}^{-0.556},$ (25)

as shown in Fig. 8. The expression Eq. (25) should not be seen as an exact
asymptotic expression, but as a power-law fit in the Stokes range
$\mathrm{St}\in[10^{2},10^{3}]$. Indeed, the energy drop would be well
described by the exponent predicted by the elastohydrodynamic soft slider
solution ($\mathrm{St}^{-4/5}$, see Eq. (24)) only when the elastohydrodynamic
parameter $\lambda^{1/5}\propto\mathrm{St}^{-1/5}$ (see Eq. (11)) is small.
Such a condition is achieved for $\mathrm{St}\gg 10^{5}$, which is outside of
the range explored numerically. Additional simulations at larger Stokes number
would allow to determine the asymptotic energy drop during contact.

Finally, the viscous adhesion phase has scales that are inertia-free, i.e.
independent of the mass of the sphere. Using the elastohydrodynamic scales
(see Appendix C), we expect an asymptotic energy drop of the form

$\frac{\Delta E_{4}}{mV_{0}^{2}}=A_{4}\,\mathrm{St}^{-1},$ (26)

that is indeed mass-free. Fitting the numerical solutions of $\Delta E_{4}$
with (26) gives a very good agreement with $A_{4}=3.886$ as a fitting
parameter. Then, the global energy budget reads
$m(V_{0}^{2}-V_{\infty}^{2})/2=\sum_{i=1}^{5}\Delta E_{i}$, which provides an
expression for the restitution coefficient as

$\frac{V_{\infty}}{V_{0}}=\sqrt{1-2\left[\mathrm{St}^{-1}\ln\left(1.71\,\mathrm{St}^{2/5}\right)+1.26\,\mathrm{St}^{-0.676}+0.358\,\mathrm{St}^{-0.556},+3.886\,\mathrm{St}^{-1}+\mathrm{St}^{-1}\ln\left(0.67\,\mathrm{St}^{2/5}\right)\right]}.$
(27)

The asymptotic prediction of the restitution coefficient with the numerical
data is excellent, as shown in Fig. 1(b), which is not surprising as it relies
on successive fitting of the energy drop through the dynamics. Most
importantly, it provides an asymptotic expression of the restitution
coefficient at large Stokes number, highlighting the physics at each phase of
the bouncing. In the experimental range of Stokes number (see Fig. 1(b)), the
full expression is necessary to obtain a good approximation of the energy
dissipated during the impact, as the energy drop in each phase is of similar
magnitude (see Fig. 8(a)).

##  VII Conclusion

In this article, we have performed numerical simulations and asymptotic
analysis of the elastohydrodynamic bouncing of a soft sphere on a rigid
surface. We have demonstrated that the lubricated film thickness has non-
trivial self-similar dynamics, that are analogous to the steady problem of a
soft sliding sphere. Interestingly, the typical scales of the lubricated film
are different in the spreading and retraction phases. The characterization of
the self-similar features of the lubrication layer allows us to find an
asymptotic expression of the restitution coefficient in elastohydrodynamic
bouncing.

More generally, this article provides a general framework to study the
coupling between the lubrication layer and interface deformations during
impacts or collisions. Our model focuses on one of the most simple system, but
many effects important in real impacts could be implemented, e.g. surface
roughness [47], viscoelasticity [48], adhesion [49], etc… An interesting
extension of this work would be to consider oblique collisions [50, 21] and
investigating the torque generated by the shear forces during the contact.
Additionally, the present framework could be extended to a large variety of
systems either changing the impactor, e.g. drop impacts [51], elastic capsules
[52, 53], or changing the impacted surface, e.g., stretched membranes [54, 55,
56], liquid pool [57, 58, 59] and so on.

Supplementary movie. A supplementary movie is available at

Acknowledgements. It is a pleasure to thank P. Chantelot, J. Eggers, A. Oratis
and J. H. Snoeijer for stimulating discussions, as well as D. Bonn, B. Gorin
and N. Ribe for sharing unpublished experimental results that motivate this
work. I also thank the anonymous referees for constructive feedback that
greatly improved this paper.

Founding.

We acknowledge financial support from NWO through the VICI Grant No.
680-47-632.

Declaration of interests.

The author reports no conflict of interest.

Data availability statement. The code developed in the present article is
permanently available on GitHub [60].

Author ORCIDs.

Vincent Bertin: https://orcid.org/0000-0002-3139-8846

##  Appendix A Dimensionless equations

The appendix expands the non-dimensionalization of the soft lubrication
equations (1)-(3). We choose to use the typical elastic scales to make the
equations dimensionless, which does not account for viscous effects. The
typical elastic deformation of a sphere with an impact speed $V_{0}$ follows

$D_{0}=\left(\frac{mV_{0}^{2}}{E^{*}\sqrt{R}}\right)^{2/5},$ (28)

obtained by balancing the kinetic energy $mV_{0}^{2}$ with the elastic energy
$E^{*}\sqrt{R}D_{0}^{5/2}$. Using the Hertz pressure profile, we introduce a
pressure scale $p_{0}$ based on the typical elastic deformation as

$p_{0}=\frac{2E^{*}}{\pi}\sqrt{\frac{D_{0}}{R}}=\frac{2E^{*}}{\pi}\left(\frac{mV_{0}^{2}}{E^{*}R^{3}}\right)^{1/5}.$
(29)

Hence, we introduce dimensionless variables with tilde as

$\begin{split}h(r,t)=D_{0}\tilde{h}(\tilde{r},&\tilde{t}),\quad\quad
w(r,t)=D_{0}\tilde{w}(\tilde{r},\tilde{t}),\quad\quad
p(r,t)=p_{0}\tilde{p}(\tilde{r},\tilde{t}),\quad\quad
F(t)=p_{0}RD_{0}\tilde{F}(\tilde{t}),\\\
D(t)&=D_{0}\tilde{D}(\tilde{t})\quad\quad
V(t)=V_{0}\tilde{V}(\tilde{t})\quad\quad r=\sqrt{RD_{0}}\tilde{r},\quad\quad
t=\frac{D_{0}}{V_{0}}\tilde{t}.\end{split}$ (30)

Substituting Eq. (30) into (1)-(3), we obtain the system of equations

$\frac{\partial\tilde{h}(\tilde{r},\tilde{t})}{\partial\tilde{t}}=\frac{\mathrm{St}}{\tilde{r}}\frac{\partial}{\partial\tilde{r}}\left(\tilde{r}\tilde{h}^{3}(\tilde{r},\tilde{t})\frac{\partial\tilde{p}(\tilde{r},\tilde{t})}{\partial\tilde{r}}\right),$
(31a)
$\tilde{h}(\tilde{r},\tilde{t})=\tilde{D}(\tilde{t})+\frac{\tilde{r}^{2}}{2}-\tilde{w}(\tilde{r},\tilde{t}),$
(31b)
$\tilde{w}(\tilde{r},\tilde{t})=-\frac{8}{\pi^{2}}\int_{0}^{\infty}\mathcal{M}(\tilde{r},\tilde{x})\tilde{p}(\tilde{x},\tilde{t})\,\mathrm{d}\tilde{x},\quad\quad\mathcal{M}(\tilde{r},\tilde{x})=\frac{\tilde{x}}{\tilde{r}+\tilde{x}}K\left(\frac{4\tilde{r}\tilde{x}}{(\tilde{r}+\tilde{x})^{2}}\right),$
(31c)
$\frac{\mathrm{d}\tilde{V}}{\mathrm{d}\tilde{t}}=\tilde{F}(\tilde{t}),\quad\quad\tilde{F}(\tilde{t})=4\int_{0}^{\infty}\tilde{p}(\tilde{r},\tilde{t})\,\tilde{r}\mathrm{d}\tilde{r},$
(31d) $\tilde{V}=\frac{\mathrm{d}\tilde{D}}{\mathrm{d}\tilde{t}}.$ (31e)

The sphere is assumed to start with an impact velocity $V_{0}$ and at a
distance $D_{0}$ for the surface, leading to the initial conditions

$\tilde{D}(0)=1,\quad\quad\tilde{V}(0)=-1.$ (32)

Initial conditions are also required for the pressure and deformation fields.
A naive choice would be to set the initial pressure and deformation to zero,
which leads to unphysical jumps at $t=0^{+}$ that affect the approach phase.
To bypass this issue, we introduce an initialization phase where the velocity
is ramped from $0$ to $-V_{0}$ during a time $t_{\mathrm{ini}}$, as suggested
in Ref. [34]. More precisely, the velocity and sphere position are set to

$\tilde{V}(\tilde{t})=-\left(1+\frac{\tilde{t}}{\tilde{t}_{\mathrm{ini}}}\right),\quad\quad\tilde{D}(\tilde{t})=1-\left(\tilde{t}+\frac{\tilde{t}^{2}}{2\tilde{t}_{\mathrm{ini}}}\right),\quad\quad-\tilde{t}_{\mathrm{ini}}\leq\tilde{t}<0,$
(33)

and the deformation and pressure fields are identically zero at
$\tilde{t}=-\tilde{t}_{\mathrm{ini}}$. Then, solving the initialization phase
with the numerical scheme discussed in Appendix B allows to obtain initial
pressure and deformation fields $\tilde{p}(\tilde{r},0)$ and
$\tilde{w}(\tilde{r},0)$. We have checked that the choice of the
initialization time has negligible effects in the overall bouncing dynamics.

##  Appendix B Finite-difference scheme

In this appendix, we detail the finite-difference scheme used to solve the
dimensionless lubrication equations (31). We follow the methodology introduced
by Liu et al. [34]. We introduce a uniform spatial grid
$\tilde{r}_{i}=i\,\Delta r$, for $i\in[0,N-1]$, where $\Delta r$ is the grid
size and $N$ the number of radial points. The temporal axis is also
discretized by using a constant time step $\Delta t$ as $\tilde{t}^{n}=n\Delta
t$. The time step is set to $10^{-4}$ or below, which is small enough to avoid
any negative values of $h$ and which ensures numerical stability. Hence, the
pressure, film thickness and deformation fields are discretized as
$\tilde{p}(\tilde{r},\tilde{t})=\tilde{p}_{i}^{n}$,
$\tilde{h}(\tilde{r},\tilde{t})=\tilde{h}_{i}^{n}$ and
$\tilde{w}(\tilde{r},\tilde{t})=\tilde{w}_{i}^{n}$ and the velocity and sphere
position $\tilde{V}=\tilde{V}^{n}$ and $\tilde{D}=\tilde{D}^{n}$ respectively.
The dimensionless thin-film equation (31a) can be expanded as

$\tilde{V}-\frac{\partial\tilde{w}(\tilde{r},\tilde{t})}{\partial
t}=\mathrm{St}\left[\frac{\tilde{h}^{3}}{\tilde{r}}\frac{\partial\tilde{p}}{\partial\tilde{r}}+3\tilde{h}^{2}\left(r-\frac{\partial\tilde{w}}{\partial\tilde{r}}\right)\frac{\partial\tilde{p}}{\partial\tilde{r}}+\tilde{h}^{3}\frac{\partial^{2}\tilde{p}}{\partial\tilde{r}^{2}}\right],$
(34)

where we used the film thickness definition Eq. (31b). The non linear term
$h^{3}$ of the thin-film equation is evaluated at the previous time, so that
we end up with a linear discrete set of equations, which greatly reduces the
computational time as compared to non linear schemes. Furthermore, we use an
implicit scheme which provides a better numerical stability as compare to
explicit methods. Hence, Eq. (34) is discretized as

$\tilde{V}^{n+1}-\frac{\tilde{w}^{n+1}_{i}-\tilde{w}_{i}^{n}}{\Delta
t}=\mathrm{St}\left[\frac{(\tilde{h}^{n}_{i})^{3}}{\tilde{r}_{i}}\left.\frac{\partial\tilde{p}}{\partial\tilde{r}}\right|^{n+1}_{i}+3(\tilde{h}^{n}_{i})^{2}\left(r_{i}-\left.\frac{\partial\tilde{w}}{\partial\tilde{r}}\right|^{n}_{i}\right)\left.\frac{\partial\tilde{p}}{\partial\tilde{r}}\right|^{n+1}_{i}+(\tilde{h}^{n}_{i})^{3}\left.\frac{\partial^{2}\tilde{p}}{\partial\tilde{r}^{2}}\right|^{n+1}_{i}\right],$
(35)

where the film thickness discretization is
$\tilde{h}^{n}_{i}=\tilde{D}^{n}+\tilde{r}_{i}^{2}/2-\tilde{w}_{i}^{n}$. The
first and second order discrete spatial derivative of Eq. (35) are evaluated
as

$\left.\frac{\partial\tilde{w}}{\partial\tilde{r}}\right|^{n}_{i}=\frac{\tilde{w}_{i+1}^{n}-\tilde{w}_{i}^{n}}{\Delta
r},\quad\quad\left.\frac{\partial\tilde{p}}{\partial\tilde{r}}\right|^{n+1}_{i}=\frac{\tilde{p}_{i+1}^{n+1}-\tilde{p}_{i}^{n+1}}{\Delta
r},\quad\quad\left.\frac{\partial^{2}\tilde{p}}{\partial\tilde{r}^{2}}\right|^{n+1}_{i}=\frac{\tilde{p}_{i+1}^{n+1}-2\tilde{p}_{i}^{n+1}+\tilde{p}_{i-1}^{n+1}}{\Delta
r^{2}}.$ (36)

The thin-film equation is a second-order differential equation in the radial
coordinate, so that we need to introduce two boundary conditions. The symmetry
condition imposes that the gradient of the pressure in the center is null,
i.e. $\frac{\partial\tilde{p}}{\partial\tilde{r}}(\tilde{r}=0)=0$.
Furthermore, at large radius, the pressure field decays rapidly as $r^{-4}$
from Eq. (31a). Hence, we impose the boundary conditions

$p_{1}^{n}-p_{0}^{n}=0,\quad\quad p_{N-1}^{n}=0.$ (37)

To evaluate the integral equation of elasticity (31c), we suppose that the
pressure field is uniform and equal to $\tilde{p}^{n}_{i}$ on each domain
$[r_{i}-\frac{\Delta r}{2},r_{i}+\frac{\Delta r}{2}]$, for $i\geq 1$, and
$p_{0}^{n}$ on the domain $[0,\frac{\Delta r}{2}]$ near the symmetry axis.
This leads to the equation

$\tilde{w}_{i}^{n}=-\frac{8}{\pi^{2}}\left[p_{0}^{n}\int_{0}^{\frac{\Delta
r}{2}}\frac{x}{x+\tilde{r}_{i}}K\left(\frac{4x\tilde{r}_{i}}{(x+\tilde{r}_{i})^{2}}\right)\mathrm{d}x+\sum_{i=1}^{N-1}p_{i}^{n}\int_{\tilde{r}_{i}-\frac{\Delta
r}{2}}^{\tilde{r}_{i}+\frac{\Delta
r}{2}}\frac{x}{x+\tilde{r}_{i}}K\left(\frac{4x\tilde{r}_{i}}{(x+\tilde{r}_{i})^{2}}\right)\mathrm{d}x\right].$
(38)

The integrals in (38) are independent of the discrete fields, and only depends
on the spatial grid. Therefore, they can be computed once and stored in a
matrix to save some computational time. Numerically, these are evaluated with
a Gaussian quadrature using the scipy integrate library in Python. Finally,
the Newton’s second law for the sphere is discretized with a backward Euler
scheme as

$\frac{\tilde{V}^{n+1}-\tilde{V}^{n}}{\Delta
t}=4\sum_{i=0}^{N-1}\tilde{p}_{i}^{n+1}\tilde{r}_{i}\Delta
r,\quad\quad\frac{\tilde{D}^{n+1}-\tilde{D}^{n}}{\Delta t}=\tilde{V}^{n+1},$
(39)

which complete the discrete system of equations. The code is available at the
link [60].

##  Appendix C Elastohydrodynamic lubrication scales in the adhesion phase and
the transition from approach to spreading.

Figure 9:  Dimensionless minimum film thickness versus the dimensionless time.
On the left (resp. right) panel, the time is shifted by $t_{1}$ (resp.
$t_{3}$). The film thickness and time axis are rescaled by
$\mathrm{St}^{-2/5}$ corresponding to the mass-free elastohydrodynamic
lubrication scales.

In this section, we rationalize the typical scales found in the main text for
both the adhesion phase (see section III.2, and inset in Fig. 2(d)) and the
transition between approach to the spreading phase (see section IV.1). The
corresponding scales are denoted below with star in superscript. Both of these
phases occur on time scales much smaller than the bouncing time $D_{0}/V_{0}$.
Hence, the kinetic energy of the sphere does not vary significantly during
these phases and the scales must be independent of the sphere mass. The
elasticity plays a role in the approach when the elastic deformations $w^{*}$
are comparable to the film thickness $h^{*}$, i.e. $h^{*}\sim w^{*}$. The film
thickness scale can be found by balancing the typical lubrication pressure
scale $\eta V_{0}R/h^{*2}$ with the elastic pressure $E^{*}w^{*}/r^{*}$, where
$r^{*}$ is the typical radial scales. The spherical shape of the film
thickness in the contact region enforces the geometric relation
$h^{*}=r^{*2}/R$. Combining these relationships, we recover the usual
elastohydrodynamic lubrication scales as

$w^{*}=h^{*}=R\left(\frac{\eta
V_{0}}{E^{*}R}\right)^{3/5}=D_{0}\mathrm{St}^{-2/5},\quad\quad
t^{*}=\frac{h^{*}}{V_{0}}=\frac{D_{0}}{V_{0}}\mathrm{St}^{-2/5},\quad\quad
F^{*}=\frac{\eta R^{2}V_{0}}{h^{*}}=F_{0}\mathrm{St}^{-3/5}.$ (40)

The figure 9 displays the minimum film thickness versus time, for times near
the transition from approach to spreading (left) and in the adhesion (right)
phase. An excellent collapse of the numerical profiles is found using $h^{*}$
and $t^{*}$ as scale which confirm that the elastohydrodynamic lubrication
scales govern the corresponding processes.

## References

  * [1] Clark R Chapman and David Morrison.  Impacts on the earth by asteroids and comets: assessing the hazard.  Nature, 367(6458):33–40, 1994.
  * [2] Xiang Cheng, Ting-Pi Sun, and Leonardo Gordillo.  Drop impact dynamics: Impact force and stress distributions.  Annual Review of Fluid Mechanics, 54:57–81, 2022.
  * [3] Rod Cross.  The bounce of a ball.  American Journal of Physics, 67(3):222–227, 1999.
  * [4] Bruno Andreotti, Yoël Forterre, and Olivier Pouliquen.  Granular media: between fluid and solid.  Cambridge University Press, 2013.
  * [5] Kenneth Langstreth Johnson.  Contact mechanics.  Cambridge University Press, 1987.
  * [6] E Falcon, C Laroche, S Fauve, and C Coste.  Behavior of one inelastic ball bouncing repeatedly off the ground.  The European Physical Journal B-Condensed Matter and Complex Systems, 3:45–57, 1998.
  * [7] Rosa Ramírez, Thorsten Pöschel, Nikolai V Brilliantov, and Thomas Schwager.  Coefficient of restitution of colliding viscoelastic spheres.  Physical review E, 60(4):4465, 1999.
  * [8] David Tabor.  A simple theory of static and dynamic hardness.  Proceedings of the Royal Society of London. Series A. Mathematical and Physical Sciences, 192(1029):247–274, 1948.
  * [9] YM Tsai and Herbert Kolsky.  A study of the fractures produced in glass blocks by impact.  Journal of the Mechanics and Physics of Solids, 15(4):263–278, 1967\.
  * [10] Stephen C Hunter.  Energy absorbed by elastic waves during impact.  Journal of the Mechanics and Physics of Solids, 5(3):162–171, 1957\.
  * [11] J Reed.  Energy losses due to elastic wave propagation during an elastic impact.  Journal of Physics D: Applied Physics, 18(12):2329, 1985.
  * [12] MG Koller and H Kolsky.  Waves produced by the elastic impact of spheres on thick plates.  International journal of solids and structures, 23(10):1387–1400, 1987.
  * [13] Luca Brandt and Filippo Coletti.  Particle-laden turbulence: progress and perspectives.  Annual Review of Fluid Mechanics, 54:159–189, 2022.
  * [14] Jorge César Brändle De Motta, W-P Breugem, Bertrand Gazanion, J-L Estivalezes, Stéphane Vincent, and Éric Climent.  Numerical modelling of finite-size particle collisions in a viscous fluid.  Physics of Fluids, 25(8):083302, 2013.
  * [15] Pedro Costa, Bendiks Jan Boersma, Jerry Westerweel, and Wim-Paul Breugem.  Collision model for fully resolved simulations of flows laden with finite-size particles.  Physical Review E, 92(5):053012, 2015.
  * [16] GG Joseph, R Zenit, ML Hunt, and AM Rosenwinkel.  Particle–wall collisions in a viscous fluid.  Journal of Fluid Mechanics, 433:329–346, 2001.
  * [17] P Gondret, M Lance, and L Petit.  Bouncing motion of spherical particles in fluids.  Physics of fluids, 14(2):643–652, 2002.
  * [18] MMA Safa and R Gohar.  Pressure distribution under a ball impacting a thin lubricant layer.  Journal of Tribology, 108:372–376.
  * [19] Guy Barnocky and Robert H Davis.  Elastohydrodynamic collision and rebound of spheres: experimental verification.  The Physics of fluids, 31(6):1324–1329, 1988.
  * [20] M Kaneta, S Ozaki, H Nishikawa, and F Guo.  Effects of impact loads on point contact elastohydrodynamic lubrication films.  Proceedings of the Institution of Mechanical Engineers, Part J: Journal of Engineering Tribology, 221(3):271–278, 2007.
  * [21] GG Joseph and ML Hunt.  Oblique particle–wall collisions in a liquid.  Journal of Fluid Mechanics, 510:71–93, 2004.
  * [22] F-L Yang and ML Hunt.  Dynamics of particle-particle collisions in a viscous liquid.  Physics of Fluids, 18(12):121506, 2006.
  * [23] Dominique Legendre, Claude Daniel, and Pascal Guiraud.  Experimental study of a drop bouncing on a wall in a liquid.  Physics of Fluids, 17(9):097105, 2005.
  * [24] Thomas Schwager and Thorsten Pöschel.  Coefficient of normal restitution of viscous particles and cooling rate of granular gases.  Physical review E, 57(1):650, 1998.
  * [25] Robert H Davis, Jean-Marc Serayssol, and EJ Hinch.  The elastohydrodynamic collision of two spheres.  Journal of Fluid Mechanics, 163:479–497, 1986.
  * [26] Fu-Ling Yang and Melany L Hunt.  A mixed contact model for an immersed collision between two solid surfaces.  Philosophical Transactions of the Royal Society A: Mathematical, Physical and Engineering Sciences, 366(1873):2205–2218, 2008.
  * [27] Guy Barnocky and Robert H Davis.  The influence of pressure-dependent density and viscosity on the elastohydrodynamic collision and rebound of two spheres.  Journal of Fluid Mechanics, 209:501–519, 1989.
  * [28] J Wang, CH Venner, and AA Lubrecht.  Central film thickness prediction for line contacts under pure impact.  Tribology international, 66:203–207, 2013.
  * [29] Cornelis H Venner, J Wang, and AA Lubrecht.  Central film thickness in ehl point contacts under pure impact revisited.  Tribology international, 100:1–6, 2016.
  * [30] Thomas Salez and Lakshminarayanan Mahadevan.  Elastohydrodynamics of a sliding, spinning and sedimenting cylinder near a soft wall.  Journal of Fluid Mechanics, 779:181–196, 2015.
  * [31] G Lian, MJ Adams, and C Thornton.  Elastohydrodynamic collisions of solid spheres.  Journal of Fluid Mechanics, 311:141–152, 1996.
  * [32] Matthew Ryan Tan, Yumo Wang, and Joelle Frechette.  Criterion for particle rebound during wet collisions on elastic coatings.  Physical Review Fluids, 4(8):084305, 2019.
  * [33] Bhargav Rallabandi.  Fluid-elastic interactions near contact at low reynolds number.  Annual Review of Fluid Mechanics, 56:491–519, 2024.
  * [34] Zezhou Liu, Hao Dong, Anand Jagota, and Chung-Yuen Hui.  Lubricated soft normal elastic contact of a sphere: a new numerical method and experiment.  Soft Matter, 18(6):1219–1227, 2022.
  * [35] Heinrich Hertz.  The contact of elastic solids.  J Reine Angew, Math, 92:156–171, 1881.
  * [36] See Supplemental Material at [link] for Supplementary Videos SV1.
  * [37] Yumo Wang, Ziyuan Feng, and Joelle Frechette.  Dynamic adhesion due to fluid infusion.  Current Opinion in Colloid & Interface Science, 50:101397, 2020\.
  * [38] Christophe Josserand and Sigurdur T Thoroddsen.  Drop impact on a solid surface.  Annual review of fluid mechanics, 48:365–391, 2016.
  * [39] Stergios G Yiantsios and Robert H Davis.  On the buoyancy-driven motion of a drop towards a rigid surface or a deformable interface.  Journal of Fluid Mechanics, 217:547–573, 1990.
  * [40] Derek YC Chan, Evert Klaseboer, and Rogerio Manica.  Film drainage and coalescence between deformable drops and bubbles.  Soft Matter, 7(6):2235–2264, 2011.
  * [41] Jacobus Hendrikus Snoeijer, J Eggers, and Cornelis H Venner.  Similarity theory of lubricated hertzian contacts.  Physics of fluids, 25(10):101705, 2013.
  * [42] EJ Bissett.  The line contact problem of elastohydrodynamic lubrication-i. asymptotic structure for low speeds.  Proceedings of the Royal Society of London. A. Mathematical and Physical Sciences, 424(1867):393–407, 1989.
  * [43] EJ Bissett and DA Spence.  The line contact problem of elastohydrodynamic lubrication-ii. numerical solutions of the integrodifferential equations in the transition and exit layers.  Proceedings of the Royal Society of London. A. Mathematical and Physical Sciences, 424(1867):409–429, 1989.
  * [44] Jacco H Snoeijer.  Analogies between elastic and capillary interfaces.  Physical review fluids, 1(6):060506, 2016.
  * [45] Francis Patton Bretherton.  The motion of long bubbles in tubes.  Journal of Fluid Mechanics, 10(2):166–188, 1961.
  * [46] Vincent Bertin, Carmen L Lee, Thomas Salez, Elie Raphaël, and Kari Dalnoki-Veress.  Capillary levelling of immiscible bilayer films.  Journal of Fluid Mechanics, 911:A13, 2021.
  * [47] Yumo Wang and Joelle Frechette.  Morphology of soft and rough contact via fluid drainage.  Soft Matter, 14(37):7605–7614, 2018.
  * [48] Anupam Pandey, Stefan Karpitschka, Cornelis H Venner, and Jacco H Snoeijer.  Lubrication of soft viscoelastic solids.  Journal of fluid mechanics, 799:433–447, 2016.
  * [49] MP Keh and LG Leal.  Adhesion and detachment of a capsule in axisymmetric flow.  Physical Review Fluids, 1(1):013201, 2016.
  * [50] R. D. Mindlin and H. Deresiewicz.  Elastic Spheres in Contact Under Varying Oblique Forces.  Journal of Applied Mechanics, 20(3):327–344, 06 2021.
  * [51] Shreyas Mandre, Madhav Mani, and Michael P Brenner.  Precursors to splashing of liquid droplets on a solid surface.  Physical review letters, 102(13):134502, 2009.
  * [52] Etienne Jambon-Puillet, Trevor J Jones, and P-T Brun.  Deformation and bursting of elastic capsules impacting a rigid wall.  Nature Physics, 16(5):585–589, 2020.
  * [53] Théophile Rémond, Vincent Dolique, Franck Vittoz, Sheedev Antony, Renaud G Rinaldi, Lionel Manin, and Jean-Christophe Géminard.  Dynamical buckling of a table-tennis ball impinging normally on a rigid target: Experimental and numerical studies.  Physical Review E, 106(1):014207, 2022.
  * [54] Laurent Courbin, Antonin Marchand, Ashkan Vaziri, Armand Ajdari, and Howard A Stone.  Impact dynamics for elastic membranes.  Physical review letters, 97(24):244301, 2006.
  * [55] R Verzicco and G Querzoli.  On the collision of a rigid sphere with a deformable membrane in a viscous fluid.  Journal of fluid mechanics, 914:A19, 2021.
  * [56] Elvis A Agüero, Luke Alventosa, Daniel M Harris, and Carlos A Galeano-Rios.  Impact of a rigid sphere onto an elastic membrane.  Proceedings of the Royal Society A, 478(2266):20220340, 2022.
  * [57] Carlos A Galeano-Rios, Radu Cimpeanu, Isabelle A Bauman, Annika MacEwen, Paul A Milewski, and Daniel M Harris.  Capillary-scale solid rebounds: experiments, modelling and simulations.  Journal of Fluid Mechanics, 912:A17, 2021.
  * [58] Thomas C Sykes, Radu Cimpeanu, Ben D Fudge, J Rafael Castrejón-Pita, and Alfonso A Castrejón-Pita.  Droplet impact dynamics on shallow pools.  Journal of Fluid Mechanics, 970:A34, 2023.
  * [59] Luke FL Alventosa, Radu Cimpeanu, and Daniel M Harris.  Inertio-capillary rebound of a droplet impacting a fluid bath.  Journal of Fluid Mechanics, 958:A24, 2023.
  * [60] Code available at https://github.com/vincent-bertin/elastohydrodynamic-bouncing.