Similarity solutions in elastohydrodynamic bouncing

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.


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 ∞ 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 St = mV 0 /(6πηR 2 ), where m, V 0 , η 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 ∼ 1000.An empirical expression V ∞ /V 0 = exp(−35/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/δ), where H is the initial separation distance and δ 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].

D(t)
< l a t e x i t s h a 1 _ b a s e 6 4 = " g x Y V O U W K M y Q e W u p y 2 C g q B z 7 U 3 Y Q = " > A A A C x 3 i c j V H L S s N A F D 3 G V 6 2 v q k s 3 w S r U T U m q o M u C L n R X w T 6 g F k n S a R u a s G 4 0 0 t e v e e i b + Z j I 1 q / c s z 8 3 w r m 9 J A 3 Z / j n M e t G t V 9 6 h a a x 5 X 6 n Y + 6 i L 2 s I 9 D m u c J 6 r h A A y 3 j / Y g n P F v n l r C k l X 2 m W o V c s 4 t v y 3 r 4 A G b w j 2 E = < / l a t e x i t > D(t) 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.
The morphology of the lubrication layer during the contact has rarely been addressed [31].The central film thickness has been identified to scale as 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.

A. 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 η.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 V (t) = F (t), where F (t) is the the vertical hydrodynamic viscous force acting on the particle, and V (t) = Ḋ(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 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) = ∞ 0 p(r, t) 2πrdr.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 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 where ) is the plane strain elastic modulus with E and ν 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, obtained by balancing the kinetic energy mV 2 0 with the Hertz elastic energy E * √ RD 5/2 0 (see below).Equations ( 1)-( 3) are made dimensionless, by using the scales of the Hertz theory, respectively D 0 , √ RD 0 , D 0 /V 0 , p 0 = 2E * D 0 /R/π 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 The Stokes number can be interpreted as the ratio between sphere inertia (mV 2 0 ) to viscous dissipation (6πηR 2 V 0 ), or alternatively the viscous dissipation time scale mD 0 /(6πη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 ∆ρR 3 gD 0 is small compare to the kinetic energy mV 2 0 [6], where ∆ρ 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 2 0 /(E * R 3 )) 2/5 ≪ 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 A. Dry limit, St → ∞: Hertz theory
In the absence of any surrounding fluid, or for St → ∞, 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 δ 3/2 (t)/3, where we introduce δ(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 which corresponds to a nonlinear mass-spring equation, where the spring stiffness increases with the indentation as √ δ.An exact solution of Eq. ( 6) can be expressed in the form of an implicit equation where 2 F 1 is the Hypergeometric function.We notice that the equivalent dynamical law for a linear spring model is arcsin(δ(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 where a(t) = δ(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 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 All these results will turn out relevant as "outer solutions" for the impact dynamics at finite Stokes number.

B. Finite Stokes number
The bouncing dynamics at finite Stokes number is illustrated in Fig. 2 (a) (see also Supp.Video 1 [36]) for 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 ∈ [1,5], that separate the phases.Initially, the sphere approaches the surface and the elastic deformation remains small, i.e., of order 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 ∞ = 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 St = 50 and 200.The specific times t i , for i ∈ [1,5], are indicated with circles only for the case 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 ∞ /V 0 < 1, due to the presence of dissipation.Figure 1 (b) displays the coefficient of restitution V ∞ /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 St −3/5 and time by 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 ηRV 0 [ηV 0 /(RE * )] −2/5 and R/V 0 [ηV 0 /(RE * )] 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 dp/dr 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).

A. Neck region: soft lubrication inlet analogy
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 → 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 = ȧ and under a load 4E * a 3 /(3R).In Refs.[41][42][43], 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 , ℓ 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 η∇ 2 u = ∇p, in the horizontal direction yields to ηu/h 2 n = p n /ℓ 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 /ℓ 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 ℓ 3/2 n /R.Combining the aforementioned expressions, one finds The dimensionless group λ is the relevant elastohydrodynamic dimensionless number, which needs to be small to enforce the hierarchy of scales h n ≪ ℓ n ≪ a ≪ R of the asymptotic expansion.Following this approach, we introduce a self-similarity ansatz of the form Substituting Eq. ( 12) into Eq.( 1), and using λ ≪ 1, one gets 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 ≈ −V 0 , such that the contact radius dynamics is governed by a(t) ≈ RV 0 (t − t 1 ) and u(t) ≈ RV 0 /[4(t − t 1 )].Hence, the elastohydrodynamic parameter, which scales as λ ∝ (u/a 4 ) ∝ (t − t 1 ) −5/2 , diverges as t → t 1 .The steady similarity solution develops after a transient regime, once λ ≪ 1.The typical scale of the transient time can be estimated as the one when λ = 1, leading to (t − t 1 ) ∼ (D 0 /V 0 )St −2/5 , which is much smaller than the bouncing time scale D 0 /V 0 in the large 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 where H 0 is an integration constant.In the limit where the length of the neck is smaller than the contact radius, i.e. ℓ n ≪ a, the elastic kernel M of Eq. ( 3) reduces to the dimensionless line force Green's function M(r, r ′ ) ∝ − 1 2 ln | ξ − ξ ′ | [5].Taking (twice) the derivative of Eq. ( 2), and combining with the line force elastic kernel, we find where we assume ℓ 2 n /(h n R) ≪ 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 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.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   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 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.
varying the Stokes number -in particular since we expect the neck to vanish in the limit St → ∞.In figure 4(a), we thus plot profiles at different 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 Hence, we find that the Stokes number scaling of the film thickness and lateral scale of the neck during the spreading phase are St −3/5 and 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 St −3/5 over the full St number range explored numerically.

B. Central region: dimple height model
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 → 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 St −1/2 (see Fig. 5(b)).In Fig. 4(c), we show the radial dependence of the film thickness profiles rescaled by St −1/2 , where the profiles indeed collapse near the center.Nevertheless, the profiles deviate away from the symmetry axis.Therefore, the D 0 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. (8a) (see Fig. 3(c)).Injecting the Hertzian pressure field in the thin film equation Eq. (1) and investigating the limit r → 0, we obtain an ordinary differential equation for h 0 (t) as ∂h(r, t) ∂t = − E * 3πηR ) 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.

Solving this equation with the initial condition h
Not only we recover the 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 ≪ D 0 /V 0 , the Hertz contact radius follows a(t) = RV 0 (t − t 1 ).In this limit, the dimple height expression simplifies to 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
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 H 0 h n that is universal and selected by the inlet profile, where H 0 = 2.478... is a numerical prefactor [41].For the soft slider, the very same scaling arguments of section IV A 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 St −3/5 scaling as in the spreading phase, but is closer to the 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 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 St scaling.For 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 ∼ D 0 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 St −1/2 , where the subscript r stands for retraction.As in section IV A, 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 = ℓ 3/2 r a 1/2 /R and p r = E * a 1/2 ℓ 1/2 r /R.Combining the latter expressions, we find the scales 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.

VI. ENERGY DISSIPATION AND COEFFICIENT OF RESTITUTION
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 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 ∆E i = E(t i ) − E(t i−1 ), for i ∈ [1,5] and where t 0 = 0, and aim at describing the asymptotic behaviour of each phase in the large 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πηR 2 V 0 /D(t).Here we suppose that the velocity is asymptotically equal to the impact velocity V (t) ≃ −V 0 , such that the distance decreases linearly in time, as D(t) ≃ 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 Integrating the dissipation rate, we obtain the amount of energy dissipated as where one needs to introduce a cut-off length δ, 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 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 2 0 St −1 ln A 1 St 2/5 , 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 ∆E 5 with mV 2 0 St −1 ln A 5 St 2/5 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 3 n,r (p n,r /ℓ n,r ) 2 (aℓ n,r )/η, and using the scales identified in sections IV A and V, the energy drop is asymptotically where A 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 St ∈ [10 2 , 10 3 ].Indeed, the energy drop would be well described by the exponent predicted by the elastohydrodynamic soft slider solution (St −4/5 , see Eq. ( 24)) only when the elastohydrodynamic parameter λ 1/5 ∝ St −1/5 (see Eq. ( 11)) is small.Such a condition is achieved for St ≫ 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 that is indeed mass-free.Fitting the numerical solutions of ∆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 2 0 −V 2 ∞ )/2 = 5 i=1 ∆E i , which provides an expression for the restitution coefficient as 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 selfsimilar 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 [21,50] 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
Appendix B: Finite-difference scheme In this appendix, we detail the finite-difference scheme used to solve the dimensionless lubrication equations (A4).We follow the methodology introduced by Liu et al. [34].We introduce a uniform spatial grid ri = i ∆r, for i ∈ [0, N − 1], where ∆r is the grid size and N the number of radial points.The temporal axis is also discretized by using a constant time step ∆t as tn = n∆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 p(r, t) = pn i , h(r, t) = hn i and w(r, t) = wn i and the velocity and sphere position Ṽ = Ṽ n and D = Dn respectively.The dimensionless thin-film equation (A4a) can be expanded as where we used the film thickness definition Eq. (A4b).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. (B1) is discretized as where the film thickness discretization is hn i = Dn + r2 i /2 − wn i .The first and second order discrete spatial derivative of Eq. (B2) are evaluated as 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. ∂ p ∂ r (r = 0) = 0. Furthermore, at large radius, the pressure field decays rapidly as r −4 from Eq. (A4a).Hence, we impose the boundary conditions (B5) The integrals in (B5) 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 Appendix C: Elastohydrodynamic lubrication scales in the adhesion phase and the transition from approach to spreading.
In this section, we rationalize the typical scales found in the main text for both the adhesion phase (see section III B, and inset in Fig. 2(d)) and the transition between approach to the spreading phase (see section IV A).The corresponding scales are denoted below with star in superscript. 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 * ∼ w * .The film thickness scale can be found by balancing the typical lubrication pressure scale η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 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.

F 5 N
J s R Q X / o B b / T P x D / Q v v D N O Q S 2 i E 5 K c O f e e M 3 P v d Z P A T 4 V l v c 4 Z 8 w u L S 8 u 5 l f z q 2 v r G Z m F r u 5 H G G f d Y 3 Y u D m L d c J 2 W B H 7 G 6 8 E X A W g l n T u g G r O k O z 2 S 8 O W I 8 9 e P o W o w T 1 g m d f u T 3 f M 8 R k j o v i c P b Q t E q W 2 q Z s 8 D W o A i 9 a n H h B T f o I o a H D C E Y I g j C A R y k 9 L R h w 0 J C X A c T 4 j g h X 8 U Z 7 p E n b U Z Z j D I c Y o f 0 7 d O ur d m I 9 t I z V W q P T g n o 5 a Q 0 c U C a m P I 4 Y X m a q e K Z c p b s b 9 4 T 5 S n v N q a / q 7 1 C Y g U G x P 6 l m 2 b +

FIG. 2 .
FIG. 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 V0t/D0 = 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 St = 50 and 200.The black dashed lines indicate the Hertz theory (7), corresponding to the absence of viscous dissipation (St = ∞).The points illustrate the five characteristic times separating the different phases for the case 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 St −2/5 and St −3/5 respectively.We also point out that the origin of time has been shifted in the inset by t3 where the force vanishes.

FIG. 3 .
FIG. 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 D0/V0) 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.

FIG. 4 .
FIG. 4. Spreading phase: Stokes number scaling.(a) Typical dimensionless film thickness as a function of the dimensionless radius at t = 1.5D0/V0 during the spreading phase.The three colors indicate different Stokes number, respectively 50, 200, 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.
t e x i t s h a 1 _ b a s e 6 4 = " C y 8 w O h t l 4 n 8 Q / a t J p3 i Z D S m w v o 8 = " > A A A C x n i c j V H L T s J A F D 3 U F + I L d e m m k Z i 4 I m 1 p A H d E N y w x y i N B Q t o y Q E N f a a c a Q k z 8 A b f 6 a c Y / 0 L / w z l g S XR C d p u 2 d c 8 8 5 M / d e O / L c h G v a e 0 5 Z W 9 / Y 3 M p v F 3 Z 2 9 / Y P i o d H n S R M Y 4 e 1 n d A L 4 5 5 t J c x z A 9 b m L v d Y L 4 q Z 5 d s e 6 9 q z K 5 H v 3 r M 4 c c P g l s 8 j N v C t S e C O X c FIG. 5.Dimple height.(a) Rescaled central film thickness height h0(t) = h(r = 0, t) as 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 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.

FIG. 6 .
FIG.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 D0/V0) 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.

FIG. 7 .
FIG. 7.Retraction phase: number scaling.(a) Typical dimensionless film thickness as a function of the dimensionless radius at t = 3.2D0/V0 during the spreading phase.The four colors indicate different Stokes number, respectively 50, 200, 1000 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.

FIG. 8 .
FIG. 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 ∆Ei, for i ∈ [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.
integral equation of elasticity (A4c), we suppose that the pressure field is uniform and equal to pn i on each domain [r i − ∆r 2 , r i + ∆r 2 ], for i ≥ 1, and p n 0 on the domain [0, ∆r 2 ] near the symmetry axis.This leads to the equation wn ri ) 2 dx .

Ṽ n+1 − Ṽ n
discrete system of equations.The code is available at the link [60].

FIG. 9 .
FIG.9.Dimensionless minimum film thickness versus the dimensionless time.On the left (resp.right) panel, the time is shifted by t1 (resp.t3).The film thickness and time axis are rescaled by St −2/5 corresponding to the mass-free elastohydrodynamic lubrication scales.
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.St ∈ [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 St ∈ [10 2 , 10 3 ], leaving the exponent as a fitting parameter.A very good agreement is found by using 1 ln 1.71 St 2/5 + 1.26 St −0.676 + 0.358 St −0.556 , +3.886 St −1 + St −1 ln 0.67 St 2/5 .