Non-wetting impact of a sphere onto a bath and its application to bouncing droplets

We present a fully predictive model for the impact of a smooth, convex and perfectly hydrophobic solid onto the free surface of an incompressible fluid bath of infinite depth in a regime where surface tension is important. During impact, we impose natural kinematic constraints along the portion of the fluid interface that is pressed by the solid. This provides a mechanism for the generation of linear surface waves and simultaneously yields the pressure applied on the impacting masses. The model compares remarkably well with data of the impact of spheres and bouncing droplet experiments, and is completely free of any of impact parametrisation.


Introduction
The study of impacts of solids onto the surface of a fluid bath is motivated by several applications, including the impacts of sea landing planes, the behaviour of projectiles entering water (Richardson 1948), slamming of boat hulls (Howison, Ockendon & Oliver 2002) and bio-locomotion mechanisms for water striders (Bush & Hu 2006). Interest in the detailed study of these phenomena was sparked by Worthington over a hundred years ago (see Worthington 1882Worthington , 1897. Since then, several works have dealt with cavity formation upon projectile entry, cavity pinch off (Truscott, Epps & Belden 2014) and the forces on the solid as it penetrates the fluid mass (Abelson 1970;Howison, Ockendon & Wilson 1991;Aristoff et al. 2010).
These impacts typically undergo three different stages, e.g. early contact, cavity formation and steady cavitating flow (Logvinovich 1969). The initial stages of penetration are of special importance, since their results influence the outcome of the following phases. Moreover, the geometric properties of the solid body and the physical properties of the fluid are determinant for the dynamics of the initial stages of these impacts (Korobkin & Pukhnachov 1988); thus, different approaches have to be developed accordingly.
A model for the early stages of the impact of a blunt body onto a bath of an incompressible fluid was developed in Wagner (1932). The Wagner problem is stated in terms of the linearised free-surface boundary conditions of an ideal fluid, subject 98 C. A. Galeano-Rios, P. A. Milewski and J.-M. Vanden-Broeck to a continuity constraint for the elevation of the fluid interface. The free surface is decomposed into two sections, one which corresponds to the part that is in contact with the solid, which is assumed to match the solid shape, and the rest of the free surface where the velocity potential is assumed to be null (Korobkin 2002).
The Wagner problem, even though it deals with a set of linear differential equations, still poses an important challenge since the domain has to be divided into two parts whose extent is to be determined (Korobkin & Pukhnachov 1988). That is to say, the regions where the different boundary conditions need to be applied are an unknown in the problem.
Wagner's approach was later expanded to cover the axisymmetric case of the entry of a solid of revolution (Schmieden 1953). Tri-dimensional cases were studied in Scolan & Korobkin (2001), Korobkin (2002) and Korobkin & Scolan (2006). All these works assume that the fluid is ideal, the flow is irrotational and that surface tension and gravitational effects are negligible. More recently, the work of Lee & Kim (2008) explored experimentally and theoretically the case of super-hydrophobic balls impacting the surface of a bath. They develop a model for the forces of an impacting ball, their model includes nonlinear free-surface elevation and hydrostatic forces. They show evidence of capturing the initial stages of impact and part of the reversed motion during the bounce off, however their model predictions fail to capture the intermediate stage between this two, during which the velocity of the ball reverses.
The present work is concerned with the study of impacts that do not break though the surface, and that might even cause the solid to bounce off, which can be thought of as the solid not surpassing the initial stages. Thus, this work has many characteristics that are similar to the models for the early phases of more general impacts, however it includes effects that are not typically accounted in these models.
We consider the solid as perfectly hydrophobic, which in practice means that we take the contact angle to be π at all times, ignoring any elaborate form of contact angle dynamics. In reality, there might be a very thin layer of air separating the two masses at early stages of 'contact'. This layer might take a finite time to be drained and pressure on the fluid and the ball could in principle depend on the air dynamics (Purvis & Smith 2004) . Here, we assume that the solid is convex and we ignore any lubrication layer dynamics of the air, i.e. we simply assume that this layer transmits the pressure from the fluid bath to the solid. We also assume smoothness of the solid to guarantee the existence of a well-defined value of curvature at every contact point.
In contrast with the works mentioned above, we do not assume that the fluid is ideal, instead we use a weakly dissipative model of free surface flow and include gravity and surface tension effects. We take special care in approximating surface tension effects as high curvatures are typical in small impacting bodies and interfacial tension plays an important role in their bounce off.
1.1. An application to bouncing droplets Impacts of droplets on a free surface have recently been the object of revived interest due to the relatively new discovery (Couder et al. 2005a,b) that when a fluid bath is subject to vertical oscillations below the Faraday threshold, it is possible to have a droplet of the same fluid sustain periodic oscillatory motion as it bounces on the free surface of the bath. The inhibition of coalescence is made possible by the sustenance of a thin layer of air that is replenished before van der Waals forces can induce coalescence (Couder et al. 2005a;Terwagne, Vandewalle & Dorbolo 2007).
Repeated droplet impacts trigger waves on the free surface that in turn influence the future impacts of the droplet. Wave and droplet synchronise to display different periodic modes and even a chaotic mode of bouncing, typically observed in the vicinity of the Faraday threshold Eddi et al. 2011;Wind-Willansen et al. 2013;Moláček & Bush 2013a;Bush 2015). The bouncing droplet and its accompanying wave field form a coherent association which has come to be known as a bouncer (Moláček & Bush 2013a).
More surprisingly, a bouncer can become unstable to lateral perturbations and, consequently, initiate a horizontal trajectory, i.e. it becomes a walker (Moláček & Bush 2013b). In this regime, the horizontal trajectory of the droplet is influenced at every bounce by the shape of the free surface that it impacts. This behaviour gives rise to a very complex and interesting dynamical system, which in fact presents many features that are reminiscent of another wave-particle association, namely that of the quantum realm (Eddi et al. 2009;Fort et al. 2010;Harris et al. 2013;Harris & Bush 2014;Bush 2015).
Several modelling approaches to describe and predict the behaviour of rebounding droplets in varied scenarios have been developed. Many of these models (Eddi et al. 2011;Moláček & Bush 2013a,b;Oza, Rosales & Bush 2013) estimate the resulting free-surface elevation by the superposition of Bessel functions, or an approximation of them, whose amplitudes depend on time. Other models solve the free-surface evolution, imposing the effects of a droplet impact by means of a pressure field on the interface through different approximations (Milewski et al. 2015;Faria 2016;Durey & Milewski 2016); however, none of these models uses the natural, geometric and kinematic, constraints to calculate the evolution of the system during impact. There is a computation of the contact area in Milewski et al. (2015), which is limited to a linear estimate based on droplet penetration and it is not required to satisfy the small-scale geometric restrictions during impact. Blanchette (2016) imposes a geometric constraint at the south pole of the droplet, from which the force value is derived; however, there is no pressure distribution or extent of the contact area explicitly calculated.
Here, we apply our impact model to the study of bouncing droplets, since in this case there is a thin layer of air separating the droplet from the bath, and droplet deformation can be neglected for small drops. Thus the drop is treated as a solid ball with tangential 'contact' with the fluid.

Formulation
We first present a linear water wave model that incorporates viscous damping effects, along the lines developed by Lamb (1895) and Dias, Dyachenko & Zakharov (2008). We then introduce a model for the generation of waves through the local deformation of the free surface due to an impacting smooth solid body.

Fluid equations
We consider the three-dimensional free-surface incompressible flow of a fluid bath of uniform density ρ, subject to gravitational and viscous forces. The domain is unbounded in all horizontal directions, is of infinite depth and the fluid is initially at rest with its free surface undisturbed. We introduce Cartesian coordinates with the positive z axis pointing upwards. We also assume that the free surface can be described by a well-defined smooth function z = η(x, y, t). Under these assumptions 100 C. A. Galeano-Rios, P. A. Milewski and J.-M. Vanden-Broeck the flow is governed by the Navier-Stokes equations for an incompressible fluid, i.e.
where u = (u 1 , u 2 , u 3 ) T and p, are the velocity and pressure field, ν is the kinematic viscosity and g = (0, 0, −g) is the acceleration due to gravity. Equations (2.1) and (2.2) are subject to the condition of no motion at infinity, i.e.
( 2.3) The interface is subject to surface tension and normal stresses due to the air pressure. The stress balance and the kinematic condition require, respectively, that (2.5) Here τ = ρν(∇u + ∇u T ), p s is the pressure of the air just above the free surface, n = [−η x , −η y , 1] T / 1 + (η x ) 2 + (η y ) 2 is the outward pointing unitary normal to the free surface, σ is the surface tension coefficient and κ[η] is the sum of the signed principal curvatures of the free surface, which is positive for convex functions. We recall that p =p − ρgz and g = ∇(gz), wherep is the dynamic pressure, and incorporate conservative volume forces in the pressure term. We linearise around the flat surface equilibrium η(x, y, t) = 0, which yields (2.6f ) We look for solutions that satisfy the Helmholtz decomposition where v = ∇φ and w = [w 1 , w 2 , w 3 ] T = ∇ × Ψ , with φ and Ψ being the scalar and vector velocity potentials, respectively. To satisfy condition (2.3), we assume that for all t 0 φ, ∇φ → 0 and Ψ , w → 0, when x 2 + y 2 + z 2 → ∞. (2.8a,b) Substituting equation (2.7) into (2.6b) and into (2.6a) we obtain respectively. We note that (2.9b) is satisfied if we choose Moreover, substituting equations (2.7) into (2.6c), (2.6d) and (2.6e) we obtain We take the x derivative of (2.11a) and the y derivative of (2.11b), add the results and use (2.10b) for i = 3 to obtain where ∆ H = ∂ xx + ∂ yy and we used w 1 x + w 2 y + w 3 z = 0. We now combine equations (2.10a), (2.9a) and (2.11c) to obtain We take ρ as the characteristic density, define the characteristic length L, velocity V and introduce the dimensionless numbers Fr = V 2 /(gL), We = ρV 2 L/σ and Re = LV/ν, (2.14a−c) i.e. the square Froude number, Weber number and Reynolds number, respectively. Then from (2.9a), (2.10b), with i = 3, (2.6f ) and (2.7), (2.13), (2.12) and condition (2.8) we obtain a closed system which is given in dimensionless form by subject to φ, ∇φ → 0 and w 3 → 0, when x 2 + y 2 + z 2 → ∞. (2.16a,b) We also note that from (2.15c) and (2.15e) we have w 3 t = 2∆ H (η t )/Re, at z = 0. Since the fluid is initially at rest and its free surface is undisturbed, this implies We substitute equations (2.7) and (2.17) into (2.6f ), and we obtain We also note that (2.17) implies that w 3 = O(Re −1 ) near the free surface, and (2.15b) implies that the boundary layer thickness scales as Re 1/2 . Thus, we rescale equation (2.15d) more properly as Finally, we recall our high-Reynolds-number assumption and disregard the term of highest order in (1/Re) in (2.19). Combining this with (2.15a), (2.18), and the condition (2.16) we obtain the dimensionless system (2.21) 2.1.1. Reducing the fluid system to the boundary We note that (2.20b) and (2.20c) require the evaluation of function φ and φ z only on the boundary plane z = 0. φ z is the normal derivative at the boundary, of the solution of the Laplace problem in the half-space with the decay conditions (2.21). We can write 22) where N is the linear Dirichlet to Neumann map, which is well defined in this domain for a smooth function φ that decays sufficiently fast at infinity. In fact, its expression is given by the principal value of a singular integral, namely where r = (x, y), s = (x , y ) and B(r; ) = {s = (x , y ), |r − s| < }.
In appendix A, we prove that under our working assumptions a function φ with a relative weak decay must satisfy equation (2.23). We also show that the singular integral in (2.23) converges wherever φ is smooth.
The considerations above reduce the problem of calculating the surface evolution to solving (2.25b) where m is the mass of the solid; A C is the normal projection, on the horizontal plane, of the portion of the solid that is in contact with the surface (see figure 1) and c f incorporates the effects of friction of the ball with the air. In general, the dependence of c f on the velocity of the solid, proximity to the free surface and other physical parameters can be very elaborate (Brenner 1961). In many cases it can be disregarded altogether or approximated simply, such as by using Stokes' drag (Moláček & Bush 2013a) for the flight of bouncing droplets. Here, for consistency with prior work, despite its effect being small as c f V/mg 1, we use Stokes' drag in the bouncing droplet case, and set c f = 0 in the solid impact case, as the sphere's velocity is imposed at impact. Similar equations can be formulated for horizontal motion of the solid, as well as for its rotation. However, at this stage, we decide to limit our study to the case of axisymmetric solids impacting on axisymmetric free surfaces. This removes the need to calculate the rotation and the horizontal motion of the solid. Moreover, the height of the lowest point of the solid body h, which lies on the symmetry axis (as a consequence of symmetry and convexity, see figure 1) is given by a vertical translation of h c .
We are imposing that the pressure applied on both impacting surfaces is the same; this means that, if there exists a thin air layer separating the two, the pressure across it does not change, consistent with lubrication theory. We will also assume that air pressure above the surface p s (r, t) is only different from zero on A C . We write equation (2.26) for h and in dimensionless form and obtain where D = c f L/(mV 2 ) and M = m/ρL 3 . We assume that the bottom part of the solid body is described by a function z s , with smooth curvature in the vicinity of the impacting region. We recall that the problem is axisymmetric, and thus we define a radial coordinate system given by r = x 2 + y 2 . Since functions of (x, y) are simply functions of r, the main constraint for the fluidsolid interaction can then be stated as (2.28) which must hold everywhere under the solid, and where we assume z s (0) = 0. We impose a second natural constraint, namely that the contact angle, at the boundary of the pressed surface S C (see figure 1), has to be π. This assumption is equivalent to stating that the effect of surface tension at the boundary of S C is exactly equal to the effect of the jump in pressure due to the curvature of S C . This can be seen using the method presented in Keller (1998), with the difference that integrations in this case need to be carried out in the contact area and not outside of it.
We ignore the dynamics of air and thus identify A C as the region of the z = 0 plane where relation (2.28) is satisfied as an equality. In practice, this means that there will be no distinction between the pressed part of the free surface and that of the solid. Figure 1 might induce the reader to think that there in fact exists a difference in height between the two pressed surface portions, however this is not the case in the model. The separation shown in figure 1 serves merely didactic purposes.
We thus wish to solve equations (2.24) and (2.27) subject to conditions (2.25), (2.28) and that the surface is tangent to the solid at the contact line. Based on our symmetry and convexity assumptions, we have that A C must be a disc. An important particularity of the system given by (2.24) and (2.27) is that we do not have a law for the evolution of A C , i.e. we do not have an a priori estimate of where to apply pressure at a future time. This scenario suggests a method that iterates on the radius of A C and verifies the constraint (2.28) and the tangency condition on the boundary of A C . Moreover pressure values, to be calculated across A C , must be such as to enforce that the condition (2.28) is satisfied in the form of an equality over A C . However, η and h are unknown; strongly suggesting the need for an implicit iterative method.
We have purposely not linearised the curvature of the free surface κ[η] in (2.24b). Under the solid we shall use the full curvature. This does not pose a mathematical difficulty, as a jump in curvature is expected at a contact line. Further, the resultant force on a wetted solid can be computed as an integral along the contact curve. This integral, as the contact angle tends to π, converges to the Young-Laplace pressure jump integrated over the wetted area (Keller 1998). We linearise the curvature on the free surface, as this removes the need to implement a nonlinear solver in the numerical scheme. In what follows, we assume κ[η] ≈ ∆ H outside of A C , and the exact curvature of the sphere (κ = 2/R o ) on the contact area.
The problem is now reduced to solving where all spatial functions are radial, r c is to be determined, R o is the radial extent of the axisymmetric solid and where we have separated the curvature operator κ into its linear (∆ H ) and its nonlinear part.

Numerical implementation for axisymmetric impacts
We look for approximate solutions of the system (2.29) on an evenly spaced radial mesh. We thus need to find discrete approximations of the ∆ H and N operators, and of the integral in (2.29c), applied to radial functions. Details of the derivation of the matrix representation of those operators on a regular mesh are presented in appendix B. We can now use an implicit Euler scheme in time to express the discrete version of the problem as where, in turn, I is the identity matrix, δt is the time step, A stands for the linear functional defined by the integral and ∆ H and N now stand for the discrete approximation of the operators they denoted above. The functions in (3.3) and (3.4) are the row vectors of discrete approximations to the functions they represent, and super-indexes indicate discrete values of time. We note that Q is a (3n r + 2) × (2n r + 2) matrix, where nr + 1 is the number of radial points in the mesh, and the outermost point is assumed to take zero value at all times as a consequence of the decay hypotheses. However, if we define k as the number of contact points of the mesh, we have p j+1 s (iδr) = 0 for i k, since we number the radial mesh points starting at the origin. This implies that in the last n r + 2 − k columns of Q, only the last 2 are relevant to the system. Moreover, we note that the first k components of vector η j+1 contain the vertical coordinates of the contact points of the solid, since the free surface matches their height at those points. Thus, we have Substituting this in (3.1), we obtain and where the super-index k on a matrix refers to the sub-matrix formed by its first k columns, and the super-index k to the sub-matrix formed its last n r − k columns. Moreover, κz s denotes the row vector of discrete curvature values of z s at the mesh points, and e k+1 is the (k + 1)th canonical row vector, which is used in (3.9) and (3.11) to account for the contribution of the outermost contact point to the computation of ∆ H η j+1 (kδr); see (B 4) in appendix B. We note that, provided we have determined the number of points of contact between the solid and the free surface, equations (3.7)-(3.9) and (3.5) yield a closed system for a recursion law for η, φ, h and h t . Thus, we can calculate the solution for k = 1, 2, . . . , k max contact points along the radius of A C , where k max is given by the extent of the solid, discard solutions that violate restriction (2.30a) on the mesh points, and among the remaining ones we can choose the one that minimises the tangency error at r = r c . This approach includes no assumptions on the temporal evolution of the S c ; however, it is reasonable to expect that the contact lines changes continuously, which suggests testing only the vicinity of the contact area used in the previous time step.
We note that our numerical method is first order in time. This was found to be sufficiently accurate for the applications presented in the following sections provided an adaptive time step is used.

Simulations
We present two instances of the implementation of the model described here. First, we study the impact of a solid sphere onto a quiescent free surface. Second, we use this model in the context of bouncing droplets, i.e. we use this interaction model to compute the repeated impact of a droplet onto a vibrating bath, assuming that the droplet does not deform nor coalesce and can therefore be well approximated by a rigid sphere. Visualisations of the results are also provided in the additional material.

Solid spheres impacting a quiescent bath
We consider a perfectly hydrophobic solid sphere of density ρ s , impacting normally onto the free surface of a quiescent bath and assume that the drag due to the air as the ball impacts the bath is negligible when compared to the other forces at work during impact, i.e. D = 0 in (2.27). Thus, defining L and V, in (2.14), as R o (radius of the solid sphere) and V 0 (incoming velocity), the four dimensionless numbers that characterise the system, are We choose the radial discretisation so as to guarantee at least 40 mesh points over the unit length of the ball's radius and a domain of computation greater than 10 such units and verify that during the numerical experiments the waves are negligible near its boundary. We use cubic interpolation for the radial function φ to estimate operator N. We initially test the method by solving for all possible discrete contact lengths at each time step, i.e. taking k = 0, 1, 2, . . . , k max , with k max equal to the integer part of R o /δr + 1. We accept as the solution, among those which satisfy condition (2.28), the one which minimises the tangency error in absolute value. Figure 2 illustrates this point; figure 2(a) corresponds to using r c too small, which typically produces a small tangency error but violates condition (2.28) (note that the fluid is overlapping with the sphere). Figure 2(b) corresponds to satisfying condition (2.28) and panel 2(c) also satisfies (2.28) but with a much larger mismatch in tangency. The minimum absolute tangency error is usually achieved by the smallest r c that satisfies condition (2.28).
We notice that successive reductions of the time step indicate that the change of r c in time is continuous; never skipping a neighbouring discrete value, provided the time step is sufficiently small. The evidenced continuous nature of r c in time, together with the fast changes in velocity observed for stronger impacts, induce the choice of an adaptive time step method. Moreover, recording the error in tangency at each tested value of r c provides evidence that the tangency error behaves monotonically with the error in r c . That is to say that, for radii that are smaller than the optimal, the surface of the fluid has a greater slope than that of the solid ball at r c , while for larger radii the situation is reversed. Figure 2 illustrates this situation. Based on the monotonicity of the tangency error, and on its continuity in time, we can test at each time step the previous value of r c , then test a neighbouring value to identify the direction of error decrease, and test values of r c simply following the gradient. If the minimum of the error is farther than δr away from the previous value, we repeat the calculation with a δt that is halved. We do this recursively until we are able to track the motion of the contact front to increases of one mesh intervals at each δt. We thus need only test at most 4 possible values of r c at each attempt to use a given time step before reducing it if necessary (two tests to determine direction of improvement and at most two more to verify that the next one in that direction is a local minimum of error). We impose an extra condition on the time steps, if the previous successful time step used was smaller than an initially prescribed value, for the next time step we will use, at most, twice its value. This allows for a relatively smooth variation of time steps, while we also require that a maximum time step is not surpassed.

Comparisons to experiments
Experimental data on solid balls covered with a superhydrophobic coating are reported in Lee & Kim (2008). Comparisons between our simulations and these experimental results support the validity of our model. Figure 4 of Lee & Kim (2008) presents a vertical tracking of the centre of a hydrophobic ball as it impacts on the free surface of water. Figure 3(a) compares our simulation results for the same configuration. The solid curves correspond to the simulation of the same experiment assuming standard values of physical properties of pure water at 25 • C (see appendix C). Figure 3(b) shows the temporal evolution of the total upward force on the droplet during impact, as well as its comparison with the fraction of this force that is due to surface tension. Figure 3(b) is consistent with the notion that inertia is the main force during the first stages of impact, as is assumed in the Wagner model, and also shows that at the scale we are working, surface tension becomes more important as the impact continues. It is, in fact, the most relevant effect during lift off, since removing it would reverse the sign of the force. Forces due to hydrostatic effects are small, in comparison, throughout the whole impact. It is important to highlight that the vertical trajectory obtained from the experiments reported in Lee & Kim (2008) deforms the surface substantially (over a full diameter), and hence one might expect a model based on a linearisation of the free surface to differ substantially from experimental results. This is certainly not the case, as shown  in figures 3(a) and 4, where we see that agreement is excellent, with the possible exception of lift off behaviour. This discrepancy could be due to a cumulative effect which results from linearising the fluid equations.

Bouncing droplets
We model the successive impacts of a bouncing droplet using the method presented here. To this end, we treat the drop as a rigid sphere, and apply at each impact the method described in the previous sections and the appendices. However, some adaptations are required. First, we follow Milewski et al. (2015) in introducing a viscosity correction to the fluid model (ν * = 0.8025ν), which essentially adjusts dissipation so as to match the experimental Faraday threshold of the low viscosity silicone oil experiments that we model. Second, characteristic length and time scales of the fluid motion are no longer set by the size of the drop, instead they are determined by the scale of the Faraday waves, which are subharmonic to the forcing. Namely, L = λ F and V = λ F f F ; wavelength and phase velocity, respectively. Here we use λ F = 0.4969 cm, which corresponds to the corrected viscosity model for this fluid equations (Milewski et al. 2015). This yields Moreover, we solve the system in the non-inertial frame of reference of the shaking bath, this implies we must introduce a time-dependent oscillatory term to the gravity field. Thus, we multiply the dimensionless number 1/Fr in (2.29b) and (2.29c) by (1 − Γ cos(ω 0 t + θ 0 )).

(4.2)
We also include friction due to air, which we approximate using Stokes' drag on a sphere, and change the size of the radial mesh intervals to a tenth of the droplet radius, as the need to calculate repeated impacts would turn our original mesh excessively time consuming. Another reason to accept this coarser mesh is simply that in this problem deflections are much smaller, consequently fewer points produce a good approximation of the surface shape. In contrast, the domain size in this problem needs to be large enough as the standing waves are spatially extended. Consequently, we define our circular domain of computation to be of 10 Faraday wavelengths in radius, this choice proved to be sufficient to produce waves that vanish at the boundary for all regimes presented in the work of Milewski et al. (2015). In order to produce simulation results that can be compared later to experimental results, we simulate the configuration of 20 cSt silicone oil and 80 Hz shaking used in Wind-Willansen et al. (2013), i.e. Fr = 0.81, We = 9.04 and Re = 61.54; more details on this configuration are given in appendix C. We run simulations that correspond to 240 forcing periods (T f ), which guaranteed attaining a steady regime. We follow Gilet & Bush (2009) in their notation of bouncing modes defined by ordered pairs (m, n), in which m describes the period of vertical motion in units of T f and n stands for the number of contacts intervals that droplet and bath sustain in one period of vertical motion. They also introduce two different (2, 1) modes, namely (2, 1) 1 and (2, 1) 2 , which are qualitatively different. Following Dorbolo et al. (2008) and Wind-Willansen et al. (2013), we also classify our droplets in terms of their vibration number Ω = ω 0 ρR 3 o /σ . In practice, Ω is a proxy for droplet radius R o , since the other variables that define Ω (fluid properties and forcing frequency) are fixed by the experimental set-up.
With this model we are able to reproduce most features of the bouncing droplet phenomena. Surface waves are triggered by the successive impacts of the droplets, and experimentally observed modes of bouncing arise spontaneously, without resorting to the parametrisation of impact used in Moláček & Bush (2013a) and Milewski et al. (2015). Figure 5 illustrates the details of a characteristic impact of the bouncing droplet on to the fluid surface. This figure shows the great difference in curvature between the free surface and the pressed surface S C . Figure 6 shows this more emphatically. In this figure, the vertical scale is exaggerated with the purpose of making small surface waves visible. On the scale of the waves the curvature of the droplet appears as a sharp cusp, see figure 6(a), this finding is consistent with observations made in previous studies, see Gilet (2014). In figure 6(b), the highlighting of the contact area illustrates the motivation to use the full nonlinear curvature where the droplet presses against the free surface, and its linear approximation everywhere else. Modes of bouncing are often identified in experimental work by taking a single pixel vertical slice of the high-speed videos and placing them side by side from left to right in order of increasing time, which facilitates somewhat the identification of contacts (Protière, Boudaoud & Couder 2006;Wind-Willansen et al. 2013;Terwagne et al. 2013). The analogue plot in this model consists of the trajectory of the south pole of the droplet and the central point of our free surface. Figure 7 shows characteristic plots for some of the bouncing modes we have found.
This model predicts the pressure field under the drop during contact. Integrating this field in space we obtain the value of the force that is applied to the droplet. Figure 8 shows the prediction of the force on the droplet throughout two forcing periods. We notice that, as expected, surface tension contributes a substantial amount to the vertical impulse on the droplet. Figure 8 also shows a comparison with the predictions of the parameterised model presented in Milewski et al. (2015), which is based on the spring model developed in Moláček & Bush (2013a). We note that there is an mistake in an equation used in Milewski et al. (2015), carried over from a typographical error in (3.7) and (A 14) of Moláček & Bush (2013a). In the denominator of the coefficient of the damping term of (2.42) in Milewski et al. (2015), the logarithm should be squared. Thus, the equation should read  . We recall this spring model was derived assuming a constant contact area underneath the droplet and without a kinematic constraint imposing that the surface deflection had to be consistent with the motion of the droplet (although the behaviour, verified a posteriori was reasonable). Both methods predict a fast rise in the forces at the initial stages of contact, followed by a slower decay until take off. We also show the decomposition of the total force on the drop, providing insight on which forces are more important at each stage of impact (inertial versus surface tension). As in the case of the fast impact of a solid ball, discussed in the previous section, both inertia and surface tension play a significant role. At the initial stages of impact, in the cases considered; inertia contributes approximately as much as surface tension, whilst at later stages capillary effects are dominant. Moreover, figure 8 provides evidence on how the impacts change as Γ increases (also seen in figure 11). The transition from the (1, 1) to the (2, 2) mode takes place when the first impact occurs later and the second impact earlier, breaking the (1, 1) mode symmetry. Continuing this process, the first impact becomes less important and the second one more until they merge, giving rise to a stronger single impact in the (2, 1) mode. Figure 9 shows the evolution of the contact area in the same cases for which the force was depicted in figure 8. Since the curvature of the sphere is constant, the contact area determines the contribution of surface tension to the vertical force on the drop (shown in solid grey lines in figure 8).  (2, 2), (c) a slightly different (2, 2) mode, (d) a (2, 1) 1 mode, (e) a (2, 1) 2 mode and ( f ) a chaotic mode.

Comparisons to experiments
We map the steady state regimes that we observe in the simulations on the Γ -Ω plane against the corresponding experimental data points found in Wind-Willansen et al. (2013), see figure 10. We include a hashed area in this figure, which corresponds to the mode shown in figure 7(c). This mode of bouncing, whilst according to the definition is a (2, 2) mode, might easily be mistaken for a (2, 1) 1 mode, since the two contacts are only briefly separated by a short time in which the droplet hovers  Damiano (2015) show that most bouncers that correspond to our hashed region are identified as being in mode (2, 1) 1 in the former and (2, 2) in the latter.
Considering the difficulty to distinguish these cases experimentally, the agreement is excellent. As mentioned above, the values of Ω were swept by changing droplet radius in each simulation. However, we also computed certain cases with a different forcing frequency, for the same fluid in the bottom right region of figure 7, which produced a (4, 2) mode, as is seen in experiments (white squares in figure 10).
The other main discrepancy between theory and experiment that can be seen in this figure is the value of Γ at which we observe the onset of chaotic bouncing is higher than experimentally observed values. Chaotic bouncing occurs for large Ω and the difference can probably be explained by the importance of droplet deformation for larger drops which is not accounted for in our simulations.
For bouncers at Ω = 0.8, we trace the dependence of the phase of shaking corresponding to the impact and release of the droplet. We can thus compare the results to the experimental observations in Damiano (2015) and the computations of Milewski et al. (2015). The background shadings in figure 11 correspond to the 95 % confidence interval of phase of impact and release found experimentally (A. P. Damiano, Personal communication, 2015). Each vertical traverse in this figure corresponds to a bouncer. At the bottom of the graph, every bouncer starts in flight and touches down as it reaches the first black curve. It then stays in contact until the next curve is reached (blue), at which moment it re-enters flight and repeats the same behaviour. As time increases, this figure is repeated periodically in the vertical direction. In the simulations; for Γ /Γ F = 0.48 we see the transition from the (1, 1) mode to the (2, 2) mode and for Γ /Γ F = 0.71 we see the (2, 2) mode to (2, 1) 1 mode transition. The experimental data in figure 11 correspond to a more recent experimental result which used more reliable techniques to determine contact. In particular, we notice that for Ω = 0.8 and 0.5 < Γ /Γ F < 0.69, Damiano (2015) reports a (2, 2) bouncing mode while Wind-Willansen et al. (2013) reports a (2, 1) 1 mode. This supports our claim that in this regime it is hard to distinguish between the two, which was the base for our choice of the introduction of a hashed region in figure 10. Our simulations show better agreement with the more recent experimental results and better agreement with experiments than the logarithmic spring model.
Finally, we are able to compare our wave field predictions to the experimental observations reported in Damiano et al. (2016), where a surface Schlieren method was used to measure wave topography. Figure 12 shows a comparison of the radial profiles obtained by measurements and simulations. We note that the agreement is reasonable, especially considering the fact that there is no fit in phase. However, numerically computed waves tend to be larger in amplitude. This discrepancy between theory and experiment probably arises from our quasi-potential approximation and the simple form of the dissipation. In this approximation, all of the forcing due to droplet impact is deposited into wave modes, whereas it is probable that in reality non-wave-like vortical motions are also generated and rapidly dissipated, resulting in a smaller total wave response. . The colour coding of the squares is the same as that of the background, with the addition of white for the (4, 2) mode and pink for the (4, 3) mode, which were not found in this sweep of physical parameters with the simulations. The hashed area corresponds to modes that are to be strictly classified as (2, 2); however, the intermediate flight that separates the two contacts is rather subtle, and thus could be easily taken to be a (2, 1) 1 . Walkers were found experimentally to the right of the red curve.

Conclusions
From the free-surface Navier-Stokes equations we derive a linear set of partial differential equations to model the motion of the free surface under a high-Reynoldsnumber approximation and couple it to the motion of the impacting hydrophobic solid through the consistent adjustment of the contact area and matching velocities. From this the pressure on the solid can be obtained, coupling back to its equation of motion. The model does not require explicitly that the impacting object be a sphere; however, calculations of the motion of the fluid and solid are largely simplified in axisymmetric cases since we do not need to calculate rotational solid motion and, provided the surface is axisymmetric initially, the net forces are vertical. In the future we wish to extend the applicability of this model weakening some of these assumptions, in particular for the case of walking droplets.
The resulting problem involves mixed boundary conditions at the surface of the fluid. In order to reduce the usual free-surface potential flow problem to a problem on the boundary alone one can define a Dirichlet-to-Neumann operator that physically corresponds to, given the tangential velocity of the surface, finding the normal velocity. The impact of a solid creates mixed conditions, as on the pressed part of the interface the constraint that the fluid surface coincides with the solid modifies the free-surface boundary conditions. We note that the domain of the impact is, a priori, unknown and needs to be found as part of the problem. Here, we present a derivation of the Neumann-to-Dirichlet map for the case of the negative half-space, which we obtain  Damiano (2015). The dashed lines correspond to the predictions obtained using the (corrected) model presented in Milewski et al. (2015).
by a modification of the case presented in Forbes (1989). We then proceed to invert this map analytically, arriving at a singular integral representation for desired operator. We prove the convergence of the singular integral in the principal value sense. The strategy we use has the advantage of not requiring the inversion of a full matrix and also yields a fast decaying kernel. Moreover; knowing the kernel explicitly, provides the chance to carefully design approximations of the operator, e.g. by integrating the kernel exactly and interpolating only the smooth argument function. We ignore the dynamics of this air cushion between the solid and the bath and simply assume that the pressure on the solid is the one necessary to produce the consistent displacement on the fluid. This is a dynamic version of the situation presented in Keller (1998); with the peculiarity that here the contact angle is assumed to be π at all times here. We could, in principle, model the impact of a non-hydrophobic solid, provided we had a suitable model of contact angle dynamics. Given the importance of surface tension due to the size of the impacting body, in calculating the pressure, we do not assume small curvature on the impacted portion of the free surface as it is determined exactly by the curvature of the wetted surface of the object.
Impact problems belong to the class of non-smooth dynamical systems since the evolution has a switch when a particular constraint (i.e. the non-overlapping of the fluid and solid) becomes active. As we enforce this constraint exactly, it implies that implicit-in-time methods are appropriate. Once we have discretised the spatial domain there results a closed system of linear equations, provided we know exactly which points are in contact. The contact area itself is part of what needs to be determined, and this is done through an adaptive local search to find the optimal contact area satisfying the constraints and contact angle conditions. As a result, we obtain a fully predictive model for the axisymmetric impact of a smooth, convex solid onto a free surface at high Reynolds numbers. Comparisons to experimental results of hydrophobic sphere impact reported in Lee & Kim (2008) show that this relatively simple modelling can produce remarkably good results. More comprehensive experiments and comparisons are currently underway in this case.
The application of this hydrophobic impact model to the problem of bouncing droplets produces the first model that yields a prediction of the evolution of the 'contact' area and pressure field in accordance with the evolution of the free surface, free of any parameterisations. In previous models (Moláček & Bush 2013a,b;Milewski et al. 2015), a nonlinear spring was assumed to model the forces between droplet and bath, which required the fitting of three parameters. Whilst, that model also produced excellent agreement with experiments it did not enforce the impact constraint we consider here, and hence the bath and the drop did not always coincide during impact. The present work, does away with the need of a prescribed spring model for the forces and therefore simultaneously reduces assumptions and matches more accurately the experiments presented in Wind-Willansen et al. (2013), Damiano (2015) and Damiano et al. (2016). Moreover, imposing consistent displacements of the impacting bodies ensures that energy is transferred between the impacting masses with no spurious loss due to the impact modelling.
The method presented here could be modified to be applied to more complex scenarios, including waves propagating in shallower water or a varying bottom topography. For simple geometries, this can be accomplished using the method of images. More interestingly a domain decomposition method can be used to solve the case of a sharp change in depth. The study of these generalisations are part of ongoing research.
It is in principle also possible to remove the axial symmetry restrictions to model the impacts of non-axisymmetric objects or to model oblique impacts. However, for the case of walking droplets, a simpler approach would be to exploit the scale separation between the impact and the waves, and to superimpose axisymmetric impacts at different locations to produce a non-axisymmetric free surface. This has been the strategy adopted by all previous models for walkers (Bush 2015), and also could be combined with a 'skidding friction' term to model horizontal drag.
Other possible improvements would be to model the air layer dynamics and resulting friction (see Brenner 1961) and modelling contact angle dynamics for the impact of non-hydrophobic solids. Lastly one could also introduce droplet deformation, which gains in importance with larger droplets. Computational efficiency might also be improved with the use of an irregular spatial mesh, since in the drop experiment it is observed that, immediately away from the droplet, only wavelengths comparable to the Faraday wavelength need to be resolved. for some δ > 0. Given any r = (x, y), we can define p r = (x, y, 0), and the fundamental solution to the three-dimensional Laplacian centred in p r , where q = (x , y , z ). We also define the compact domain of integration Ω r shown in figure 13. Ω r is bounded by the negative half-spheres D 1 and D 2 with centre at p r and radii r 1 and r 2 , respectively; and by the plane annulus D 3 . In turn, D 3 is bounded by the circles C 1 and C 2 of radii r 1 and r 2 . We can apply Green's second identity in Ω r to functions φ and ζ r and consider the limits of r 2 → ∞ and r 1 → 0; which, in view of the decay conditions (A 1) yield The singularity in this problem is in fact integrable, meaning that the principal value in question is just a regular integral. However, this particular form is useful to perform the inversion we desire. We note that φ z is also a solution to the Laplace problem and we assume for φ z (and ∇φ z ) the same decay conditions as those given in (A 1) for φ. Thus, φ z is subject to the exact same considerations that lead to (A 3). This yields where B(r; ) = {(x , y ); r − (x , y ) < }, and since the differential relations of the fluid equations are assumed to hold even on the boundary, we have We can rewrite equation (A 5) as We can now apply Green's second identity to φ and ζ r on the plane annulus D 3 , while observing (A 7). This yields (A 8) where r 1 , r 2 , C 1 , C 2 and D 3 are as shown in figure 13, andn is the outward pointing unit normal to C 1 and C 2 . We note that where ∂n(1/|p r − q|) is positive on C 1 and negative on C 2 . We thus have where |p r − q| 2 dl(q).
(A 11a,b) We note now that J 2 and K 2 converge to zero as r 2 → ∞ on account of the decay of φ and its derivatives. In addition, we can easily see that J 1 is equals to where e(q) → 0 when r 1 → 0. The first integral on the right-hand side of (A 12) is identically zero, due to symmetry of the normal along C 1 , the second integral is bounded and thus J 1 → 0 when r 1 → 0. We now rewrite K 1 as Note now that the first integral on the right-hand side of (A 13) can be reformulated as where ξ (q − p r ) → 0 as r 1 → 0. The first integral in (A 14) is identically null by the same argument used for the first integral for (A 12), and the second integral converges to zero when r 1 → 0. We now note that the second term on the right-hand side of (A 13) can be restated as φ(p r ) and we note that given any radius r 0 > 0 arbitrarily small, for any bounded φ of class C 2 we have P.V. φ(p r ) − φ(q) Here L 2 is bounded, and we can prove that L 1 is bounded. First we note that φ(p r ) − φ(q) = −∇φ(p r )(q − p r ) − 1 2 (q − (p r )) T H(p r )(q − p r ) − ξ (q − p r )|q − p r | 2 , (A 23) where H is the Hessian matrix of φ and ξ (q − p r ) → 0 as q → p r . We then define where (r, θ ) is (q − p r ) in polar coordinates, ξ (q) → 0 as q → p r . We note that the inner integral in N 1 is identically null, and that N 2 and N 3 are bounded. Moreover, by (A 23), L 1 = N 1 + N 2 + N 3 , and thus the boundedness of L 1 follows from all the N i being bounded and the convergence of the singular integral of the N operator follows from the convergence of L 1 and L 2 . l = 1, 2, 3, . . . , n i r and m = 1, 2, 3, . . . , n θ ; where θ = 2π/n θ , and n i r is chosen to guarantee that each polar mesh covers the radial mesh completely. In particular, for all examples presented in the previous sections we use cubic interpolation. The limit in (2.23) is approximated simply by avoiding the central point of each polar mesh, thus effectively integrating everywhere outside a circle of radius r/2. We treat the value of φ estimated at each point of our polar mesh as a constant in the region limited radially by l r − r/2 and l r + r/2, and angularly by θ i − θ /2 and θ i + θ /2, see figure 14(b). We note that r is chosen to be smaller than the δr spacing of the radial mesh and that matrix N is computed once and used at every time step.
Finally, we approximate the integral of the pressure p s in the radial mesh simply using linear interpolation of the nearest mesh points.
B.1. Convergence study of the operator N We test the convergence of operator N as a function of the refinement of the radial mesh. To this end we build a radial decaying oscillatory test function that resembles the waves in our problem, namely J o (k F r)e −r 2 /(3λ 2 F ) , whose Dirichlet-to-Neumann transform we obtain using spectral methods for a fine mesh (n r = 4096). Our operator N is constructed using polar meshes with radial spacing r = δr/10 and an angular spacing that satisfies θ/2 < r/10. We compute the result of the evaluation of the N operator for different resolutions of our radial test function and we plot the L 1 norm of the error in figure 15. The least squares linear regression of the resulting logarithmic plot yields a slope m = −2.55. This suggests a convergence of at least second order. (C 9) ω 0 = 80 · 2π s −1 , (C 10) θ 0 = 0 rad. (C 11)