1. Introduction
Droplet impacts occur frequently in both natural and industrial settings. Rain drops impacting on leaves have been shown to be a primary mechanism for pathogen transport among plants (Kim et al. Reference Kim, Park, Gruszewski, Schmale and Jung2019) and birds with superhydrophobic feathers stay warmer in a cold rain due to a reduced droplet contact time (Shiri & Bird Reference Shiri and Bird2017). Spray cooling devices have attracted the attention of researchers due to the large heat transfer rates and high uniformity of heat transfer (Kim Reference Kim2007). Wet scrubbing of exhaust gases relies on the inertial impaction of small particles and aerosols on the surface of freely falling droplets (Park et al. Reference Park, Jung, Jung, Lee and Lee2005). Various other drop impact phenomena, such as splashing, were experimentally documented by Worthington at the start of the 20th century (Worthington Reference Worthington1908). More recently, droplets bouncing repeatedly on a vertically oscillated bath have received considerable interest as a macroscopic pilot-wave system capable of reproducing some behaviours reminiscent of quantum particles (Couder et al. Reference Couder, Protiere, Fort and Boudaoud2005b; Bush & Oza Reference Bush and Oza2020). For instance, bouncing droplets confined to submerged cavities can exhibit wave-like statistical behaviour analogous to electrons in quantum corrals (Harris et al. Reference Harris, Moukhtar, Fort, Couder and Bush2013; Sáenz, Cristea-Platon & Bush Reference Sáenz, Cristea-Platon and Bush2018). Droplet impact onto solid surfaces is also an extremely well-studied field (Yarin Reference Yarin2006; Josserand & Thoroddsen Reference Josserand and Thoroddsen2016), with the combination of high quality experiments and direct numerical simulation (DNS) leading to a deep understanding of the multiscale dynamics.
The problem of droplet coalescence onto a bath of the same fluid has also been studied extensively over the last century and a half, beginning with Rayleigh (Reference Rayleigh1879) and Thomson & Newall (Reference Thomson and Newall1886). These early works included sketches of drop–interface coalescence, as well as a detailed description of the vortices that are formed in the fluid bath. Preceding coalescence, the thin gas film that forms between the two interfaces drains until van der Waals forces act to initiate coalescence. Coalescence of a drop into a bath occurs when the film is of the order of 
$100$ nm thick (Couder et al. Reference Couder, Fort, Gautier and Boudaoud2005a; Yarin Reference Yarin2006; de Ruiter et al. Reference de Ruiter, Oh, van den Ende and Mugele2012; Kavehpour Reference Kavehpour2015). The development of accessible high-speed photography and high-performance computing has ushered in a rapid expansion in quantity and quality of data on these free surface problems. Thoroddsen & Takehara (Reference Thoroddsen and Takehara2000) used a high-precision and high-speed visualization set-up to quantify the coalescence time of droplets on an air–liquid interface. Tang et al. (Reference Tang, Saha, Law and Sun2019) studied the dynamics of the gas layer on a liquid bath whose depth was similar to that of the droplet radius. The rich class of outcomes and dynamics that arise from such a simple interaction between droplet, surrounding gas and interface proves that these fundamental problems merit considerable attention.
During contact, the combined effects of inertia, surface tension, gravity and viscosity govern the hydrodynamic interaction between the droplet and the interface, and the complex balance of forces within this regime creates a variety of distinct phenomena. The Weber number 
${We}= \rho V_0^2 R /\sigma$, the Bond number 
${Bo}= \rho R^2 g/\sigma$ and the Ohnesorge number 
${Oh} = \mu /\sqrt {\sigma R \rho }$ are often used to describe these capillary-scale dynamics. In this work, 
$R$ represents the undeformed droplet radius, 
$\rho$, 
$\sigma$, 
$\mu$ are the density, surface tension and viscosity of the fluid in both the droplet and bath, 
$V_0$ is the impact velocity of the droplet and 
$g$ is the gravitational acceleration. The present work focuses on the inertio-capillary regime, where fluid inertia and surface tension dominate viscous and gravitational effects (specifically, 
${Oh}\ll 1$ and 
${Bo}\ll 1$). During impact, a thin gas film develops between the free interface and surface of the droplet. The drainage of this thin film plays a crucial role in determining the fate of the droplet: specifically, whether it rebounds from or coalesces with the underlying bath (de Ruiter et al. Reference de Ruiter, Oh, van den Ende and Mugele2012). At sufficiently low 
${We}$, the droplet and the interface never come into physical contact and remain separated by a stable air film. The droplet then levitates on this thin film and can eventually rebound due to the relaxation of the bath and droplet interface. Droplets bouncing on a free interface were first documented by Reynolds in 1881, when he noted that drops can ‘float’ on a bath of the same liquid if the impact velocity is sufficiently small (Reynolds Reference Reynolds1881).
In cases of droplet–bath impact, as well as droplet–droplet impact, there exists a parameter regime where droplets bounce completely (figure 1). The bouncing–coalescence threshold is often characterized by a critical 
${We}$ that depends sensitively on all parameters in the problem (Tang et al. Reference Tang, Saha, Law and Sun2019). Droplet bouncing on an undisturbed interface at variable impacting angles was first studied in detail by Jayaratne & Mason (Reference Jayaratne and Mason1964) experimentally. They were able to determine a relationship between the drop radius, impact speed and impact angle at which the bouncing–coalescence threshold occurs between uncharged drops. Building on the work of Gopinath & Koch (Reference Gopinath and Koch2001), Bach, Koch & Gopinath (Reference Bach, Koch and Gopinath2004) studied the droplet impact of small (
$R \leq 50\,\mathrm {\mu }\mathrm {m}$) aerosol droplets impacting a fluid bath. They developed a rarefied gas model to describe the dynamics of the gas layer separating the droplet and bath, and used an inviscid potential flow model to describe the transfer of energy from the droplet to the bath during impact. The authors determined that the criterion for drop bouncing is more sensitive to gas mean-free path and gas viscosity than to the Weber number itself. Zou et al. (Reference Zou, Wang, Zhang, Fu and Ruan2011) investigated water droplets bouncing on an air–water interface, and examined the role of bath depth in bounce-back behaviour. They determined that the contact time was independent of the impact velocity for a large range of Bond numbers. Wu et al. (Reference Wu, Hao, Lu, Xu, Hu and Floryan2020) used a drop-on-demand generator to study the bouncing of water droplets, developed a model for the maximum penetration depth, and compared it with their experimental study. They varied the droplet diameter, and found that the maximum rebound height decreased with increasing diameter. An experimental work utilizing three different fluids was completed by Zhao, Brunsvold & Munkejord (Reference Zhao, Brunsvold and Munkejord2011). They chose water, 1-propanol and ethanol as the working fluids and found good agreement in measured contact times with Jayaratne & Mason (Reference Jayaratne and Mason1964). Also, they determined that the contact time of the droplet was relatively independent of the impact velocity, similar to that found by Richard, Clanet & Quéré (Reference Richard, Clanet and Quéré2002) for a droplet impacting a non-wetting, dry surface. In the variety of experimental work on this particular problem, the scaling for contact time 
$t_c$ of the droplet appears to be mostly independent of 
${We}$, except at very low 
${We}$ (Zhao et al. Reference Zhao, Brunsvold and Munkejord2011; Zou et al. Reference Zou, Wang, Zhang, Fu and Ruan2011; Wu et al. Reference Wu, Hao, Lu, Xu, Hu and Floryan2020). Additionally, numerous papers report a saturation of translational energy recovery by the droplet at intermediate 
${We}$, as measured by the coefficient of restitution 
$\alpha$ (Jayaratne & Mason Reference Jayaratne and Mason1964; Bach et al. Reference Bach, Koch and Gopinath2004; Zhao et al. Reference Zhao, Brunsvold and Munkejord2011; Zou et al. Reference Zou, Wang, Zhang, Fu and Ruan2011; Wu et al. Reference Wu, Hao, Lu, Xu, Hu and Floryan2020). These observations have neither yet been fully explained nor their parametric dependencies clearly elucidated to the best of the current authors’ knowledge. This motivates the development of a first principles model that can accurately and efficiently describe the dynamics of both the droplet and the fluid bath over the physically relevant parameter regime.
Figure 1. A small water droplet (
$R\approx 0.4$ mm) rebounds from a bath of the same fluid.
The multiscale hydrodynamics present in these impact problems creates significant challenges for numerical simulations, and all but eliminates analytical solutions to these problems. Wagner (Reference Wagner1932) proposed the first theoretical study of an object impacting on an inviscid, incompressible fluid, utilizing linearized free-surface kinematic and dynamic conditions to develop a theory that decomposed the fluid domain into two parts, one where the applied pressure is unknown but the interface shape is known, and vice versa. The so-called Wagner theory was extended to a solid of revolution by Schmieden (Reference Schmieden1953) and eventually to three dimensions by Scolan & Korobkin (Reference Scolan and Korobkin2001). These models assume that the working fluid is ideal, and thus any waves generated upon impact are not subject to viscous dissipation. Dias, Dyachenko & Zakharov (Reference Dias, Dyachenko and Zakharov2008) derived a theory to include the effects of weak damping in free surface problems, which appear as leading-order corrections in the free surface boundary conditions. The inclusion of damping in this method provides a mechanism for the waves generated by impact to decay in time. The Dias et al. (Reference Dias, Dyachenko and Zakharov2008) theory is valid in the weakly viscous regime and represents a rigorous derivation of a linearized free surface model first proposed by Lamb (Reference Lamb1895).
More recently, Galeano-Rios, Milewski & Vanden-Broeck (Reference Galeano-Rios, Milewski and Vanden-Broeck2017), Galeano-Rios, Milewski & Vanden-Broeck (Reference Galeano-Rios, Milewski and Vanden-Broeck2019) and Galeano-Rios et al. (Reference Galeano-Rios, Cimpeanu, Bauman, MacEwen, Milewski and Harris2021) applied the quasipotential model of Dias et al. (Reference Dias, Dyachenko and Zakharov2008) to free surface impact problems and solve the problem of the unknown, time evolving contact region through the use of a so-called ‘kinematic match’. In the kinematic match framework, the free surface shape within the region of contact is determined by the geometry of the problem, and the extent of this region can be computed with the use of a tangency boundary condition. The model in Galeano-Rios et al. (Reference Galeano-Rios, Milewski and Vanden-Broeck2017) worked well in determining the trajectory of the droplet in some cases; however, it neglected any deformations of the droplet. Blanchette (Reference Blanchette2016, Reference Blanchette2017) modelled the impact of a droplet onto a still and oscillating bath, where a simplified version of the kinematic match concept was used by assuming that the shape and radial extent of the pressure distribution in the contact region are known a priori. Additionally, droplet deformations were modelled as a vertical spring or as an octahedral network of springs and masses. For still bath impacts, only very limited direct comparison with experimental measurements were made, with mixed success. Moláček & Bush (Reference Moláček and Bush2012) developed a quasistatic model for a droplet impacting on a non-wetting rigid solid surface with fixed curvature. They compared this quasistatic model with a dynamic model that described the droplet–air interface using spherical harmonics derived from a balance of surface, kinetic and potential energies and found good agreement between the two at low 
${We}$ numbers, as compared with experiments and the model of Okumura et al. (Reference Okumura, Chevy, Richard, Quéré and Clanet2003). However, these models do not predict the energy transfer and time dependent waves on a fluid bath. Terwagne et al. (Reference Terwagne, Ludewig, Vandewalle and Dorbolo2013) wrote a linear mass–spring–damper model for a bouncing droplet on a vertically oscillated bath. Similarly, this model assumed that the bath surface was non-deformable. Additionally, Moláček & Bush (Reference Moláček and Bush2013) studied silicone oil droplets bouncing on a vertically oscillated bath and developed linear and logarithmic spring models to classify bouncing dynamics. While efficient to solve, these models require the input of free parameters determined from experimental data and thus cannot independently predict bouncing metrics such as the coefficient of restitution or contact time. Other linear spring-type models have been proposed in the literature, but again such models generally rely on fitting parameters obtained from experimental data or DNS (Sanjay et al. Reference Sanjay, Lakshman, Chantelot, Snoeijer and Lohse2022a). DNSs of free surface impact problems have been completed in other recent works (Pan & Law Reference Pan and Law2007; He, Liu & Qiao Reference He, Liu and Qiao2015; Sharma & Dixit Reference Sharma and Dixit2020; Fudge, Cimpeanu & Castrejón-Pita Reference Fudge, Cimpeanu and Castrejón-Pita2021), and provide very good results, even in regimes presently inaccessible to experiments. From these simulations the droplet shape, trajectory and waves, as well as the flow within the droplet and bath, can be captured and analysed in detail. However, due to the high computational cost of these free surface flow problems, the vast parameter space encompassed by this problem renders large sweeps impractical for direct simulation, and leaves much to be understood about the overall dynamics over a more complete space.
In this work, we develop an efficient model that accurately predicts the trajectory of the impacting droplet, the instantaneous droplet shape, and the transient waves generated on the bath interface by impact, without any adjustable parameters. First, we use the Navier–Stokes equations with linearized free surface conditions and include viscosity as leading-order corrections to these boundary conditions, which holds in the limit of large 
$Re$. We then derive a set of ordinary differential equations (ODEs) to describe the motion of the bath interface. The droplet shape is modelled by another set of ODEs that govern the weakly damped oscillation of individual modes on the droplet interface that hold for small 
${Oh}$. Both the bath and droplet models are the result of linearizing about their undeformed states, and thus we anticipate best agreement when deformations are small. The bath and drop models are coupled using a single-point kinematic match condition and evolved simultaneously in time. We validate this model with new experimental data as well as DNSs. We then apply the validated model over a wide range of parameters where the relative influence of the hydrodynamic, surface tension, and gravitational forces on the rebound behaviour of the bouncing droplet will be elucidated.
2. Experimental methods
2.1. Experimental set-up
A series of droplet impact experiments were conducted utilizing two working fluids: deionized water and silicone oil with viscosity of 
$5$ cSt. Ranges of experimental parameters are summarized in table 1. A drop-on-demand generator is used to reliably produce droplets with a maximum variation in the diameter of less than 
$1\,\%$ (Ionkin & Harris Reference Ionkin and Harris2018). This device, along with a schematic of the experimental set-up, is shown in figure 2(a). The drop generator is entirely three-dimensionally printed (3D-printed), with the exception of a small piezoelectric disk, hardware and connective tubing. The deformation of the piezoelectric disk due to an applied voltage pulse acts to expel fluid through a small nozzle. As the fluid exits the nozzle, the piezoelectric disk relaxes, initiating pinch off of the droplets. The droplets then fall under the action of gravity towards the bath. Two visualizations of droplet impact and rebound are shown in figure 2(b,c). The drop generator is mounted on a 3D-printed translation stage, allowing for repeatable changes to the impacting velocity via height increases of the droplet generator. Directly underneath the drop generator is a 3D-printed fluid bath. The bath is 
$70$ mm in width and length, and 
$50$ mm deep. The impact location was 
$25$ mm from the front wall of the bath. This impact location allowed for consistent focus above and below the free surface yet was still sufficiently far from the front panel that the waves created during impact do not have time to reflect and interact with the droplet during contact. For the water experiments, the front and rear walls are constructed using polystyrene that has an equilibrium contact angle of 
$87.4^{\circ }$ (Ellison & Zisman Reference Ellison and Zisman1954). Being close to 
$90^{\circ }$, this creates a negligibly small meniscus that allows for detailed photography of the impact from the side (Galeano-Rios et al. Reference Galeano-Rios, Cimpeanu, Bauman, MacEwen, Milewski and Harris2021). For the silicone oil experiments, we use a shorter bath window panel, constructed of extruded acrylic and a thin transparent plastic sheet. The bath was brim-filled to the height of the acrylic window panel, such that the contact line was pinned with angle held at approximately 
$90^{\circ }$. The drops are imaged using a high-speed camera (Phantom Miro LC 311) and illuminated by a Phlox LED-W back light. Movie data is taken at 10 000 frames-per-second (f.p.s.) with an exposure time of 
$99.6\,\mathrm {\mu }\mathrm {s}$.
Table 1. Relevant parameters and their range of values in our experimental study.
Figure 2. (a) A rendering of the experimental set-up. (b) Experimental montage of impact of a deionized water droplet on a bath of the same fluid. Images are spaced 0.7 ms apart. (c) Spatiotemporal diagram of a deionized water droplet bouncing. The image is constructed by taking a single-pixel-wide stripe of the raw movie footage, and plotting time along the 
$x$-axis. Panels (b,c) correspond to an impact of deionized water on a bath of the same fluid with 
${We}=0.7$, 
${Bo} = 0.017$ and 
${Oh} = 0.006$. (d) Schematic of the problem.
2.2. Experimental procedure
Care must be taken to ensure that both the fluid interface and the fluid in the reservoir of the drop generator are contaminant-free, as dust or surfactants can modify the physics involved. Prior to each experiment, the bath and tubing are thoroughly cleaned with an isopropyl alcohol solution, flushed with deionized water, and then left to dry in a fume hood with particulate filtering for 30 min. The drop generator and nozzle are cleansed with an ethanol solution, and then flushed with deionized water for five minutes. Gloves are worn at all times to minimize contamination. The drop generator is controlled by an Arduino Uno board, with a simple H-bridge circuit initiating the voltage switching of a DC power supply (Ionkin & Harris Reference Ionkin and Harris2018). The fluid within the bath is periodically flushed, approximately after every 15 droplet impacts to reduce surface contamination (Kou & Saylor Reference Kou and Saylor2008). Overflow from the flushing is caught by a small lip in the bottom of the bath, which is then drained to a waste container. There are two syringes connected to the bath, which allow for fine adjustments of the equilibrium bath depth after flushing.
We collect experimental data for the top and bottom of the droplet during free flight. During contact, we track the height of the top of the droplet and that of the centre of the deformed free surface. Since the air layer that separates the droplet and the interface is negligibly thin relative to the scale of the droplet, we assume that this point is also effectively the location of the droplet's south pole. Just after the drop rebounds off the surface, the axisymmetric surface wave created by the impact is partially in the line of sight of the camera, and obscures the bottom of the droplet for a brief period during take-off. These data points have been omitted from the bottom trajectory when reported. The raw movie data are postprocessed using a custom Canny edge detection software implemented in MATLAB 2021b, which quantifies the droplet trajectory, and then computes impact parameters and bouncing metrics (Galeano-Rios et al. Reference Galeano-Rios, Cimpeanu, Bauman, MacEwen, Milewski and Harris2021).
There are several metrics of interest in our study, which we define in what follows. The maximum penetration depth, 
$\delta$, of a bounce is defined as the position of the bottom of the droplet at the lowest point in the trajectory (relative to the undisturbed interface height). In our experiment, the contact time, 
$t_c$, is defined as the time duration from which the top of the droplet crosses the height 
$z=2R$ to the time the top of the droplet returns to that height. Due to the nature of visualization set-up it was impossible to determine precisely when the droplets lost physical contact with the fluid; however, this always occurred before the top of the drop returned to the level 
$z=2R$. Each bounce was also characterized by its coefficient of restitution, 
$\alpha$, which is defined here as the negative of the normal exit velocity, 
$V_e$, divided by the normal impact velocity, 
$V_0$. This parameter ranges between 0 and 1, and is related to the momentum exchange during impact.
In order to determine the contact time (
$t_c$) and coefficient of restitution (
$\alpha$), a parabola was fitted using a least-squares method to the incoming and outgoing trajectories, separately, with at least 30 data points prior to impact and at least 40 data points following rebound. The analytical form of the parabolic fit was then used to extrapolate the time at which the sphere crosses the still air–water interface height (which corresponds to a root of the parabolic function). The derivative of the parabolic fit function was then computed analytically and its value evaluated at these times in order to calculate the impact speed, 
$V_0$, and exit speed, 
$V_e$ (Galeano-Rios et al. Reference Galeano-Rios, Cimpeanu, Bauman, MacEwen, Milewski and Harris2021).
3. Linearized quasipotential fluid model
In this section, we develop a model for droplet impact on a flat fluid interface from first principles. First, we use the linearized Navier–Stokes equations to model the flow within the bath and the shape of the bath interface. Then, we use an orthogonal function decomposition of the bath model to derive a single set of linear ODEs that govern the bath mode amplitudes. We then write a similar model for the droplet interface, which reduces to another set of linear ODEs governing droplet mode amplitudes. Finally, we propose a model for the pressure distribution and its extent acting on the bath and droplet during contact, couple the two sets of ODEs together using a single-point kinematic match condition, and solve the system implicitly using standard numerical integration techniques. A schematic of the problem is illustrated in figure 2(d).
3.1. Bath interface model
The present work models the bath interface dynamics using a linearized, quasipotential flow model following the work of Galeano-Rios et al. (Reference Galeano-Rios, Milewski and Vanden-Broeck2017). For the problem of a droplet impacting on a free interface, the Navier–Stokes equations govern the flow generated by the bath–droplet interaction. Assuming the flow to be incompressible, isothermal and Newtonian, we can define the fluid velocity vector 
$\boldsymbol {u} = [u,v,w]^{\rm T}= \boldsymbol {\nabla } \phi + \boldsymbol {\nabla } \times \varPsi$ and the bath interface shape 
$\eta =\eta (r,\varTheta,t)$, Here, 
$\phi$ is the scalar potential and 
$\varPsi$ is the vector stream function. We then linearize the governing equations and boundary conditions about the undisturbed free surface 
$z=0$. Utilizing the arguments presented in Galeano-Rios et al. (Reference Galeano-Rios, Milewski and Vanden-Broeck2017) and Dias et al. (Reference Dias, Dyachenko and Zakharov2008), we can recast the governing equations to be
$$\begin{gather} \nabla^2 \phi {+ \partial_z^2 \phi} = 0,\quad z \leq 0 , \end{gather}$$
$$\begin{gather}\partial_t \eta = \partial_z \phi + 2 \nu {\nabla^2 \eta} ,\quad z = 0 , \end{gather}$$
$$\begin{gather}\partial_t \phi ={-} g\eta -2\nu \partial_z^2 \phi + \frac{\sigma}{\rho}\kappa - \frac{p_s}{\rho} ,\quad z = 0 , \end{gather}$$
$$\begin{gather}\partial_z \phi = 0,\quad {\rm at}\ z=h_0 . \end{gather}$$
Here, 
$p_s(r,\varTheta,t)$ is the contact pressure, 
$g$ is the gravitational acceleration, 
$\kappa (r,\varTheta,t)=\nabla ^2 \eta$ is twice the linearized mean curvature of the interface, 
$h_0$ is the depth of the undisturbed bath and 
$\rho$, 
$\sigma$, 
$\nu = \mu /\rho$ are the fluid density, surface tension and kinematic viscosity, respectively. In this notation, 
$\partial _{( )}$ denotes partial differentiation with respect to the variable given in the parenthesis and 
$\nabla ^2 = \partial _{r}^2 + (1/r)\partial _r + (1/r^2)\partial _{\varTheta }^2$. The tangential stress boundary conditions are automatically satisfied in these approximations. As detailed in Galeano-Rios et al. (Reference Galeano-Rios, Milewski and Vanden-Broeck2017), this leading-order theory is valid in the weakly viscous limit when 
$Re=\sqrt {We}/Oh \gg 1$. A similar bath model was used by Blanchette (Reference Blanchette2016, Reference Blanchette2017), although the viscous correction term was not included in the dynamic boundary condition (3.3).
We assume that the impact occurs in a bath of some viscous fluid which is subject to two boundary conditions, 
$\partial _n \phi = 0$ on the walls of the bath and 
$\partial _z \phi =0$ on the bottom of the bath, where 
$n$ is the outward facing normal of the walls of the bath (Benjamin & Ursell Reference Benjamin and Ursell1954). The former condition implies that 
$\partial _n \eta =0$ on the walls of the container. These conditions correspond physically to a bath where the working fluid maintains a constant contact angle of 
$90^{\circ }$ at the walls, and has no-flux boundary conditions along the walls and bottom. Applying these conditions to the governing system of equations, an orthogonal function expansion for the unknowns 
$\eta$, 
$\phi$ and their derivatives can be explicitly written.
The orthogonal basis functions 
$S_{j,m}(r,\varTheta )$ ultimately must satisfy
\begin{equation} (\nabla^2 + k_{j,m}^2) S_{j,m} = 0 , \end{equation}
inside of any bounding curve 
$K$ that exists on the border of the bath and the boundary condition of 
$\partial _n S_{j,m} = 0$ on 
$K$. Here 
$k_{j,m}$ are the eigenvalues of the system, and depend on the choice of the physical domain of the problem. If the boundary curve 
$K$ is a circle, the functions become 
$S_{j,m} = {\rm J}_j(k_{j,m} r) \cos {j \varTheta }$. Here 
${\rm J}_j$ are Bessel functions of the first kind and 
$k_{j,m}$ are the solutions to 
${\rm J}'_j(k_{j,m}b) = 0$, where 
$b$ is the bath radius. If we further assume the problem to be axisymmetric, then we can choose 
$j=0$, and write 
$S_{0,m} = {\rm J}_0(k_{0,m} r).$ For convenience, we define 
$k_{0,m}=k_m$ henceforth. We then can express the free surface elevation as
\begin{equation} \eta(r,t) = \sum_{m=0}^{\infty} a_m(t) {\rm J}_0(k_{m} r) . \end{equation}
We then rewrite all of the unknowns of the axisymmetric bath problem, using (3.1) and (3.2), as a function of the time varying amplitude coefficients, 
$a_m(t)$:
$$\begin{gather} \eta(r,t) = \sum_{m=0}^{\infty} a_m(t) {\rm J}_0(k_m r) , \end{gather}$$
$$\begin{gather}\kappa(r,t)=\nabla^2 \eta ={-}\sum_{m=0}^{\infty} k_m^2 a_m {\rm J}_0(k_m r) , \end{gather}$$
$$\begin{gather}\phi(r,z,t) = \sum_{m=0}^{\infty} \left( \frac{{\rm d} a_m}{{\rm d}t} + 2 \nu k_m^2 a_m \right) {\rm J}_0(k_m r) \frac{\cosh{k_m(h_0+z)}}{k_m\sinh{k_m h_0}} . \end{gather}$$In order to arrive at the final equations of motion for the free surface, we take the decompositions ((3.7)–(3.9)) and substitute them into the dynamic boundary condition (3.3). Rearranging, we find
\begin{align} \sum_{m=0}^{\infty} \left[\frac{{\rm d}^2 a_m}{{\rm d}t^2} + 4 \nu k_m^2 \frac{{\rm d}a_m}{{\rm d}t} \,{+} \left(\frac{\sigma k_m^2}{\rho} + g\right)k_m\tanh{(k_m h_0)}a_m \right ] {\rm J}_0(k_mr) ={-}\frac{p_s}{\rho}k_m\tanh{(k_m h_0)}. \end{align}Each wave mode in the bath is described by a forced, damped harmonic oscillator equation.
3.2. Droplet interface model
Additionally, we wish to recover a similar set of equations that describe the gravity-capillary waves present in the droplet, and then couple these equations to the motion of the bath. The full derivation of the droplet oscillation model can be found throughout prior works (Lamb Reference Lamb1924; Tsamopoulos & Brown Reference Tsamopoulos and Brown1983; Courty, Lagubeau & Tixier Reference Courty, Lagubeau and Tixier2006; Chevy et al. Reference Chevy, Chepelianskii, Quéré and Raphaël2012; Balla, Tripathi & Sahu Reference Balla, Tripathi and Sahu2019) and is briefly summarized below.
We begin by utilizing spherical harmonics to decompose the droplet radius in a spherical domain,
\begin{equation} \xi(\theta,t) = R + \sum_{l=1}^{\infty} \sum_{n={-}l}^l \beta_l^n(t) Y_l^n(\cos{\theta},\phi) , \end{equation}
with 
$Y_l^n = P_l^n(\cos {\theta }) e^{i n \phi }$. Due to the axisymmetry of the problem, we set 
$n=0$ and the spherical harmonics reduce to associated Legendre polynomials, 
$P_l^n (\cos (\theta ))$. For convenience, we write 
$\beta _l^n=\beta _l$ henceforth. We assume that the velocity potential takes the same form as the decomposition of the interface (Lamb Reference Lamb1924; Balla et al. Reference Balla, Tripathi and Sahu2019). We then turn to an energy conservation equation of the form
\begin{equation} \frac{{\rm d} T}{{\rm d}t} = \dot{W} - \epsilon , \end{equation}
with 
$T = K + G$, 
$\dot {W}$ and 
$\epsilon$, as the total energy (sum of the kinetic energy 
$K$ and potential energy 
$G$) of the drop, the rate of work done on the droplet interface, and the viscous dissipation, respectively. We can express 
$T$, 
$\dot {W}$ and 
$\epsilon$ using the decomposition (3.11), and substitute these expressions into the conservation of energy equation. Then, utilizing the linearized kinematic boundary condition yields a set of forced, damped harmonic oscillators that describe the amplitude of each individual spherical mode,
\begin{align} \sum_{l=1}^{\infty} \left[ \frac{{\rm d}^2 \beta_l}{{\rm d}t^2} + 2 \alpha_{l} \frac{{\rm d} \beta_l}{{\rm d}t} + \omega_{l}^2 \beta_l \right ] = \sum_{l=1}^{\infty} \left[ -\frac{(2l+1)l}{2\rho R} \int_0^{\rm \pi} p_s(\theta) \sin{\theta} P_l(\cos{\theta}) \,{\rm d}\theta + g \delta_{1l} \right ] . \end{align}We drop the sums, and arrive at the result
\begin{equation} \frac{{\rm d}^2 \beta_l}{{\rm d}t^2} + 2 \alpha_{l} \frac{{\rm d} \beta_l}{{\rm d}t} + \omega_{l}^2 \beta_l ={-}\frac{(2l+1)l}{2\rho R} \int_0^{\rm \pi} p_s(\theta) \sin{\theta} P_l(\cos{\theta}) \,{\rm d}\theta + g \delta_{1l}, \end{equation}with
$$\begin{gather} \alpha_{l} = (2l+1)(l-1)\frac{\mu}{\rho R^2} , \end{gather}$$
$$\begin{gather}\omega_{l}^2 = l(l-1)(l+2) \frac{\sigma}{\rho R^3} \end{gather}$$
and 
$\delta _{1l}$ is the Kronecker delta function. This model is valid in the weakly viscous limit, when 
$Oh\ll 1$. An extension of the free droplet model to arbitrary 
${Oh}$ can be found in other prior works (Miller & Scriven Reference Miller and Scriven1968; Moláček & Bush Reference Moláček and Bush2012; Chandrasekhar Reference Chandrasekhar2013).
3.3. Pressure forcing during impact
There is still an additional unknown in the bath mode (3.10) and drop mode (3.14) equations: the applied pressure distribution 
$p_s(r,t)$. This is generally a function of the properties of the fluid medium, the impacting speed of the object, the shape of the impacting object and the motion of the gas that surrounds the fluid. However, in this work, we will assume that the viscosity of the ambient gas is small relative to the fluid bath such that the flow within the small air film is negligible, and the pressure acts solely to apply an upwards hydrodynamic force on the droplet. In non-dimensional terms, we can construct two additional restrictions for our model, 
$\rho _g / \rho \ll 1$ and 
${Oh}_g = \mu _g / \sqrt {\rho \sigma R} \ll {Oh}$ following prior work (Moláček & Bush Reference Moláček and Bush2012). In the current experimental and DNS work, 
$\rho _g/ \rho \sim O(10^{-3})$ and 
${Oh}_g \sim O(10^{-4})$, thus the influence of these two additional parameters is indeed negligible. We have verified this for our typical experimental parameters through DNS, and find that both a four-fold increase and decrease in the ambient density and viscosity from air properties at standard temperature and pressure (STP) produces negligible changes to the trajectory, instantaneous shape of the droplet, and free surface shape throughout the interaction of the droplet and bath (Appendix A).
The radial extent of the pressure distribution is generally unknown for impact problems, and constitutes an additional problem that we must solve. In the present work, we assume that this unknown pressure distribution takes the form
\begin{equation} p_s(r,t) = \frac{F(t)}{{\rm \pi} R^2} H_r(r/r_c(t)), \end{equation}
where 
$F$ is the instantaneous magnitude of the contact force, evaluated at 
$r=0$ and 
$H_r$ is an assumed spatial profile of the pressure in the contact region. For this distribution, we can use a function that resembles the true shape of the pressure distribution during contact. The contact region, 
$A_c$, will be assumed to a simply connected disk, following Galeano-Rios et al. (Reference Galeano-Rios, Milewski and Vanden-Broeck2017) and Korobkin (Reference Korobkin1995). This allows us to write a single unknown 
$r_c(t)$ to fully describe the temporal evolution of this region of contact. Blanchette (Reference Blanchette2016, Reference Blanchette2017) used a fixed parabolic pressure shape function
\begin{equation} H_r(r) = \begin{cases} C \left( 1 - \left(\dfrac{r}{R}\right)^2 \right), & r\leq R ,\\ 0, & r > R . \end{cases} \end{equation}
Here, 
$C$ is the constant magnitude of the pressure at 
$r=0$ and 
$R$ is the undeformed radius of the droplet. Blanchette (Reference Blanchette2016, Reference Blanchette2017) chose the value of the magnitude 
$C$ such that 
$\int _0^b p_s r \,{\rm d}r= {\rm \pi}R^2$. Thus, the pressure acting on the bath interface in the respective models had a constant pressure shape function 
$H_r$ for all times during contact. However, simulation results from Galeano-Rios et al. (Reference Galeano-Rios, Milewski and Vanden-Broeck2017, Reference Galeano-Rios, Cimpeanu, Bauman, MacEwen, Milewski and Harris2021) show that the shape of the pressure distribution at the surface of a fluid bath due to an impacting, non-wetting sphere is flatter and more similar to a top-hat function for most times, and that the spatial extent of the distribution changes continuously with time during impact. Additionally, for a droplet impacting on a solid surface, the pressure in the air film has been inferred by de Ruiter et al. (Reference de Ruiter, Lagraauw, Van Den Ende and Mugele2015). The air film thickness during a bounce was measured using interferometry and the pressure estimated using a lubrication model. The film pressure in both the impacting and rebounding regimes is approximately uniform, with deviations from uniformity only near the edge of the film. In related work, the impact pressure between a droplet and a wettable solid substrate has been studied extensively by Mandre, Mani & Brenner (Reference Mandre, Mani and Brenner2009) and Mani, Mandre & Brenner (Reference Mani, Mandre and Brenner2010), and their results indicate that the impact pressure increases sharply near the contact line, likely a consequence of the decreased air film thickness in that region. For the case of droplets bouncing on a deep pool, Tang et al. (Reference Tang, Saha, Law and Sun2019) measured the air film thickness and found the film thickness to be significantly more uniform in both impacting and rebounding stages for 
${We}$ values similar to those explored in the present work, presumably as a result of the deformability of the substrate and impactor. Our predictions from DNS (presented and discussed in § 4) similarly suggest a more uniform air film thickness for the present problem, and a nearly uniform pressure profile during all stages of rebound.
We will use a simple polynomial that resembles a smoothed top hat in this work, with
\begin{equation} H_r(r/r_c(t)) = \begin{cases} C \left( 1 - \left(\dfrac{r}{r_c(t)}\right)^6 \right), & r\leq r_c(t), \\ 0, & r > r_c(t) . \end{cases} \end{equation} In order to remain consistent with our linearization, we do not allow 
$r_c$ to exceed 
$R$. Requiring that the integral of the pressure over the contact area is 
$F(t)$, we find
\begin{equation} \int_0^{b} H_r(r/r_c) r \, {\rm d}r = \frac{R^2}{2} , \end{equation}
which sets the constant 
$C$ in the pressure shape function 
$H_r(r/r_c)$. Our bath model relies on the decomposition of the fluid motion into a linear superposition of infinitely many waves with wavenumbers 
$k_m$. Therefore, we apply a similar decomposition to this pressure function 
$p_s=\sum _{m=0}^{\infty } \,{\rm d}_m {\rm J}_0(k_m r)$. Since we are working in a cylindrical domain, we will choose the zeroth-order Bessel function of the first kind as our orthogonal basis function, and thus the 
$d_m$ are the Fourier–Bessel coefficients of the function 
$H_r$:
\begin{equation} d_m = \frac{2}{(b {\rm J}_1(k_m))^2}\int_0^b H_r(r) r {\rm J}_0(k_m r) \,{\rm d}r, \end{equation}
with the domain extending from 
$r=[0,b]$. The reconstruction of the top-hat function in Fourier–Bessel space converges too slowly to be of practical use (Storey Reference Storey1968), also noted by Blanchette (Reference Blanchette2016), and as such we use a polynomial expression that resembles a smoothed top hat. Additionally, we tested higher-order polynomials (corresponding to a larger flat region), and found increasingly poor convergence behaviour, similar to that of the top hat (see Appendix B for a case study on the sensitivity of the results to the choice of shape function). The ultimate choice of shape function used here thus represents a practical compromise.
Substituting in the definition of the pressure (3.17) into (3.10), performing the Fourier–Bessel decomposition, we find
\begin{align} \frac{{\rm d}^2 a_m}{{\rm d}t^2} + 4 \nu k_m^2 \frac{{\rm d}a_m}{{\rm d}t} + \left(\frac{\sigma k_m^2}{\rho} + g\right)k_m\tanh{(k_m h_0)}a_m ={-}\frac{F}{\rho {\rm \pi}R^2} d_{m} k_m\tanh{(k_m h_0)} , \end{align}
which govern the evolution of bath wave modes 
$m$. Similarly, substituting the pressure (3.17) into (3.14), we write
\begin{equation} \frac{{\rm d}^2 \beta_l}{{\rm d}t^2} + 2 \alpha_{l,0} \frac{{\rm d} \beta_l}{{\rm d}t} + \omega_{l,0}^2 \beta_l ={-}\frac{F}{2 {\rm \pi}\rho R^3} c_{l} (2l+1) l + g \delta_{1l} . \end{equation}
The coefficients 
$c_{l}$ result from the mode decomposition of the projection of the pressure into spherical space,
\begin{equation} c_{l} = \int_0^{\rm \pi} H_r(\theta) \sin(\theta) P_l(\cos{\theta}) \,{\rm d}\theta , \end{equation}
which naturally arise in the derivation of (3.14). Additionally, the definition of the pressure (3.17) reduces the governing equation of the 
$l=1$ ‘translational’ mode to
\begin{equation} \frac{{\rm d}^2 \beta_1}{{\rm d}t^2} ={-}\frac{F}{m} + g , \end{equation}
which clearly governs the droplet centre of mass motion 
$\beta _1=-z_{cm}$. Evidence from the simulations of Galeano-Rios et al. (Reference Galeano-Rios, Milewski and Vanden-Broeck2017) indicates that impact trajectory is very sensitive to the instantaneous size of the contact area. Utilizing a constant pressed area for the pressure, as done in Blanchette (Reference Blanchette2016), does not produce results that compare well with experiment (Appendix B), particularly for cases of small 
${We}$. Appendix B details how the choice of this pressure shape function modifies the predicted trajectory of the droplet for typical experimental parameters. The trajectory is largely insensitive to the choice of pressure shape function, but incorporating a time-dependent contact radius is essential for agreement. The method for determining both 
$F(t)$ and 
$r_c(t)$ are discussed in the next section.
3.4. Modelling contact
The contact force 
$F(t)$ is determined through the use of a ‘1-Point’ kinematic match (1PKM) condition. Essentially, we enforce contact only at a single point; the centre of our axisymmetric domain. Thus, the additional constraint can be written as
\begin{equation} \eta(r=0,t) = z_{cm}(t) - \xi(\theta=0,t) = \sum_{m=0}^{\infty} a_m(t) = z_{cm}(t) -\left(R + \sum_{l=2}^{\infty}\beta_l(t)\right) . \end{equation}
This additional constraint allows us to determine the unknown contact force 
$F(t)$. Contact between the droplet and the bath ends when the magnitude of the contact force as predicted by the kinematic match becomes negative. We note that this 1PKM model is a significant simplification of the full kinematic match successfully used to study related impact problems (Galeano-Rios et al. Reference Galeano-Rios, Milewski and Vanden-Broeck2017, Reference Galeano-Rios, Milewski and Vanden-Broeck2019, Reference Galeano-Rios, Cimpeanu, Bauman, MacEwen, Milewski and Harris2021). The full kinematic match predicts the evolution of the contact area and the contact pressure distribution (without requiring an assumption for 
$H_r$) by imposing natural geometric and kinematic constraints, essentially considering additional equations to solve at each time step. The algorithm requires iteration at each time step, and the minimization of a tangency boundary condition is used to determine the correct contact area and pressure shape.
Lastly we turn to the unknown contact radius 
$r_c(t)$. By not restricting the deformations of the bath and droplet interface with the use of additional tangency and distributed kinematic match conditions, we find that the results of our simulation consistently produce interfacial shapes that cross each other. The amount of overlap between the two interfaces is generally small, for the typical experimental parameters in our study the maximum overlap is less than 
$0.05R$. However, we can use the predictions from both interfacial models to determine the exact location where the two interfaces cross and separate, and use this as the instantaneous radius of contact 
$r_c(t)$. Thus, at each time, contact between the bath and drop is ensured at both 
$r=0$ and 
$r=r_c(t)$. In order to enforce contact within the entirety of the contact region, a full kinematic match would be required – this circumvents the need for any assumptions on the pressure profile shape, but is substantially more computationally expensive. Our contact radius criterion is similar to that of the numerical model presented in prior work on droplet rebound from solid substrates (Moláček & Bush Reference Moláček and Bush2012). While this method is unphysical, it yields accurate predictions for the contact radius as compared with DNS as demonstrated in § 5.
3.5. Summary
Choosing a time scale of 
$t_{\sigma }=\sqrt {\rho R^3/\sigma }$ (with 
$\tau = t/t_{\sigma }$), a length scale of 
$R$, and a force scale of 
$2 {\rm \pi}\sigma R$ (with 
$f = F / 2 {\rm \pi}\sigma R )$, we recast the governing equations in non-dimensional form as
$$\begin{gather} \eta(r,\tau) = \sum_{m=0}^{M} a_m(\tau) {\rm J}_0(k_m r) , \end{gather}$$
$$\begin{gather}\xi(\theta,\tau) = 1 + \sum_{l=2}^{L}\beta_l(\tau) P_l(\cos{\theta}) , \end{gather}$$
$$\begin{gather}\frac{{\rm d}^2 {a_m}}{{\rm d}\tau^2} + 4 {Oh} {k_m}^2 \frac{{\rm d} {a_m}}{{\rm d}\tau} + ({k_m}^2 + {Bo}){k_m}\tanh{({k_m} {h_0})}{a_m} ={-}2{f} d_{m,0} {k_m} \tanh{({k_m} {h_0})} , \end{gather}$$
$$\begin{gather}\frac{{\rm d}^2 {\beta_l}}{{\rm d}\tau^2} + 2 {Oh} (2l+1)(l-1) \frac{{\rm d} {\beta_l}}{{\rm d}\tau} + l(l-1)(l+2) {\beta_l} ={-} {f} c_{l} (2l+1)l + {Bo}\delta_{1l}, \end{gather}$$
$$\begin{gather}\frac{{\rm d}^2 z_{cm}}{{\rm d} \tau^2} = \frac{3}{2} f - {Bo} , \end{gather}$$
$$\begin{gather}\eta(0,\tau) = {z}_{cm}(\tau) - {\xi}(\theta=0,\tau) . \end{gather}$$
Equations (3.29) and (3.30) describe the evolution of the bath and droplet oscillation modes, respectively. Equation (3.31) governs the vertical motion of the droplet's centre of mass. Equation (3.32) couples these equations all together, with 
$r=0$ as the single point of ‘contact’ enforced between the droplet and the bath and allows for determination of the unknown 
$f(\tau )$. These equations are solved using standard ODE numerical integration techniques. The shape of the bath and droplet can be reconstructed at any time 
$t$ via the sums in (3.27) and (3.28), respectively.
The complete model is valid when 
$Re = \sqrt {We}/Oh \gg 1$ and 
$Oh\ll 1$. Also, since the model is linearized about the undeformed state, we anticipate it to hold when deformations remain small, further suggesting 
$Bo\ll 1$ and 
$We\ll 1$. However, we later demonstrate through direct comparison with experiment and DNS that the model remains predictive even for moderate 
${We}$.
3.6. Numerical methods
We solve these equations using a backward Euler method, ensuring a minimum of 100 time steps within the inertio-capillary time 
$t_{\sigma }$. An implicit method was chosen, following Galeano-Rios et al. (Reference Galeano-Rios, Cimpeanu, Bauman, MacEwen, Milewski and Harris2021), as the instantaneous size of the pressure distribution acting on both the droplet and bath at the next time step is unknown. Treating the pressure explicitly on either the droplet or the bath can lead to non-physical behaviour in the system. We used 
$M=150$ modes for the bath interface and 
$L=55$ modes for the droplet interface. These values were determined by running simulations of a 
${We} =0.7$, 
${Bo} =0.017$, 
${Oh} =0.006$ impact and assessing convergence as described in what follows. First, we kept the number of droplet modes fixed at 
$L=15$ and increased the number of bath modes from 
$30$ to 
$500$ in increments of 
$25$. Then, the simulation was run again, fixing the number of bath modes at 
$75$ and increasing droplet modes from 
$15$ to 
$200$. Sufficient convergence was determined if the maximum absolute value of the difference in centre of mass trajectories during contact changed by less than 
$1\,\%$. Finally, both the droplet and bath number of modes were increased simultaneously, and convergence was still observed. These values are similar to comparable to those found in Blanchette (Reference Blanchette2016) (
$M=200$, using sine functions as the basis functions in a square bath), and Moláček & Bush (Reference Moláček and Bush2012) (
$L=150$, but found good accuracy in comparison with experimental data at 
$L=20$). Additionally, once mode convergence was determined, we decreased the time step of the simulation in increments until time step convergence was similarly reached using the same criterion. Unexpectedly, in a fully converged simulation, there is a time step threshold below which the algorithm results in unstable oscillations in the magnitude of 
$F(t)$. This time step threshold is typically at least three orders of magnitude smaller than 
$t_{\sigma }$. This apparent instability deserves awareness and future attention, but does not affect the results presented in the present work. The bath size was set to 
$b=25 R$ which was determined to be sufficiently large such that reflected waves did not influence the droplet during impact. In order to find the instantaneous contact radius, we take two line segments from the reconstruction of the droplet and bath interface shapes. From these we can write a linear system of equations for four unknowns: the 
$(r,z)$-pair of the intersection location, the normalized distance from the starting point of the first line segment to the intersection, and the normalized distance from the starting point of the second line segment to the intersection (Schwarz Reference Schwarz2022). We then loop over every line segment to find every intersection. We take the largest of this set as 
$r_c(t)$. In the reconstructions of the interfaces at each time step, at least 5000 points are used in both 
$\theta$ and 
$r$ to ensure that error is minimized. All code associated with the implementation of the model is available at https://github.com/harrislab-brown/BouncingDroplets.
4. Direct numerical simulation
Complementing the previous investigative tools in terms of both experimental and modelling capabilities, we build a dedicated computational framework within the open-source solver Basilisk. This has the dual aim of validation and exploration of quantities of interest outside the reach of previous methodologies. Basilisk (Popinet Reference Popinet2015) and its predecessor Gerris (Popinet Reference Popinet2003, Reference Popinet2009) have been widely adopted by the computational fluid dynamics community over the past two decades. Its second-order accuracy in both space and time, alongside adaptive mesh refinement and parallelization features, have led to successful studies of multiscale fluid systems such as the scenario here. The equations for conservation of momentum and mass are solved as part of a one-fluid formulation within the volume-of-fluid (VOF) methodology, with the Bell–Collela–Glaz advection (Bell, Colella & Glaz Reference Bell, Colella and Glaz1989) employed for advective terms using a Courant–Friedrichs–Lewy (CFL)-limited time-stepping strategy, a well-balanced surface tension implementation (Popinet Reference Popinet2018) and an implicit treatment of the viscous terms.
The non-dimensionalization described in table 1 is retained, with the drop diameter 
$R$ and initial impact velocity 
$V_0$ providing the reference length scale and velocity scale in the system. The additional gas region, fully captured in the DNS, has physical properties modelled from typical air values at room temperature. In particular we highlight the value of density ratio 
$\rho _g/\rho \approx 0.0012$ and 
${Oh}_g = O(10^{-4}) \ll {Oh}$, as outlined in table 1, both consistent with the modelling assumptions described in § 3.3. The axisymmetric computational domain is constructed as a 
$20R \times 20R$ square, with 
$20R$ proving sufficient to avoid any artifacts from boundary reflections, while describing the dynamics of the impact (Galeano-Rios et al. Reference Galeano-Rios, Cimpeanu, Bauman, MacEwen, Milewski and Harris2021). The pool height is set to 
$z=10R$, rendering bottom effects negligible, while also allowing a generous region occupied by air for the droplet bounce to be quantified. Given that scales range from microns (in the gas layer between drop and pool) to centimetres (for the full domain size), this is a challenging set-up which takes full advantage of the numerical capabilities available. The adaptive mesh refinement strategy prioritizes fluid–fluid interfacial location and changes in the magnitude of the velocity components in order to concentrate resources where needed. Extensive validation efforts have determined that a minimum grid cell size of approximately 
$3\,\mathrm {\mu }{\rm m}$ (or just above 
$100$ cells per radius) is sufficient to ensure mesh independence for the observed metrics, with additional refinement having been further tested in sensitive regimes and yielding no substantial benefit. This leads to a grid cell count of 
$20\,000\unicode{x2013}80\,000$ over the duration of a simulation, with the workload typically distributed across 4–8 cores over approximately 
$48$ CPU core hours per run. A subset of the computational domain, including the grid cell distribution, is shown in figure 3(a). We note that, with a uniform mesh, a similarly resolved computation would require in excess of 
$4.2 \times 10^6$ grid cells, rendering full parameter studies intractable. Similar deployment of resources has previously proven successful in the investigation of multifluid systems in the context of impact regimes ranging from small (Galeano-Rios et al. Reference Galeano-Rios, Cimpeanu, Bauman, MacEwen, Milewski and Harris2021) to moderate (Fudge et al. Reference Fudge, Cimpeanu and Castrejón-Pita2021) and finally large (Cimpeanu & Moore Reference Cimpeanu and Moore2018) velocities, resulting in tools and improved insight into physical phenomena such as bouncing, coalescence and splashing. The implementation described above is made available to interested users at https://github.com/rcsc-group/BouncingDroplets.
Figure 3. The DNS detail of the rebound of a 
$R=0.35$ mm deionized water droplet from a bath of the same fluid corresponding to 
${We} =0.7$, 
${Bo} =0.017$ and 
${Oh}= 0.006$. (a) Adaptive mesh refinement strategy overlaid to highlighted drop, pool and gas regions (i); 
$2$-norm of the velocity field 
$\|\boldsymbol {u}\|_2$ and fluid–fluid interfaces (ii) at dimensionless 
$t=1.5$. (b) Thickness of the air gap entrapped between the impacting drop and the pool, with three highlighted time steps during descent, maximal spread and rebound, presented in dimensional form. (c) Raw pressure data extracted at the centre of the air gap, along the contact radius, for the highlighted time steps in panel (b). Movie supplementary material available at https://doi.org/10.1017/jfm.2023.88 containing additional detail is also provided.
A particularity of our set-up lies in the integration of the recent functionality developed for non-coalescence scenarios, as previously employed by Ramírez-Soto et al. (Reference Ramírez-Soto, Sanjay, Lohse, Pham and Vollmer2020) and recently by Sanjay et al. (Reference Sanjay, Sen, Kant and Lohse2022b). With different symbolic definitions for underlying colour functions representing the droplet and the bath in the VOF framework, the algorithm ensures that numerically induced coalescence is avoided, and the entrapped gas region between the impacting drop and the pool is well resolved throughout the studied motion. A sufficiently high-resolution level is still required for 
$O(1)\, \mathrm {\mu }\text {m}$ thicknesses to be maintained (Tang et al. Reference Tang, Saha, Law and Sun2019), as illustrated in figure 3(b), with the non-coalescence package allowing suitable subgrid cell level length scales to be reached. The region highlighted in colour represents the approximate contact area as a function of space and time, with the contact radius limit defined by a gap of twice the typical width of the gas film trapped between the drop and pool being reached as one navigates radially outwards (a definition consistent with the mathematical model). The typical gas film profile in space is given by a nearly uniform value in the bulk of the drop–pool contact area, with a thinning observed towards the very end, and the smallest scales reached therein approaching 
$0.25\, \mathrm {\mu }{\rm m}$. While a direct comparison with existing experimental data would require a more dedicated set-up and is beyond the scope of this study, similar length scales were recently observed by Tang et al. (Reference Tang, Saha, Law and Sun2019). Finally, this study also constituted one of our most stringent convergence tests, with more refined mesh levels producing only negligible changes in this observed space–time map. The time variation of the contact radius is a clear hallmark of this behaviour, with an initially rapid increase to almost 
$85\,\%$ of the initial droplet radius being followed by a slower decrease as restorative forces push the droplet back up during the ascent stage, consolidating this as a key assumption to be retained as part of the mathematical model presented previously. While the above discussion focuses on a representative set of parameters, the generated outcomes and insight apply to the much wider parameter set we explore.
Figure 3 (accompanied by movie supplementary material) also outlines velocity field information inside each of the fluid phases (figure 3aii), as well as the aforementioned gas film thickness (figure 3b) and pressure distribution (figure 3c) across relevant flow stages. Manipulating raw pressure data in projection methods is a known challenge (Philippi, Lagrée & Antkowiak Reference Philippi, Lagrée and Antkowiak2016; Negus et al. Reference Negus, Moore, Oliver and Cimpeanu2021), with the oscillatory behaviour observed in figure 3(c) linked to the VOF approximation on a structured grid, and the underlying colour function sampled across various approximation points inside the cells in near vicinity of interfaces with large value contrasts in the quantities of interest. We underline in particular that the aforementioned variations oscillate around well-defined means within the contact region. The pressure in the air gap across the evolving contact radius may thus be robustly approximated by a top-hat function at every flow stage investigated, reinforcing previous modelling choices in § 3.3 and Appendix B.
After having undergone stringent numerical verification procedures, a typical study in our parameter regime of interest is showcased in figure 4 and described in detail in the following section. The agreement between experiment, model and simulation is very encouraging, with the developed resources in an excellent position to bridge one another and comprehensively explain the target rebound metrics.
Figure 4. Rebound of a 
$R=0.35$ mm deionized water droplet from a bath of the same fluid corresponding to 
${We}=0.7$, 
${Bo}=0.017$ and 
${Oh} = 0.006$. (a) Trajectory comparison between the experiment (red dashed line with typical variation shown as dotted red lines), DNS (black dashed line) and quasipotential model (blue solid line). (b) Instantaneous contact radius, normalized by the undeformed radius 
$R$, as a function of time for the quasipotential model and DNS. (c) Maximum width of the droplet 
$w$, as a function of time for the quasipotential model and DNS. (d) Comparison of droplet shape between experiment, DNS and quasipotential model. Supplementary movies are available.
5. Results
In this section, we first present the results of a direct comparison between experiment, quasipotential model and DNS for a single impact 
${We}$. Then, we vary 
${We}$ for two working fluids, and compare the results of the three different impact metrics between the experiment, DNS and model. Having validated the model and DNS, we then run sweeps over 
${Bo}$ and 
${Oh}$ to deduce the effect that these non-dimensional constants have on droplet rebound metrics, and compare predictions with existing experimental data sets available in the literature.
5.1. Comparison with experiment
We first consider a single impact of a deionized water droplet with 
$R=0.35$ mm. A direct comparison of the trajectory results of the model, DNS and the experiments is depicted in figure 4(a). We find excellent agreement for the top, centre of mass and bottom trajectories between the experiment, quasipotential model and DNS. Additionally, we make a direct comparison between the quasipotential model and the DNS for the prediction of the radius of contact in figure 4(b). During impact, the quasipotential model accurately predicts the instantaneous contact area, as well as the maximum contact area. The DNS, which accurately resolves the air film, shows that a finite region of contact is already developed prior to impact. At 
$t/t_{\sigma } > 3$, the two models deviate from each other, and this is most likely due to suction effects as the droplet pulls away from the surface. These effects cannot be predicted while neglecting the motion of the air, as in the current quasipotential model, however, do not appear consequential to the overall droplet trajectory. A comparison between the two models for the maximum width of the droplet is depicted in figure 4(c). As is the case with the contact radius, the deformation of the droplet in the DNS occurs slightly before that of the quasipotential model, as the pressure in the film is already building prior to contact. This leads to an overall phase shift of the oscillation between the models, although the maximum value for the deformation between both models remains in very good agreement. We also compare the droplet and interface shape between the experiments, DNS, and linearized model can be seen in the panels of figure 4(d). The impact is depicted for six different instants during contact. As the droplet deforms the interface, a capillary wave travels from the impact location to the north pole of the droplet. The collapse of this wave onto itself occurs just before the time of maximum deformation of the bath, and corresponds to the maximum deformation of the droplet. After this time, surface tension in the bath begins to act to restore the equilibrium, having redistributed some of the initial impact energy in the form of interfacial waves. The droplet remains mostly spherical as the bath relaxes, until contact is lost. During free flight, the droplet oscillates as an underdamped harmonic oscillator, dissipating additional energy through viscosity. Both the DNS and experiment show slightly stronger oscillations in the top of the droplet as compared with the model, to the point which the instantaneous slope at the top is occasionally close to zero. Overall, the DNS and quasipotential model are in excellent agreement and predict the bath shape, droplet shape and droplet trajectory with high accuracy for these parameters.
As we begin to explore the larger parameter space, we consider three different output parameters for the rebounds: coefficient of restitution (
$\alpha$); contact time (
$t_c$); and maximum surface deflection (
$\delta$). As mentioned in § 2.2, given the experimental difficulty of accurately determining the time of surface detachment of the droplet, contact time, 
$t_c$, is defined in the experiment as the interval between the two instances when the north pole of the droplet crosses level 
$z=2R$ and the coefficient of restitution, 
$\alpha$, is defined as minus the ratio of the vertical velocities at those times. For the model and DNS, we define the metrics in the same way, but when the centre of mass of the droplet crosses 
$z=R$. This is chosen because a measure on the centre of mass more accurately describes the total translational energy transfer from the droplet. However, in comparing the results from the model and DNS using both top (measured at 
$z=2R$) and centre of mass (measured at 
$z=R$), we found a typical difference of 
$2\,\%$ for 
$\alpha$ and 
$t_c/t_{\sigma }$ in the silicone oil experiments, and 
$5\,\%$ for the same parameters in the deionized water experiments.
Figure 5 outlines a variation of impact 
${We}$ for a deionized water droplet onto a water bath. In this parameter sweep, the coefficient of restitution generally decreases as the 
${We}$ number is increased, eventually saturating at a value just below 
$0.3$, and remaining nearly independent of the 
${We}$ number thereafter. The contact time also decreases as the 
${We}$ number is increased, but remains relatively independent of the 
${We}$ number at an earlier value, consistent with results found for impact on a solid surface (Richard et al. Reference Richard, Clanet and Quéré2002) and for previous experiments for impact on a deep pool (Zhao et al. Reference Zhao, Brunsvold and Munkejord2011). The maximum penetration depth increases monotonically with 
${We}$. We find good agreement between the model, DNS and experiments with regard to the restitution coefficient, contact time and maximum surface deflection for these experiments. Additionally, the quasipotential model is able to predict 
$\alpha$ for the entire range of experiments. Both the model and DNS slightly underpredict the dimensionless contact time at intermediate 
${We}$, yet do agree with the experimental data for 
${We}$ 
$\leq 1$. The quasipotential model also underpredicts the penetration depth and contact time at moderate 
${We}$. Nonlinear effects associated with larger deformation have been neglected in this model, and are most likely the cause of this discrepancy.
Figure 5. (a) Coefficient of restitution, (b) contact time and (c) maximum penetration depth for a 
$R=0.35$ mm deionized water droplet rebounding from a bath of the same fluid as a function of 
${We}$ (with 
${Bo} = 0.017$ and 
${Oh} = 0.006$). Error bars on experimental data points are quantified as the standard deviation of at least five experimental trials. Predictions of the quasipotential model are shown as blue solid lines, DNS as black dashed lines.
Figure 6 depicts a 
${We}$ number variation using 
$5$ cSt silicone oil. In non-dimensional terms, this represents an increase in both the 
${Bo}$ and 
${Oh}$ numbers. Similar trends in 
$t_c/t_{\sigma }$ and 
$\delta$ are observed for the silicone oil as compared with the deionized water. However, 
$\alpha$ tends to generally increase with 
${We}$ in this case, as opposed to water which showed the opposite trend. The quasipotential model accurately predicts 
$\alpha$ for almost the full range of 
${We}$, with slight underprediction for 
${We}$ 
$<2$. Similar to the water case, the quasipotential model underpredicts 
$t_c/t_{\sigma }$ and 
$\delta$ at intermediate 
${We}$, with the DNS capturing these metrics more accurately. Although not verified experimentally, the model and DNS predict that droplets with very, very low 
${We}$ numbers cease to return to their original height at all. Note that 
${We} \leq 0.5$ is challenging to explore experimentally for these parameters, as pinch off of the droplets from the generator induces oscillations that need to dampen out prior to impact. The short free flight time and low viscosity of these extremely low 
${We}$ cases mean that there is still oscillation present at impact, which has been shown in prior work to influence rebound dynamics in related droplet impact problems (Biance et al. Reference Biance, Chevy, Clanet, Lagubeau and Quéré2006; Yun Reference Yun2018).
Figure 6. (a) Coefficient of restitution, (b) contact time, and (c) maximum penetration depth for a 
$R=0.35$ mm 5 cSt oil droplet rebounding from a bath of the same fluid as a function of 
${We}$ (with 
${Bo} = 0.056$ and 
${Oh} = 0.058$). Error bars on experimental data points are quantified as the standard deviation of at least five experimental trials. Predictions of the quasipotential model are shown as blue solid lines, DNS as black dashed lines.
5.2. Inertio-capillary limit
The validated quasipotential model can now be used to explore other sets of parameters. In particular, based on the assumptions of the model, we expect the model to remain accurate (and even perhaps improve) for cases of even smaller 
${Bo}$ and 
${Oh}$ than achieved in experiments. As a grounding point, we first turn our attention to the pure inertio-capillary limit, where both gravitational and viscous effects are ignored (i.e. 
${Bo}=0$ and 
${Oh}=0$). This case reduces the number of dimensionless parameters that describe the physical problem to just one: the 
${We}$ number. The results in the inertio-capillary limit are presented in figure 7, along with the droplet and bath shape predictions for various 
${We}$. The penetration depth 
$\delta ^*$ increases monotonically with increasing 
${We}$, as found for both the water and oil experiments. Additionally, the contact time 
$t_c^*$ decreases before becoming mostly independent of 
${We}$. However, in this limiting case, the coefficient of restitution 
$\alpha ^*$ monotonically decreases with 
${We}$ and does not have a local maximum in the restitution coefficient, and droplets can even retain a majority of their impacting energy at sufficiently low 
${We}$. Furthermore, the coefficient of restitution then remains nearly independent of 
${We}$ above 
${We}$ 
$>1.75$, and is predicted to saturate to a value of approximately 
$\alpha _{s}^* = 0.31$.
Figure 7. (a–c) Rebound parameters of a droplet in the inertio-capillary regime (as denoted by the asterisk). (d) Droplet and bath shape in the pure inertio-capillary regime (
${Bo} = 0$, 
${Oh} = 0$) for three different 
${We}$.
5.3. Influence of viscosity and gravity
We then individually probe the parameter space by increasing either 
${Bo}$ or 
${Oh}$. These variations are presented in figures 8 and 9. The result of increasing 
${Bo}$, while keeping 
${Oh}$ constant (yet negligibly small), is depicted in figure 8. At low values 
${Bo}$, the curves converge to the inertio-capillary limit as presented in the prior section. For a given 
${We}$, the coefficient of restitution decreases monotonically with 
${Bo}$, until eventually ceasing to return to its original height. At intermediate values of 
${Bo}$ the curves exhibits an interesting non-monotonic dependence on 
${We}$, with a local maxima at finite 
${We}$. Furthermore, the contact time of the droplets is predicted to increase with increasing 
${Bo}$. These qualitative trends are consistently reproduced by the DNS, with satisfactory quantitative agreement between the quasipotential model and DNS for 
${Bo} \lesssim 0.1$. For larger 
${Bo}$, and for the highest 
${We}$ cases explored, the model generally underpredicts the coefficient of restitution as compared with the DNS.
Figure 8. (a) Coefficient of restitution and (b) contact time as a function of the Bond number 
${Bo}$. Viscous effects are set to be finite but negligible, with 
${Oh}$ 
$=6\times 10^{-4}$. Predictions of the quasipotential model are shown as solid lines, DNS as individual markers. The vertical dashed lines in panel (b) reference the critical 
${We}$ for each 
${Bo}$ below which droplets do not bounce.
Figure 9. (a) Coefficient of restitution and (b) contact time as a function of the Ohnesorge number 
${Oh}$. Gravitational effects are set to be finite but negligible, with 
${Bo}= 1 \times 10^{-3}$. Predictions of the quasipotential model are shown as solid lines, DNS as individual markers.
Increasing 
${Oh}$ while keeping 
${Bo}$ negligibly small is shown in figure 9, and also predicts a monotonic decrease in the restitution coefficient, with curves converging to the inertio-capillary limit for small 
${Oh}$. Unlike the 
${Bo}$ variation, the shape of the curve remains relatively unchanged. The non-dimensionalized contact time changes only marginally, even over an order of magnitude increase in 
${Oh}$. Very similar trends are predicted by DNS, with models diverging quantitatively beyond 
${Oh} \gtrsim 0.1$, consistent with a breakdown of the weakly viscous modelling assumptions. For larger 
${Oh}$, the quasipotential model overpredicts the coefficient of restitution and underpredicts the contact time, as compared with the DNS.
5.4. Scaling analysis
In this subsection we present scaling arguments to rationalize the dependence of the coefficient of restitution on 
${Bo}$ and 
${Oh}$ detailed in figures 8(a) and 9(a), respectively. As revealed in the prior section, at a fixed 
${We}$, the coefficient of restitution decreases monotonically from the inertio-capillary limit (figure 7a) as either 
${Bo}$ or 
${Oh}$ is increased. Due to the number of parameters involved, in order to proceed, we will assume that this additional energy transfer (or loss) due to weak gravitational (or viscous) effects is independent of the baseline energy transferred (
$\Delta E^*$) in the inertio-capillary limit. Mathematically, this assumption can be expressed in the form
\begin{equation} \alpha^2 = \frac{E_o - \Delta E^* - \Delta E_{g,\mu}}{E_o} = \alpha^{*2}-\frac{\Delta E_{g,\mu}}{E_o}, \end{equation}or
\begin{equation} \alpha^{*2}-\alpha^2 = \frac{\Delta E_{g,\mu}}{E_o},\end{equation}
where 
$E_o=({2{\rm \pi} }/{3})\rho R^3 U^2$ is the initial droplet kinetic energy, and 
$\Delta E_{g,\mu }$ is the supplemental energy transferred or lost due to gravity or viscosity, respectively. In what follows, we propose scalings for these energies.
Since gravity is a conservative force, increases in 
${Bo}$ must lead to an overall increase in the deformation (and subsequent oscillation) of the bath and droplet. In the capillary-dominated regime, the gravity-induced deformation can be estimated to scale with 
$\rho R^3 g/\sigma$, corresponding to an additional surface energy scaling as 
$\Delta E_g \sim \rho ^2 R^6 g^2 / \sigma$. Normalizing this estimate by the incident kinetic energy, we find an additional fractional energy transfer to droplet–bath deformations that scales as
\begin{equation} \frac{\Delta E_g}{E_o} \sim \frac{\rho R^3 g^2}{\sigma U^2}=\frac{{Bo}^2}{We}. \end{equation}
Motivated by (5.2) and (5.3), we replot the data from figure 8(a) in 10(a), and find a satisfactory collapse. In particular, the scaling in (5.3) correctly predicts that bounces at lower 
${We}$ will be more heavily penalized in terms of their coefficient of restitution. This inequity rationalizes the observed non-monotonic behaviour of 
$\alpha$ with 
${We}$ predicted for intermediate 
${Bo}$ (figure 8a). Furthermore, (5.2) and (5.3) suggest a scaling for the ‘critical’ Weber number 
${We}_c$ below which the droplet ceases to bounce (i.e. 
$\alpha ^2<0$):
\begin{equation} We_c \sim \frac{Bo^2}{\alpha^{*2}}. \end{equation}
The data for the critical Weber number as predicted by the quasipotential model is shown in figure 10(b), and follows the proposed scaling in (5.4). We note that Blanchette (Reference Blanchette2016) observed a parabolic scaling for the critical Weber number with 
${Bo}$ in prior work, also finding this threshold to be largely independent of 
${Oh}$.
Figure 10. (a) Coefficient of restitution predictions from the quasipotential model for 
${Oh} =6\times 10^{-4}$ (figure 8a) replotted as informed by the scaling in (5.3). (b) Critical 
${We}$ as a function of 
${Bo}$; the points are predictions from the quasipotential model, and the dotted curve is a parabolic fit (
${We}_c = 11.43 {Bo}^2/\alpha ^{*2}$) consistent with the scaling in (5.4). (c) Coefficient of restitution predictions from the quasipotential model for 
${Bo} = 1 \times 10^{-3}$ (figure 9a) replotted as informed by the scaling in (5.5).
Upon the inclusion of viscosity, there is viscous dissipation in the drop and bath that now occurs during contact. The rate of viscous energy dissipation (per unit volume) can be estimated to scale as 
$\mu (\boldsymbol {\nabla } \boldsymbol {u})^2\sim \mu U^2 / R^2$. Assuming a characteristic fluid volume 
$R^3$, we find a viscous energy dissipation rate that scales like 
$\mu U^2 R$. As demonstrated in the prior sections, the contact time 
$t_c \sim t_\sigma$, and thus we may estimate the additional fractional energy loss during contact as
\begin{equation} \frac{\Delta E_\mu}{E_o} \sim \frac{\mu}{\sqrt{\sigma \rho R}}=Oh. \end{equation}
Replotting the data from figure 9(a) in 10(c) shows that apart from the very lowest 
${We}$ cases considered, the curves collapse to a single line, confirming the proposed scaling. For all 
${We}$, the curves are approximately linear in 
${Oh}$ (consistent with (5.5)) with the slope evidently depending on 
${We}$ for the smallest 
${We}$ cases.
In summary, our models predict an additional energy transfer (or loss) over the pure inertio-capillary limit when gravitational (or viscous) effects are introduced. Our scaling suggests that gravity leads to additional deformations in the system, coming at an additional energetic cost. Additionally, viscosity provides a mechanism for energy dissipation, occurring over the finite contact time of the droplet. Computing the various energies directly in DNSs (such as in Sanjay et al. (Reference Sanjay, Lakshman, Chantelot, Snoeijer and Lohse2022a)) may provide additional insight to the remaining subtleties present in the data.
5.5. Comparison with prior literature data
In the works of Jayaratne & Mason (Reference Jayaratne and Mason1964), Bach et al. (Reference Bach, Koch and Gopinath2004), Zhao et al. (Reference Zhao, Brunsvold and Munkejord2011), Zou et al. (Reference Zou, Wang, Zhang, Fu and Ruan2011) and Moláček & Bush (Reference Moláček and Bush2013) there is a reported saturation in the energy transfer from the drop during rebound at modest 
${We} > 1$ and low 
${Oh}$, as measured by the coefficient of restitution. The exact value of the saturation restitution coefficient does seem to vary, however, from 
$0.2$ in Moláček & Bush (Reference Moláček and Bush2012) for more viscous drops, to 
$0.3$ in Bach et al. (Reference Bach, Koch and Gopinath2004) and Zhao et al. (Reference Zhao, Brunsvold and Munkejord2011), to as low as 
$0.1$ (Zou et al. Reference Zou, Wang, Zhang, Fu and Ruan2011) for large 
${Bo}$ impact scenarios. Remarkably, recent experiments on rebound of liquid metal droplets in viscous media also showed similar values of the coefficient of restitution with 
$\alpha =0.27$ (McGuan, Candler & Kavehpour Reference McGuan, Candler and Kavehpour2022). A similar saturation is also observed in our experimental results presented herein, with water droplets generally bouncing higher that the 5 cSt silicone oil droplets. In figure 11, we overlay existing available experimental data for 
$\alpha$ from numerous sources and find that the prediction from the quasipotential model accurately captures much of this data. The grey line represents the extrapolation of data from non-normal impacts by (Jayaratne & Mason Reference Jayaratne and Mason1964) over the range of 
${We}$ reported in their work. Data from the experiments completed in this work utilizing 
$5$ cSt silicone oil and deionized water are included with error bars. The historical data generally indicates a decrease in 
$\alpha$ with 
${Bo}$, as captured by the present model. Despite the relatively large variation in 
${Oh}$, the experimental data appear to match the results of the quasipotential model quite well.
Figure 11. Comparison of model 
${Bo}$ predictions (figure 8a) with the existing literature data where 
${We}$ was reported. The grey markers represent data points where the exact value of 
${Bo}$ is unknown, with only ranges reported. The grey dot–dashed line represents the extrapolation of oblique impact data to normal impacts by Jayaratne & Mason (Reference Jayaratne and Mason1964). Data from the experiments completed in the present work are included with error bars. Studies included here are summarized in table 2.
Table 2. Relevant publications and ranges of droplet parameters for experiments presented in figure 11.
Furthermore, the existing experimental data (apart from a small number of outlying points) is well bounded by the inertio-capillary limit presented earlier. This curve thus appears to define a universal upper bound on the coefficient of restitution for a droplet rebounding from a deep pool of the same fluid as a function of 
${We}$, regardless of any other parameters.
There is substantially less data available on contact times for low 
${Oh}$ impacts, but when reported they generally take values within the range of 
$4\unicode{x2013}6 t_{\sigma }$ (Zhao et al. Reference Zhao, Brunsvold and Munkejord2011; Moláček & Bush Reference Moláček and Bush2013; Wu et al. Reference Wu, Hao, Lu, Xu, Hu and Floryan2020).
6. Discussion
In this work we have used a combination of custom experiments, state-of-the-art DNS techniques and a new quasipotential model to study the bouncing of millimetric droplets on a deep bath of the same fluid in the inertio-capillary regime. Weakly viscous models for the bath and droplet interfaces are coupled to one another through the use of a simplified kinematic matching condition, and allow us to make accurate predictions for the droplet trajectory and time-dependent droplet and bath shapes. Furthermore, the quasipotential model is relatively efficient to compute and uses only standard off-the-shelf algorithms, resolving multiple bounces in less than five minutes on a standard desktop computer. Both the quasipotential model and DNS are demonstrated to accurately predict the coefficient of restitution, contact time and maximum surface deflection, validated by direction comparison with experiments with two different working fluids.
Starting from the inertio-capillary limit (where both gravitational and viscous are negligible) in the model, as we increase 
${Bo}$, we see a decrease in the coefficient of restitution and an increase in the contact time. Additionally, a local maximum develops in 
$\alpha$ at low 
${We}$ as 
${Bo}$ increases. As 
${Bo}$ increases further, droplets cease to return to their original height. Furthermore, as 
${Oh}$ is increased away from the inertio-capillary limit, a simpler, monotonic decrease in 
$\alpha$ is observed, with the contact time remaining almost unchanged for the much of the 
${Oh}$ range explored in this work. By further comparison with DNS, the complete model is shown to hold in the limit of small 
${Bo}$ and 
${Oh}$, and up to intermediate 
${We}$, provided that the influence of the intervening gas layer that inhibits coalescence on the overall dynamics is minimal. These trends can be rationalized using simple scaling arguments, with gravity resulting in additional droplet–bath deformations (and thus energy transfer), and viscosity providing a mechanism for energy dissipation in the fluid during contact. Additionally, the model can be used to connect much of the existing experimental data on this particular topic. In particular, the inertio-capillary limit appears to define an upper bound on the possible coefficient of restitution for droplet–bath impact, with the value depending only on the Weber number, and saturating to a near constant value at intermediate 
${We}$.
The related problem of a droplet impactor rebounding off a solid surface has been considered in numerous previous works (Anders, Roth & Frohn Reference Anders, Roth and Frohn1993; Richard & Quéré Reference Richard and Quéré2000; Richard et al. Reference Richard, Clanet and Quéré2002; Gilet & Bush Reference Gilet and Bush2012). The dependence of the coefficient of restitution on the Weber number has previously been reported (Biance et al. Reference Biance, Chevy, Clanet, Lagubeau and Quéré2006; Aussillous & Quéré Reference Aussillous and Quéré2006; Gilet & Bush Reference Gilet and Bush2012), and the trend observed in these studies is similar to what is found in this work for low 
${We}$ and low 
${Bo}$; however, the typical values of restitution coefficients in these studies are significantly larger. This is likely due to the fact that a large portion of the initial droplet energy in the present case is carried away by surface waves excited in the fluid bath. Our general findings also have many similarities with the investigation of Galeano-Rios et al. (Reference Galeano-Rios, Cimpeanu, Bauman, MacEwen, Milewski and Harris2021), in which non-wetting spheres impact and rebound from a water bath. In particular, the general trends for maximum penetration depths and contact times are consistent with the present work. However, spheres with density most similar to that of water show a consistent monotonic increase in the coefficient of restitution with increasing 
${We}$ rather than saturating to a near constant value for the case of droplet–bath rebound. Furthermore, at intermediate 
${We}$, coefficients of restitution can take values as high as 
$\alpha \approx 0.5$, distinctly greater than otherwise equivalent droplet–bath rebounds considered in the present work. Evidently, the nature of the impactor and substrate influences the subtle energy transfer mechanisms across these different capillary rebound problems.
The model presented here is highly versatile, with only a single embodiment thereof considered here in terms of target canonical physical scenario. Future work will consider the effect of relative surface tension and viscosity, where the droplets are composed of a different fluid than the bath. Also, in the present work we have specifically selected the size of the bath to be much greater than that of the droplet radius to eliminate any possible wave reflection and interaction effects. In Zou et al. (Reference Zou, Ren, Ji, Ruan and Fu2013), experimental results indicate that reducing the size of the bath (such that the impacting wave on the bath has time to travel to the edge and return to the impact point) can increase the coefficient of restitution with the droplet recovering energy initially lost to waves. Furthermore, the effects of incident droplet or bath deformations could be readily studied, and has been shown to influence bouncing behaviour in similar systems (Biance et al. Reference Biance, Chevy, Clanet, Lagubeau and Quéré2006; Yun Reference Yun2018).
There are a number of possible additions to the existing model that could expand its reach to other related problems. For instance, the model for the droplet deformation can be extended into a regime where the dynamics of the gas layer does matter, and the gas layer dynamics coupled to the droplet deformation through the use of lubrication equations. The model can also be adapted to non-axisymmetric domains, or to droplet impacts at varying angles of incidence. However, in these cases the full kinematic match should be utilized, as the shape of the pressure distribution would likely change at each time step. Moreover, numerous authors have studied the variety of phenomena that occur when a droplet impacts another droplet (Qian & Law Reference Qian and Law1997; Tang, Zhang & Law Reference Tang, Zhang and Law2012). Droplet–droplet collisions are of extreme importance in combustion science (Jiang, Umemura & Law Reference Jiang, Umemura and Law1992) and cloud formation (Grabowski & Wang Reference Grabowski and Wang2013), for instance, the general effect of cloud turbulence acts to increase droplet–droplet interaction, and droplet impact and coalescence is postulated as the primary mechanism by which warm rain forms (Grabowski & Wang Reference Grabowski and Wang2013). With such motivation in mind, the present model could be readily extended to cases where equal and unequal sized droplets impact and rebound from one another. Overall, the quasipotential model developed in this work has the potential to continue to inspire and inform the rich subject of capillary rebounds.
Supplementary movies
Supplementary material (experimental, model and simulation animations) are made available to the interested reader.
Supplementary movies are available at https://doi.org/10.1017/jfm.2023.88.
Acknowledgements
The authors would like to acknowledge C. Galeano-Rios for his support and guidance throughout the project, A.H. Kumar for helpful discussions, F. Blanchette for sharing code associated with prior work and J. Moláček for sharing experimental data from prior work. The Imperial College London Research Computing Service and the Scientific Computing Research Technology Platform at Warwick have both supported this work with computational resources. All authors would like to thank the referees for their constructive suggestions.
Funding
The authors gratefully acknowledge the financial support of the National Science Foundation (NSF CBET-2123371) and the Engineering and Physical Sciences Research Council (EPSRC EP/W016036/1). R.C. also thanks the London Mathematical Society for early stage funding through a Scheme 4 Research in Pairs Grant (Ref. 41813).
Declaration of interests
The authors report no conflict of interest.
Code accessibility
Codes and associated documentation can be found at https://github.com/harrislab-brown/BouncingDroplets for the mathematical model implementation (MATLAB) and https://github.com/rcsc-group/BouncingDroplets for the DNS framework (Basilisk C) used as part of this work.
Appendix A. Influence of ambient gas properties
In order to verify our assumption that the flow within the air layer is negligible to the droplet and bath dynamics in our parameter regime of interest, we run DNSs where the ambient gas density and viscosity are varied. First, the gas viscosity is held fixed for air at 
$21\,^{\circ }$C and 1 atm and the density is increased by a factor of four, and then decreased by a factor of four, respectively. Then the variation process is repeated for the gas viscosity, with density held fixed. The results of these simulations are presented in figure 12 and are nearly indistinguishable, particularly during contact.
Figure 12. The DNS predictions of the trajectory of a deionized water droplet with 
$R=0.35$ mm in air at 
${We}$ 
$=0.7$ (with 
${Bo}$ 
$=0.017$, 
${Oh}$ 
$=0.006$). The ambient gas density and viscosity are increased and decreased independently by a factor of four. These simulations are compared with the case with the reference case of air properties at STP.
Appendix B. Influence of pressure shape function
We tested several pressure shape functions 
$H_r$ to assess the relative influence of this choice. In the simulations shown in figure 13 a parabola, fourth-order and sixth-order polynomial with time varying contact areas set as described in the modelling section are compared. A parabola with fixed contact radius 
$r_c(t) = R$ is also included in the comparison. These results all correspond to 
${Bo} = 0.017$ and 
${Oh} = 0.006$ impacts. The parabola with the fixed contact area performs the most poorly, especially at lower impact 
${We}$. At higher impact 
${We}$, the fixed radius 
$r_c(t)=R$ parabolic distribution prediction becomes similar to the time evolving parabolic case. For the simulations presented in this work, the contact radius is defined to have a maximum value of 
$R$, as the projection of the pressure distribution onto the undisturbed spherical surface is no longer well defined for 
$r_c > R$. The value for 
$r_c(t)$ quickly saturates to 
$R$ at higher impact 
${We}$, and the agreement between the constant contact radius and the contact radius model used in the present work improves. Overall, inclusion of a time-varying contact radius appears necessary to capture the correct trends in 
$\alpha$ and 
$t_c/t_{\sigma }$ over the range of 
${We}$ presented in this work. In addition, as the order of the polynomial increases, the shape of the pressure function more closely resembles a top hat, and the predictions of the model improves as compared with the corresponding DNS. The sixth-order polynomial was ultimately chosen for the present work, as higher polynomials (corresponding to broader flat regions) converged more slowly with only marginal changes in the quantitative predictions.
Figure 13. Quasipotential model and DNS predictions for a deionized water droplet with 
$R=0.35$ mm corresponding to 
${Bo} =0.017$ and 
${Oh} =0.006$. Droplet centre of mass trajectories for (a) 
${We} =0.07$, (b) 
${We} =0.36$, and (c) 
${We} =1.36$. (d) Coefficient of restitution 
$\alpha$ and contact time 
$t_c/t_{\sigma }$ as a function of 
${We}$. (e) Plots of the pressure shape functions 
$H_r(r/r_c)$ tested in this figure, shown for reference.
Open access
