Appendix A. Governing equations and dimensionless groups

A solid sphere approaches with velocity $(0,0,-V)$ a rigid flat surface whose normal is $(0,0,1)$ in the axisymmetric coordinate system $(r, \theta, z)$, with a thin fluid layer in between; see figure 1(a). The surface is located at $z=0$, and its distance to the impactor (gap) can be expressed as

(A1)\begin{equation} h(r,t) = z_i(t) + \frac{r^2 }{2 R} + w(r,t) + \mathcal{O}\left[\left(\frac{r}{R}\right)^4\right] , \end{equation}

where the first term corresponds to the central line solid position of the impactor in the absence of any solid deformation, and $w(r, t)$ is the vertical deformation of the impactor surface.

The solid momentum balance in the sphere is expressed as

(A2)\begin{equation} \rho_i\,\frac{\partial \boldsymbol{v} }{\partial t} = \boldsymbol{\nabla}\boldsymbol{\cdot}\boldsymbol{\sigma} , \end{equation}

where $\boldsymbol {u}$ is the displacement of the material points, and $\boldsymbol {v} = \partial \boldsymbol {u} / \partial \boldsymbol {t}$ is their velocity, while the Cauchy stress tensor is given by

(A3)\begin{equation} \boldsymbol{\sigma} = 2 G \boldsymbol{\varepsilon} +\left(K-\tfrac{2 }{3} G\right) \mathrm{tr}(\boldsymbol{\varepsilon})\,\boldsymbol{\mathsf{I}} , \end{equation}

for a linear-elastic material in the small-strain limit, where $G$ is the shear modulus, $K$ is the bulk modulus, and $\boldsymbol {\varepsilon } = (\boldsymbol {\nabla }\boldsymbol {u} + \boldsymbol {\nabla }\boldsymbol {u}^{\rm T})/2$ is the infinitesimal-strain tensor, with $\boldsymbol {\nabla }\boldsymbol {u}$ the spatial displacement gradient. As the sphere is surrounded by a fluid, when $h$ is sufficiently small, the flow between the leading edge and the flat surface can be modelled with Reynolds lubrication equations (Szeri Reference Szeri2010)

(A4)\begin{equation} r\,\frac{\partial (\rho_m h) }{\partial t} = \frac{\partial }{\partial r} \left( \frac{r \rho_m h^3 }{\mu_m'}\,\frac{\partial p_m }{\partial r} \right) - \frac{\partial }{\partial r} \left( \frac{r \rho_m h }{2}\,v_r \right) , \end{equation}

where $p_m$ is pressure, $\rho _m$ is density, and $v_r = v_r(r,t)$ is the radial velocity of the material points at the leading edge of the sphere. The lubrication pressure is directly related to the stress in the solid impactor, $\boldsymbol {\sigma }$, via the boundary conditions of stress equilibrium at the interface (normal vector $\boldsymbol {n}$, tangential vector $\boldsymbol {t}$):

(A5a)$$\begin{gather} (\boldsymbol{\sigma} \boldsymbol{n}) \boldsymbol{\cdot} \boldsymbol{n} = p_m , \end{gather}$$

(A5b)$$\begin{gather}(\boldsymbol{\sigma} \boldsymbol{n}) \boldsymbol{\cdot}\boldsymbol{t} = \frac{h }{2}\,\frac{\partial p_m }{\partial r} + \frac{\mu_m v_r }{h} , \end{gather}$$

where the first and second terms on the right-hand side of (A5b) are, respectively, the Poiseuille and Couette flow contributions (Sprittles Reference Sprittles2024). The scaling presented in § 2 and Appendix D can be obtained alternatively by identifying the governing dimensionless groups of the physical phenomenon.

The ideal gas state equation for a polytropic process reads

(A6)\begin{equation} \frac{p }{P_0} = \left(\frac{\rho_m }{\rho_{m,0}} \right)^{\varGamma}, \end{equation}

with a certain constant $\varGamma$ that depends on the specific process. Finally, from (A1), the appearance of a dimple at location $r=r^*(t)$ is obtained by setting $\partial h(r,t) / \partial r = 0$, leading to

(A7)\begin{equation} \frac{r^*(t) }{R} + \frac{\partial w }{\partial r} = 0 , \end{equation}

which is where the geometry of the problem enters.

We can perform dimensional analysis by appropriately scaling the physical dimensions around the tip of the impactor, where most of the dynamics is supposed to take place (see figure 1b). For lubrication equations to be valid, we require $H/\mathcal {L} \ll 1$, i.e. the height of the gap to be much smaller than its width. This results in the radial velocity in the fluid being much greater than the vertical one, to satisfy mass conservation. We also assume that the gap height is set by the scale of solid vertical deformation $w(r,t)$. The resulting scaling is

(A8a)$$\begin{gather} t =\tau \tilde{t}, \end{gather}$$

(A8b)$$\begin{gather}\boldsymbol{v} =V \tilde{\boldsymbol{v}} , \end{gather}$$

(A8c)$$\begin{gather}\left[r, r^*(t) \right] = \mathcal{L} \left[\tilde{r},\tilde{r}^*(\tilde{t})\right] , \end{gather}$$

(A8d)$$\begin{gather}\left[ h, w \right]= H \left[\tilde{h},\tilde{w} \right] , \end{gather}$$

(A8e)$$\begin{gather}p(r,t) = P\,\tilde{p}(\tilde{r}, \tilde{t}) , \end{gather}$$

(A8f)$$\begin{gather}\boldsymbol{\sigma}(r,w,t)=P\,\tilde{\boldsymbol{\sigma}}(\tilde{r}, \tilde{w}, \tilde{t}) , \end{gather}$$

where we set the same pressure for the fluid and the solid stress. To choose the appropriate scale for $v_r$, we make a distinction between the elastic and inertial regimes. In the quasi-static limit, Hertzian analysis tells us that the maximum radial displacement is at the contact radius (Johnson Reference Johnson1985). Assuming that the contact radius coincides with the lateral scale of the dimple $\mathcal {L}$, and the pressure is that of the elastic regime, $P_{el}$, we find $U_r \sim {P_{el} \mathcal {L} }/{E} \sim ( {G }/{E } )H \sim H$, where we assume the Young's and shear moduli to be of the same order. Taking into account the time scale of the problem $\tau$, it follows that $v_r^{el} \sim H / \tau _{el}$. Conversely, in the inertial regime, we estimate in Appendix D that $\varepsilon _{zz}\sim V / c_p$. In the nearly incompressible limit, we have $\varepsilon _{rr} \sim \varepsilon _{zz}$, resulting in a velocity $v_r^{in} \sim \mathcal {L} V/ \tau _{in} c_p$. We proceed by plugging in the rescaled quantities in the equations discussed formerly. Equation (A2) yields

(A9a,b)\begin{equation} \mathcal{G}_i\,\frac{\partial \tilde{\boldsymbol{v}} }{\partial \tilde{t}} = \tilde{\boldsymbol{\nabla}} \boldsymbol{\cdot} \tilde{\boldsymbol{\sigma}} , \quad \mathcal{G}_i = \frac{\rho_i V\mathcal{L}}{\tau P}, \end{equation}

which is the ratio of the solid inertia pressure scale over the pressure scale of the problem. From (A4), plugging in the scale $v_r^{el}$ gives

(A10a,b)\begin{align} \mathcal{G}_m\,\frac{\partial (\tilde{p}^{1/\varGamma} \tilde{h}) }{\partial \tilde{t}} - \frac{1 }{\tilde{r}} \left( \tilde{r} \tilde{p}^{1/\varGamma} \tilde{h}^3\,\frac{\partial \tilde{p} }{\partial \tilde{r}} \right) ={-} \mathcal{G}_m \delta_{el}\,\frac{1 }{\tilde{r}} \left( \tilde{r} \tilde{p}^{1/\varGamma} \tilde{h}\,\frac{\tilde{V}_r }{2} \right) , \quad \mathcal{G}_m = \frac{\mu_m' \mathcal{L}^2 }{P \tau H^2} , \end{align}

where $\delta _{el} = (\mu '_m V / G R)^{1/5} \sim H_{el} / \mathcal {L}_{el} \ll 1$, and $\mathcal {G}_m$ is the lubrication dimensionless group. We conclude that the Couette term is an order $\delta _{el}$ smaller than the others, and can be neglected. Following the same logic, using $v_r^{in}$ in place of $v_r^{el}$, we find that the Couette contribution has a prefactor $\mathcal {G}_m V/c_p$, and it is negligible insofar as $V/ c_p \ll 1$.

The dimple geometry (A7) leads to

(A11a,b)\begin{equation} \mathcal{G}_d \tilde{r}^* + \frac{\partial \tilde{w}}{\partial \tilde{r}} = 0 , \quad \mathcal{G}_d = \frac{\mathcal{L}^2 }{R H} . \end{equation}

In (A3), there are several possibilities to gauge the value of strain. If we assume the solid to deform over the scale of the dimple, then strain can be scaled as $\boldsymbol {\varepsilon }=(H/\mathcal {L})\tilde {\boldsymbol {\varepsilon }}$, which leads to

(A12a,b)\begin{equation} \tilde{\boldsymbol{\sigma}} = \mathcal{G}_e \left( 2 \tilde{\boldsymbol{\varepsilon}} + \left(\frac{K}{G} - \frac{2 }{3}\right) {\rm tr}(\tilde{\boldsymbol{\varepsilon}})\,\boldsymbol{\mathsf{I}} \right) , \quad \mathcal{G}_e = \frac{G H }{P \mathcal{L}} . \end{equation}

Alternatively, for stresses manifesting where pressure waves propagate, we can write $\boldsymbol {\varepsilon }=({H }/{c_p \tau } )\tilde {\boldsymbol {\varepsilon }}$ and

(A13a,b)\begin{equation} \tilde{\boldsymbol{\sigma}}\approx \mathcal{G}_M \tilde{\boldsymbol{\varepsilon}} , \quad \mathcal{G}_M = \frac{H M}{c_p \tau P } , \end{equation}

where we assumed a uniaxial strain scenario, and $M= 2 G (1-\nu )/(1-2 \nu )$ is the pressure wave (P-wave) modulus of the solid, and $\nu$ is its Poisson's ratio. Additionally, gas compressibility effects (A6) are described by the dimensionless group $\mathcal {G}_c = (P/P_0)^{1/\varGamma }$.

Assuming dimple formation and thin-film flow in effect corresponds to imposing $\mathcal {G}_d =1 \implies \mathcal {L}= \sqrt {R H}$, and $\mathcal {G}_m =1$. The elastic, inertial and solid-compressible regimes are obtained by enforcing $\mathcal {G}_e$, $\mathcal {G}_i$ and $\mathcal {G}_M$ being of order one, respectively. By doing so, we obtain the values of $\mathcal {L}, H, \tau, P$ in each regime, and we can substitute them in the dimensionless groups to assess their relative importance.

An analogous way of reasoning can be used to scale the droplet impact. The full incompressible Navier–Stokes equation is

(A14)\begin{equation} \rho_i\,\frac{\partial \boldsymbol{v} }{\partial t} = \boldsymbol{\nabla}\boldsymbol{\cdot}\boldsymbol{\sigma} -\rho_i (\boldsymbol{v} \boldsymbol{\cdot}\boldsymbol{\nabla} \boldsymbol{v}) , \end{equation}

when we used the Newtonian stress for the fluid

(A15)\begin{equation} \boldsymbol{\sigma} = p_i \boldsymbol{\mathsf{I}} + \mu_i (\boldsymbol{\nabla} \boldsymbol{v} + \boldsymbol{\nabla} \boldsymbol{v}^\mathrm{T} ) . \end{equation}

The normal boundary condition at the droplet–gas interface is

(A16)\begin{equation} (p-\boldsymbol{\sigma} \boldsymbol{\cdot}\boldsymbol{n})= \gamma \kappa , \end{equation}

where the mean curvature $\kappa$ characteristic value can be approximated as $H/\mathcal {L}^2$; see § 4.3. The gas state equation and the dimple condition remain unchanged.

By performing the same scaling, we obtain from (A14) that

(A17a,b)\begin{equation} \mathcal{G}_i\,\frac{\partial \tilde{\boldsymbol{v}} }{\partial \tilde{t}} = \tilde{\boldsymbol{\nabla}} \boldsymbol{\cdot}\tilde{\boldsymbol{\sigma}} - \frac{\rho_i V^2 }{P} \, (\tilde{\boldsymbol{v}} \boldsymbol{\cdot} \tilde{\boldsymbol{\nabla}} \tilde{\boldsymbol{v}}), \quad \mathcal{G}_i = \frac{\rho_i V\mathcal{L}}{\tau P} \end{equation}

and $(\rho _i V^2/P)/\mathcal {G}_i = H/\mathcal {L} \ll 1$, which entails that the nonlinear term can be neglected. Equation (A16) yields

(A18a–c)\begin{equation} \tilde{p} - \frac{p_i }{P}\,(\boldsymbol{\mathsf{I}} \boldsymbol{\cdot} \boldsymbol{n}) + \mathcal{G}_v ((\tilde{\boldsymbol{\nabla}} \tilde{\boldsymbol{v}} + \tilde{\boldsymbol{\nabla}} \tilde{\boldsymbol{v}}^\mathrm{T} ) \boldsymbol{\cdot} \boldsymbol{n}) = \mathcal{G}_s\,\frac{\tilde{h} }{\tilde{r}^{*2} },\quad \mathcal{G}_v = \frac{\mu_i V }{\mathcal{L} P}, \quad \mathcal{G}_s = \frac{\gamma H }{\mathcal{L}^2 P}, \end{equation}

where we defined the inertial and the capillary dimensionless groups. Assuming $p_i/P \ll 1$, the viscous, inertial and capillary regimes are obtained by enforcing $\mathcal {G}_v$, $\mathcal {G}_i$ and $\mathcal {G}_s$ being of order one, respectively.

Appendix B. Auxiliary parameters to assess fluid compressibility

The effects of fluid compressibility arise when the pressure becomes of the order of the atmospheric one. The critical height (characteristic value) when the compressibility effects become relevant is $H^* = (R \mu _m' V /P_0)^{1/2}$, where $P_0$ is the ambient pressure; unlike Mandre et al. (Reference Mandre, Mani and Brenner2009), we used $\mu _m'$ in place of $\mu _m$. By comparing the dominant pressure scale to the ambient one, we obtain so-called ‘compressibility parameters’ for each regime:

(B1a)$$\begin{gather} \epsilon_{el}=\frac{P_0 }{P_{el}} \approx \frac{P_0 }{(R^{{-}1} \mu_m' G^4 V )^{1/5}} , \end{gather}$$

(B1b)$$\begin{gather}\epsilon_{in}= \frac{P_0 }{P_{in}} \approx \frac{P_0 }{(R \mu_m'^{{-}1} V^7 \rho_i^4)^{1/3}}, \end{gather}$$

where (B1b) is the well-known compressibility parameter for high-speed droplets (Mandre et al. Reference Mandre, Mani and Brenner2009), and (B1a) is relevant in the elasticity-dominated scenario.

Appendix C. Simulation methodology

Our approach is conceptually similar to that used by Chubynsky et al. (Reference Chubynsky, Belousov, Lockerby and Sprittles2020) for isothermal impacts of droplets, by Chakraborty, Chubynsky & Sprittles (Reference Chakraborty, Chubynsky and Sprittles2022) for soft solids Leidenfrost levitation, and by Habchi et al. (Reference Habchi, Eyheramendy, Vergne and Morales-Espejel2008) for hard solids: the complete impactor is represented in the model, unlike some studies that consider its leading edge as a half-space (e.g. Davis et al. Reference Davis, Serayssol and Hinch1986; Hicks & Purvis Reference Hicks and Purvis2010; Mani et al. Reference Mani, Mandre and Brenner2010). We recall the methodology details in this appendix. We model the interaction between the impacting solid particle and the fluid surrounding it by means of monolithic finite-element simulations for the fluid–structure interaction problem. To simplify the computation, we study the axisymmetric case as we are interested primarily in normal impacts. The schematic of the simulation is illustrated in figure 1(b), and we define ‘$\partial \varOmega ^-$’ as the boundary of the impactor leading edge (the bottom half of the sphere in the figure).

For the solid mechanics, we employ the neo-Hookean model, which has been widely used to model polymers and rubber-like materials undergoing large deformations. In this setting, $\boldsymbol {x} = (r, \theta, z)$ are the spatial coordinates, $\boldsymbol {u} = (u_r(r,z),0, u_z(r,z))$ is the displacement field, and the deformation gradient is defined as $\boldsymbol{\mathsf{F}} = \boldsymbol{\mathsf{I}} + \boldsymbol {\nabla } \boldsymbol {u}$.

Consequently, denoting by $(\cdot )^{\rm T}$ the transpose operator, the Green–Lagrange tensor is written as in Marsden & Hughes (Reference Marsden and Hughes1994):

(C1)\begin{equation} \boldsymbol{\mathsf{E}} = \tfrac{1 }{2} \left( \boldsymbol{\mathsf{F}}^{\rm T} \boldsymbol{\mathsf{F}} - \boldsymbol{\mathsf{I}} \right) . \end{equation}

The corresponding work-conjugate is the second Piola–Kirchhoff stress defined as

(C2)\begin{equation} \boldsymbol{\mathsf{S}} = J \boldsymbol{\mathsf{F}}^{{-}1} \boldsymbol{\sigma} \boldsymbol{\mathsf{F}}^{-{\rm T}} , \end{equation}

where $J = \det (\boldsymbol{\mathsf{F}})$, and $\boldsymbol {\sigma }$ is the Cauchy stress tensor. The neo-Hookean strain energy density is defined by

(C3)\begin{equation} \psi(\boldsymbol{\mathsf{E}}) = G \left(\mathrm{tr}(\boldsymbol{\mathsf{E}}) - \ln J \right) + \frac{\lambda }{2}\, (\ln J )^2 , \end{equation}

where $\lambda, G$ are the Lamé constants of the impactor material. These are related to the Young's modulus $E$ and the Poisson's ratio $\nu$, via the usual relations:

(C4a,b)\begin{equation} E = \frac{ G (3\lambda + 2 G) }{\lambda + G} , \quad \nu = \frac{\lambda }{2 (\lambda + G)} . \end{equation}

The neo-Hookean hyperelastic constitutive relationship is obtained by deriving the strain energy density with respect to the Green–Lagrange strain:

(C5)\begin{equation} \boldsymbol{\mathsf{S}} = \frac{\partial \psi }{\partial \boldsymbol{\mathsf{E}}} . \end{equation}

In the absence of volume forces, the linear momentum balance equation in the reference configuration for a solid undergoing finite deformations is

(C6)\begin{equation} \rho_0\,\frac{\partial^2 \boldsymbol{u} }{\partial t^2} = \boldsymbol{\nabla} \boldsymbol{\cdot} (\boldsymbol{\mathsf{F}}\boldsymbol{\mathsf{S}})^\mathrm{T} , \end{equation}

where $\rho _0$ is the density in the initial configuration.

For the cushioning fluid, we use the lubrication approximation in incompressible axisymmetric form (Szeri Reference Szeri2010), (A4), where $\rho _m$ can be simplified. The solid impactor is assigned a homogeneous velocity field $(0,0,-V)$, while on $\partial \varOmega ^-$ we enforce the stress continuity (A5). The rest of the boundary is stress-free. Finally, the boundary condition $p_m(r=R) = 0$ is enforced, effectively neglecting the contribution of ambient pressure.

All the simulations are performed using the commercial software COMSOL Multiphysics 6.0 and the built-in modules ‘Solid mechanics’ and ‘Thin-film flow, edge’ for the solid and fluid mechanics, respectively. An implicit backward differentiation scheme is used, with second-order accuracy. Quadratic elements are used for the structural mechanics, while the lubrication equation is solved on the one-dimensional boundary $\partial \varOmega ^-$. Quadratic interpolation functions are used for both the fluid flow and the structural mechanics. A damped Newton method algorithm is adopted, with minimum relative tolerance $10^{-5}$ on the weighted residual. The Parallel Sparse Direct Solver (PARDISO) is used to solve the resulting linear algebraic system of equations.

The initial gap, $h_0 \approx 12\ \mathrm {\mu } \mathrm {m}$, is set such that no considerable deformation takes place in the early stages of the simulation (see the elasticity parameter (Davis et al. Reference Davis, Serayssol and Hinch1986) $\epsilon = 4 \mu _m V R^{3/2}(1-\nu _i^2)/E h_0^{5/2}\approx 0.004\unicode{x2013}0.35$, where the highest value is obtained in the case $V=5 \ \mathrm {m}\ \mathrm {s}^{-1}$ for the soft impactor). A consistent initialization step is performed at the beginning of the simulation, and we verify that the initial pressure is much smaller than the final one. At the same time, in order to neglect inertial effects in the thin film, we require $\rho _i V h_0/ \mu _m \ll 1$ (Davis et al. Reference Davis, Serayssol and Hinch1986), therefore we should not increase $h_0$ significantly. Simulations with the same time step showed limited discrepancy with values of the initial gap approximately $6\ \mathrm {\mu } \mathrm {m}$ and $24\ \mathrm {\mu } \mathrm {m}$ (see figure 3), except near the transition, where the small numerical discrepancies determine whether the solid glides on the fluid film or locally pierces it. The influence of the parameter $h_0$ on the simulations is due to the fact that we need to account explicitly for wave propagation in the solid, unlike, for example, Davis et al. (Reference Davis, Serayssol and Hinch1986) and Bertin (Reference Bertin2024). In our dynamic scenario, it would take a wave ($c_p\sim 20 \ \mathrm {m}\ \mathrm {s}^{-1}$) approximately ${\sim }10^{-3}\ \mathrm {s}$ to propagate through the impactor. The impactor ($V \sim 1\ \mathrm {m}\ \mathrm {s}^{-1}$) would have to move a distance $h_{fall} \sim 1 \mathrm {mm}$ for a steady-state solution to be established. At that distance from the wall, the lubrication approximation is inadequate ($\rho _m V h_{fall} /\mu _m\sim 100$) and would not reflect the actual physics, therefore making the quasi-static approximation incorrect for our case. Performing simulations with the full Navier–Stokes equations for the fluid, along with elastodynamics in the solid, should yield results independent of $h_0$ for large enough values, since lubrication stresses would become dominant at the correct distance from the wall, and the elastodynamics would evolve naturally. Note that also for droplets, the equations being solved in the impactor are elliptic in nature (Hicks & Purvis Reference Hicks and Purvis2010; Mani et al. Reference Mani, Mandre and Brenner2010), unlike our hyperbolic case. That is why droplet simulations can be initialized from any (high enough) height: the boundary of the droplet adapts locally to the evolving boundary conditions, and no wave propagation has to be resolved in the bulk. In conclusion, the chosen value of $h_0$ is consistent with the thin-film assumption and the elastodynamic response of the impactor, and yields results for the lateral size of the trapped air pocket in good agreement with the experiments.

Mesh sizes ranging from $5\ \mathrm {\mu } \mathrm {m}$ to $40\ \mathrm {\mu } \mathrm {m}$ on the region where contact is expected to take place give similar results when $\phi \gg 1$ and $\phi \ll 1$. The same occurs for the time step size: defining $T=h_0/V$ as the time it would take for the impactor to make contact if no lubricant was present, time steps of the order $T/10$ and $T/100$ yield similar results in the slow and fast impact limits. However, in the case $\phi \gtrsim 1$, a highly loaded thin film of a few hundred nanometres leads to gliding (see figure 4b). In this case, we observe earlier a numerical instability similar to the one reported also by Stupkiewicz (Reference Stupkiewicz2009) in steady-state simulations of soft lubrication. To the best of the authors’ knowledge, the origin of this instability for soft solids is still unclear, and to investigate it thoroughly may require the addition of extra physics in the thin film, similarly to the work of Chubynsky et al. (Reference Chubynsky, Belousov, Lockerby and Sprittles2020), which included van der Waals forces for droplets. The final choice is to adopt a mesh size $20\ \mathrm {\mu } \mathrm {m}$ and a time step $\Delta t = T/30$, which reduces the error tracked in the simulation, renders smooth stable results when the nanometric film height makes the conditioning of the finite element matrix poor, and furthermore gives results consistent with experiments; recall figure 3. The enhanced stability is due to the numerical damping of the implicit scheme, which decreases with finer time steps. Ultimately, we trust the simulation outcomes because they are able to reproduce the experimental measurements in all cases, even when the gliding phenomenon takes place.

The simulations are terminated when the impactor first makes contact with the flat surface, i.e. $h_{min} \approx 0$. Reaching this value is a numerical artefact, since refining mesh and time step allows us to get ever closer to touchdown (except for the simulations displaying gliding, as mentioned in the previous paragraph), but this leads eventually in all cases to a singularity of stress when the gap is almost closed, featuring the aforementioned numerical oscillations. Nevertheless, the data reported are within $8\,\%$ error from data recorded at the first instance any point of the impactor reaches the $100$ nm threshold, which is approximately where kinetic effects would become relevant in the air (Batchelor Reference Batchelor2000; Mani et al. Reference Mani, Mandre and Brenner2010). This choice allows us to compare our data directly to those of Zheng et al. (Reference Zheng, Dillavou and Kolinski2021) at touchdown, whose experimental error is similar to what we achieved with simulations. For values $\phi <0.1$, the solver is unable to converge as a singularity of stress is reached at approximately $400 \ \mathrm {nm}$, so we do not report those values, even if the observed pressure and gap profiles are consistent with those in the elastic regime; see figure 4(a).

Appendix D. Solid compressibility effects

Figure 6(a) shows how releasing the nearly incompressible assumption for the solid completely changes the strain distribution within it. It is speculated that the lubrication pressure can be balanced by volumetric stresses in the case of compressible solids. We estimate the value of these stresses ($P_M$) by considering a uniaxial strain scenario so that $P_M \sim M \varepsilon _{zz}$, where $M= 2 G (1-\nu )/(1-2 \nu )$ is the P-wave modulus of the solid, $\nu$ is its Poisson's ratio, and $\varepsilon _{zz}$ is the uniaxial strain. The vertical displacement at the tip is of the order ${\sim }H$, while the part of the impactor that undergoes deformation is within a radius ${\sim }c_p \tau$ from the tip, where $c_p = \sqrt {M/\rho _i}$ is the P-wave speed within the solid. Therefore, we can estimate $P_M \sim \rho _i c_p V$. Assuming that this pressure balances $P_{lub}$, we obtain the following scales associated with impactor deformation:

(D1a,b)\begin{equation} \mathcal{L}_{M} = \left( \frac{R^3 \mu_m' }{\rho_i c_p} \right)^{1/4} , \quad H_{M} = \left( \frac{R \mu_m' }{\rho_i c_p} \right)^{1/2} , \end{equation}

which do not depend on the impact velocity. Unlike the scales shown in § 2, these scales have not been reported previously. We can also analyse the conditions for the inertial regime to feature meaningful volumetric deformation by setting $P_M \sim P_{in}$. Pursuing reasoning similar to that which led to the definition of $\phi$ (§ 2.1), we evaluate the ratio between $P_M$ and $P_{in}$ both when inertial scales (2.1a,b) are dominant and when those in (D1a,b) are. Thus we obtain a second transition parameter by enforcing them to be of the same magnitude at the said transition:

(D2)\begin{equation} \left( \frac{P_{in}[\mathcal{L}_{in}, H_{in}] }{P_{M}[\mathcal{L}_{in}, H_{in}]} \right)^{3/4} = \frac{P_{in}[\mathcal{L}_{M},H_{M}]}{P_{M}[\mathcal{L}_{M}, H_{M}] } = \psi = \frac{V\,St^{1/3}}{c_p } = \frac{\mathcal{L}_{in} / c_p }{H_{in} / V} = \frac{\tau'_{propagation} }{\tau_{impact}} , \end{equation}

which again can be interpreted as the ratio of the time it takes compressive waves to propagate in the solid ($\tau '_{propagation}$) over the lubrication time scale in the thin film. In addition, $\psi$ is a function of $\phi$:

(D3)\begin{equation} \psi = \frac{c_s }{c_p}\,\phi = \sqrt\frac{(1 - 2 \nu) }{2 (1-\nu)}\,\phi . \end{equation}

The fact that $\psi < \phi$ implies that solid compressibility can play a role only in the inertial regime. Moreover, the transition (the point at which $\psi \sim 1$) depends on the Poisson's ratio value: the closer to incompressibility ($\nu = 1/2$), the greater the value of $\phi$ (e.g. higher impact velocities) for which the solid compressibility would become relevant.

We performed an additional set of simulations for Poisson's ratio spanning the interval $[0.1,0.4]$, and velocities in $0.1\unicode{x2013}10\ \mathrm {m}\ \mathrm {s}^{-1}$, while keeping the other parameters fixed. The simulations’ outcome of the dimple radius as a function of $\psi$, as defined in (D2), can be examined in figure 9.

Figure 9.

Numerically computed dimple radius as a function of the transition parameter $\psi$ for different values of Poisson's ratio. For almost incompressible impactors ($\nu = 0.47$), the sharp drop is evident, but it fades away as the solid becomes more easily compressible; by $\nu = 0.2$, the drop disappears. The later scale obtained by balancing the lubrication pressure in the thin film with the volumetric stress in the solid proves correct for $\psi >3$.

Remarkably, for $\psi > 3$, $r_{contact} \approx 2.9 \, \mathcal {L}_M$ independently of compressibility. For low values of the Poisson's ratio, the bubble radius evolves smoothly and there is no sharp drop in the transition between regimes. We have posited that the absence of pressure spikes, which would be released immediately due to volumetric compressibility, prevents the impactor from gliding (figure 10a). Looking at the pressure and free-surface profiles for $\psi = 20$, where compressibility in the solid is fully dominant (see figure 10b), the free-surface profiles are qualitatively similar to those observed for the inertial regime (figures 4c); however, the pressure profile appears much flatter in the late stages of impact. The lengths, times and pressure scales are consistent with the scaling that we proposed in this appendix.

Figure 10.

Evolution of the impactor free-surface (top plots) and the pressure profile (bottom plots) during early stages of impact for values $\nu =0.1$ ($\tau _{in} = H_{in}/V$ and $\tau _{M} = H_{M}/V$). In the top plots, we use the same value of $\phi$ as in figure 4(b), which corresponds to a higher value of $\psi$. Compressibility prevents the formation of the pressure spikes that would cause gliding. The time of impact is therefore shorter. In the bottom plots, the profiles are in the solid compressible regime. Compared to the inertial case, there is flattening of the pressure in the final stages. Here, (a) $\phi =4.0$, $\psi =2.7$, $V = 1.1 \ \mathrm {m}\ \mathrm {s}^{-1}$, and (b) $\phi =30$, $\psi =20$, $V = 5.2 \ \mathrm {m}\ \mathrm {s}^{-1}$.

Finally, with regard to the time span as a function of compressibility, the case $\nu = 0.1$ lasts until $t/\tau _{in} \approx 5$ (figure 10a), while $\nu = 0.47$ reaches $t/\tau _{in} \approx 10$ (figure 4b). The time scale $\tau _{in} = H_{in}/V$ is the same in both cases, as $H_{in}$ in (2.1a,b) does not depend on $\nu$, which is the only value that changes between these simulations, thus a direct comparison of times in absolute terms is warranted. Hence these results showcase how near incompressibility alone can significantly extend the duration of the approach phase, leading to larger entrapped air bubbles. Solid compressibility adds another layer of complexity compared to fluid droplets, and another dimensionless parameter to consider, $\psi$.

Appendix E. Verification of characteristic values from simulation data

See figure 11.

Figure 11.

(a,c,e) Dimple length and (b,d,f) height for the (a,b) elastic, (c,d) inertial, and (e,f) solid compressibility regimes, made dimensionless with the respective parameters.