1. Introduction
Collisions of liquid drops are key to a wide variety of processes, including cloud formation (Grabowski & Wang Reference Grabowski and Wang2013), combustion engines (Zhang et al. Reference Zhang, Zhang, Yuan and Zhang2016), disease transmission (Gralton et al. Reference Gralton, Tovey, McLaws and Rawlinson2011) and spray drying of suspensions such as milk (Finotello et al. Reference Finotello, Kooiman, Padding, Buist, Jongsma, Innings and Kuipers2017). During the initial stage of the collision, a layer of gas is trapped between the drops, and consequently the drops can either bounce off each other, permanently coalesce/merge or coalesce then separate. Central to the understanding of these various processes is the question of whether a pair of colliding drops will merge or bounce.
Drop collisions have been studied as far back as the 19th century, with Rayleigh (Reference Rayleigh1899) discovering that drops in a jet can bounce off each other. Early experiments focused on water drops (e.g. Brazier-Smith, Jennings & Latham Reference Brazier-Smith, Jennings and Latham1972; Ashgriz & Poo Reference Ashgriz and Poo1990), suggesting that head-on collisions always led to merging, and focused on whether drops that merged would permanently merge or separate. More recently, experiments with hydrocarbons, starting with Jiang, Umemura & Law (Reference Jiang, Umemura and Law1992), identified that drops colliding head-on could actually bounce. They found two transitions in collision behaviour with increasing Weber number, 
$\textit{We}=\rho V^2R/\gamma$, for drop radius 
$R$, relative impact speed 
$V$, liquid density 
$\rho$ and surface tension 
$\gamma$. These are a low-speed transition from merging to bouncing (termed the ‘soft transition’, and can be thought of as less ‘violent’, observed at critical Weber number 
$\textit {We}_S$), and a higher-speed transition from bouncing back to merging (the ‘hard transition’, i.e. more ‘violent’, at 
$\textit {We}_H$). Qian & Law (Reference Qian and Law1997) then provided an extensive study of the process, including identifying the effects of changing the ambient gas and its pressure, to isolate the influence of the gas film. They discovered strong dependencies of the critical 
$\textit {We}$ numbers on the gas properties. Notably, they discovered that water drops colliding head-on could bounce off each other if the ambient gas pressure was raised sufficiently high. Conversely, Willis & Orme (Reference Willis and Orme2000) found that in a vacuum the bouncing regime disappears and the drops always merge.
Early work focused on characterising collision events based purely on 
$\textit {We}$, with different 
$\textit {We}_S$ and 
$\textit {We}_H$ for different liquid and gas properties. The aforementioned analyses and experiments controlled 
$\textit {We}$ by varying both the relative impact speed and the drop radius. However, as identified in simulations by Li (Reference Li2016), for a fixed value of 
$\textit {We}$ drops can be made to either bounce or merge by altering the drop radius. This was confirmed experimentally by Huang & Pan (Reference Huang and Pan2021), where experiments were performed to determine 
$\textit {We}_S$ and 
$\textit {We}_H$ for a range of different drop radii. They found that water drops colliding head-on could bounce if their radius was large enough, and also discovered for all liquids considered that at sufficiently small radius drops always merge.
A process related to head-on drop collisions is the impact of a liquid drop onto a smooth solid surface (here, drop–drop events are ‘collisions’ and drop–solid ones are ‘impacts’). The study of drop-on-solid impacts has the advantage over drop–drop collisions that high-speed imaging techniques (Josserand & Thoroddsen Reference Josserand and Thoroddsen2016) can measure the thickness of the air film from beneath the solid (Driscoll & Nagel Reference Driscoll and Nagel2011; Kolinski et al. Reference Kolinski, Rubinstein, Mandre, Brenner, Weitz and Mahadevan2012; van der Veen et al. Reference van der Veen, Tran, Lohse and Sun2012; de Ruiter, Mugele & van den Ende Reference de Ruiter, Mugele and van den Ende2015b; Langley, Li & Thoroddsen Reference Langley, Li and Thoroddsen2017; Lo, Liu & Xu Reference Lo, Liu and Xu2017; Pack et al. Reference Pack, Hu, Kim, Zheng, Stone and Sun2017), with thicknesses as small as tens of nanometres observed. Relatively recently, in Kolinski, Mahadevan & Rubinstein (Reference Kolinski, Mahadevan and Rubinstein2014) and de Ruiter et al. (Reference de Ruiter, Lagraauw, van den Ende and Mugele2015a), it was discovered that the trapped gas allows drops to bounce without contact (in contrast to contact-driven bouncing from super hydrophobic surfaces) on sufficiently smooth surfaces, with similar transitions between bouncing and contact to those seen in drop–drop collisions identified. Due to the difficulty of producing surfaces with sub-nanometric roughness, experiments have also considered coating surfaces with a very thin liquid film (Lo et al. Reference Lo, Liu and Xu2017; Zhang et al. Reference Zhang, Soori, Rokoni, Kaminski and Sun2021a), in order to mimic impact on a smooth solid. A review of these works can be found in Sprittles (Reference Sprittles2024).
In impact experiments where the gas film's dynamics can be measured, different types of contact have been identified. For example, Zhang et al. (Reference Zhang, Soori, Rokoni, Kaminski and Sun2021a) identified that contact can occur either at a ‘kink’ at the edge of the gas film, at a point mid-way along the film, or at a ‘dimple’ at the centre of the film. The mode of contact is related to the time at which merging occurs and can be manipulated in a non-trivial way by changing material parameters and impact speed. In contrast, much less is known about contact modes for drop–drop collisions, as the trapped gas films cannot, at present, be probed experimentally.
It is now well established that, to capture transitions between bouncing and merging in drop–solid impacts, unconventional models for the gas’ dynamics are necessary (Sprittles Reference Sprittles2024). This is because nanoscale gas films are not accurately described by the Navier–Stokes equations with no-slip boundary conditions, as the length scales involved are comparable to, or even less than, the mean free path 
$\lambda$ in the gas. This is characterised by the Knudsen number, 
$\textit {Kn}=\lambda /h$, where 
$h$ is a characteristic length scale; here, the height of the film. Notably, the classical Navier–Stokes with no-slip boundary conditions are accurate for 
$\textit {Kn}\ll 1$.
Gas-kinetic corrections that depend on the mean free path were first introduced to model the approach of rigid spheres using a slip boundary condition, valid for small Knudsen number (Davis Reference Davis1972; Hocking Reference Hocking1973). As we shall see, they can crudely be thought of as modifying the ‘effective viscosity’ of the film and act to enhance the flow rate, for a given pressure gradient, as 
$\textit {Kn}$ becomes larger. Bach, Koch & Gopinath (Reference Bach, Koch and Gopinath2004) used solutions of the Boltzmann equation for the gas film between weakly deformable drops, building on their earlier work using potential flow to solve for the flow within the drops (Gopinath & Koch Reference Gopinath and Koch2002). This involved a transfer of knowledge from well-known applications of gas lubrication in narrow regions with fixed shape to free-surface flow dynamics (Sundararajakumar & Koch Reference Sundararajakumar and Koch1996). The effect of 
$\textit {Kn}$ on splashing and bubble entrapment in drop–solid impacts was considered by Duchemin & Josserand (Reference Duchemin and Josserand2012). The corrections used in our work are those derived in Chubynsky et al. (Reference Chubynsky, Belousov, Lockerby and Sprittles2020), which are based on asymptotic results for the large and small Knudsen number limits with a fit to solutions of the Boltzmann equation for intermediate values.
Any model using the incompressible Navier–Stokes equations for the gas film cannot capture experimentally observed dependencies on gas ambient pressure, as only derivatives of gas pressure appear in the equations. At higher-speed impacts, it could be that compressibility is important (Mandre, Mani & Brenner Reference Mandre, Mani and Brenner2009; Hicks & Purvis Reference Hicks and Purvis2011; Josserand & Thoroddsen Reference Josserand and Thoroddsen2016), but at the speeds we are interested in, where transitions between bouncing and merging occur, the pressure deviations observed are small enough that incompressibility is a reasonable assumption (Sprittles Reference Sprittles2024). Our model will utilise gas-kinetic corrections derived from the Boltzmann equation that are dependent on the mean free path, and hence on the ambient pressure, even when an incompressible flow model is used. In particular, the mean free path of the gas varies inversely proportionally to the ambient pressure, thus providing a connection of the model with experimental observations.
In addition to gas-kinetic corrections, the effect of the van der Waals force between the drops must also be incorporated into the model (Pan et al. Reference Pan, Law and Zhou2008). This takes the form of a disjoining pressure applied to the surface of the drops, and depends on the distance between the drops (i.e. the gas film height). Without this effect, which acts to pull interfaces together which are 
${\sim }$tens of nanometres apart, the drops will always bounce (Sprittles Reference Sprittles2024).
The drop-collision phenomenon is typically a high-deformation flow in which computational methods become crucial, particularly to predict the hard transition. As reviewed in Sprittles (Reference Sprittles2024), numerical methods have typically failed to reliably reproduce transitions between bouncing and merging here, due to a failure to incorporate the required physics – errors are further compounded by inabilities to resolve the strongly multiscale flow (the drop is often millimetre sized, whilst the gas film can be nanometres thick).
Our model builds on the work of Li (Reference Li2016) for drop–drop collisions, who identified the need to include additional nanoscale physics in the form of gas-kinetic effects and the van der Waals driven disjoining pressure, but considered only a couple of isolated experimental cases, and Chubynsky et al. (Reference Chubynsky, Belousov, Lockerby and Sprittles2020), who implemented these additional effects in a lubrication framework, but focused on drop-on-solid impacts. Chubynsky et al. (Reference Chubynsky, Belousov, Lockerby and Sprittles2020) used drop–drop collisions as a benchmark against the simulations of Li (Reference Li2016), and briefly considered the contact modes, but did not investigate the bounce–merge transitions. Both used a finite element approach combined with an arbitrary-Lagrangian–Eulerian mesh design (Anthony et al. Reference Anthony2023) to accurately capture the interfacial dynamics. We have incorporated the (previously closed-source) model of Chubynsky et al. (Reference Chubynsky, Belousov, Lockerby and Sprittles2020) into a robust open-source framework, using the finite element package oomph-lib (Heil, Hazel & Matharu Reference Heil, Hazel and Matharu2022), and for the first time applied the lubrication methodology to investigate the transitions between bouncing and merging in drop–drop collisions, quantitatively comparing our computed transitions in drop-collision behaviour with a variety of experimental results.
In contrast, volume-of-fluid methods using an Eulerian mesh have the advantage that they can easily capture changes in topology, but are at a disadvantage because they require a very fine mesh to resolve the dynamics and position of the gas layer. In recent years, major progress has been made towards achieving this aim (Liu & Bothe Reference Liu and Bothe2019; Sharma & Dixit Reference Sharma and Dixit2021), but the incorporation of necessary additional physics remains a challenge (Sprittles Reference Sprittles2024). Our method does not require us to resolve the variation of velocity and pressure across the gas film, due to the lubrication approach; only variations along it laterally are solved for, which is achieved by tracking the interface, so that all scales in the model can be resolved and, as we shall see, the additional physics fits neatly into the lubrication framework.
Notably, our model contains no parameters fitted to the experiments of drop collisions, in contrast to simulations that merge drops at a predetermined time or that adjust the Hamaker constant to match experimental results (as discussed in Sprittles Reference Sprittles2024).
In this article, in § 2, we describe the model and computational method. Then, in § 3, we quantitatively compare this computational model with experimental values obtained for the critical Weber numbers of the transitions between bouncing and merging found in Huang & Pan (Reference Huang and Pan2021). In § 4, the dependence on gas pressure observed for collisions of tetradecane drops in Qian & Law (Reference Qian and Law1997) is considered, and a regime diagram contour map is developed to easily identify both the critical Weber numbers and, where merging is observed, the type of contact mode. These modes are discussed in § 5. In § 6, we focus on the dependence of the transitions on the drop radius, which is often overlooked. Then, in § 7, the collision of water drops at different pressures is simulated and compared with experiments. Our results lead us to probe very-low-speed collisions in § 8, comparing with theoretical predictions in Bach et al. (Reference Bach, Koch and Gopinath2004), to identify an additional ‘halting transition’ not yet observed experimentally. Finally, we conclude in § 9, and in § 10, future extensions to this work are discussed.
2. Methodology
Consider the head-on collision of two initially spherical liquid drops of radius 
$R$ in a cylindrical coordinate system 
$(r,z,\theta )$ with a plane of symmetry at 
$z=0$ and axisymmetry assumed, so that only one drop needs to be simulated in the 
$(r,z)$-plane. The drops initially move at speed 
$U$ towards the 
$z=0$ plane, so that their relative speed of impact is 
$V=2U$. The base of the upper drop starts from height 
$z=0.1R$, which is far enough for the lubrication pressure to be negligible (this is discussed in Appendix B). Gravity is negligible during the collision process and experiments never consider any dependence on it.
The Navier–Stokes equations are solved for the velocity, 
$\boldsymbol {u}=(u_r,u_z)$, and pressure, 
$p$, in the liquid drop
$$\begin{gather} \rho_l \frac{\mathrm{D}\boldsymbol{u}}{\mathrm{D}t}=\boldsymbol{\nabla} \boldsymbol{\cdot}\boldsymbol{\mathsf{T}}, \end{gather}$$
$$\begin{gather}\boldsymbol{\mathsf{T}}={-} p\boldsymbol{\mathsf{I}}+ \mu_l(\boldsymbol{\nabla}\boldsymbol{u}+\boldsymbol{\nabla}\boldsymbol{u}^\top), \end{gather}$$
$$\begin{gather}0=\boldsymbol{\nabla} \boldsymbol{\cdot}\boldsymbol{u}. \end{gather}$$
Here, the liquid viscosity is 
$\mu _l$ and the stress tensor in the drop is denoted by 
$\boldsymbol{\mathsf{T}}$. Initially, the drop is a sphere of radius 
$R$, centre at 
$(r,z)=(0,1.1R)$, with uniform velocity 
$\boldsymbol {u}(t=0)=-U\boldsymbol {e}_z$.
At 
$r=0$, we apply the symmetry boundary conditions 
$u_r=0$ and 
${\partial u_z}/{\partial r}=0$. On the free surface of the drop, a stress balance boundary condition is applied
\begin{equation} \boldsymbol{n} \boldsymbol{\cdot} \boldsymbol{\mathsf{T}}=\boldsymbol{n} \boldsymbol{\cdot} \boldsymbol{\mathsf{T}}_g-\gamma(\boldsymbol{\nabla}\boldsymbol{\cdot} \boldsymbol{n})\boldsymbol{n}+\frac{A}{6{\rm \pi} h^3}\boldsymbol{n}. \end{equation}
Here, 
$\boldsymbol{\mathsf{T}}_g$ is the stress tensor in the gas and 
$A$ is the Hamaker constant. The entire gas film thickness is 
$h(r,t)$, so that, due to the symmetry at 
$z=0$, the upper free surface is then located at 
$z=h/2$.
The disjoining pressure term in (2.4), which contains the Hamaker constant, accounts for the nanoscale van der Waals attraction between the drops (Pan et al. Reference Pan, Law and Zhou2008). Without this, in this model the drops always bounce (Sprittles Reference Sprittles2024). The Hamaker constant, 
$A$, depends on the liquid considered but can be measured or computed independently of drop-collision experiments (i.e. we do not ‘fit it’). We have used Hamaker constants for each liquid as calculated by Huang & Pan (Reference Huang and Pan2021) using Lifshitz theory (Israelachvili Reference Israelachvili2011).
The liquid's dynamics is coupled to a lubrication equation that solves for the gas pressure, 
$p_g(r,t)$, following the method introduced in Chubynsky et al. (Reference Chubynsky, Belousov, Lockerby and Sprittles2020). The lubrication model is composed of an equation enforcing incompressibility of an evolving gas film and one from the momentum equation for the volume flux, 
$Q_r$, giving
$$\begin{gather} 0=\frac{\partial h}{\partial t}+\frac{1}{r}\frac{\partial (rQ_r)}{\partial r}, \end{gather}$$
$$\begin{gather}Q_r={-}\frac{\varDelta_Ph^3}{12\mu_{g}}\frac{\partial p_g}{\partial r}+hu^d_r. \end{gather}$$
Here, gas viscosity is denoted by 
$\mu _g$, gas-kinetic factor by 
$\varDelta _P(h(r,t))$ (a function of 
$h$, discussed below) and 
$u^d_r=u_r(r,z=h(r,t)/2,t)$ is the component of fluid velocity (on the drop surface) in the radial direction. The lubrication flow has two components: a Poiseuille flow due to the pressure gradient, and a plug flow due to the horizontal motion of the drop surfaces.
The lubrication equation is solved on the lower surface of the drop, and provides normal and tangential stresses that are applied as boundary conditions to the Navier–Stokes equations in (2.4)
\begin{equation} \boldsymbol{n} \boldsymbol{\cdot} \boldsymbol{\mathsf{T}}_g={-}p_g\boldsymbol{n}+ \frac{h}{2}\frac{\partial p_g}{\partial r}\boldsymbol{e}_r. \end{equation}The second term on the right-hand side of (2.7) is the tangential stress, and comes from the Poiseuille flow component of the lubrication flow.
For (2.5) and (2.6), at 
$r=0$, we apply the symmetry boundary condition 
${\partial p_g}/{\partial r}=0$, and at the edge of the film the gas pressure, 
$p_g$, is fixed to a reference gas pressure 
$p_0$. As described in Chubynsky et al. (Reference Chubynsky, Belousov, Lockerby and Sprittles2020), and confirmed by our own tests, the exact position where this edge is set does not affect the results: away from the thin gas film, as 
$h$ becomes large, 
${\partial p_g}/{\partial r}$ will tend to 0, and 
$p_g$ will tend to 
$p_0$ before reaching the exact location of the boundary condition. We set it based on where the gradient 
${\partial h}/{\partial r}$ first exceeds a threshold of 
$10$.
The lubrication model assumes that the gas is incompressible. Analysis of our results shows this is acceptable, with excess pressures usually at most 10 % of atmospheric pressure (this is only violated for the smallest drops at the highest 
$\textit {We}$, which can peak at 25 %). This is in agreement with the estimates for excess pressure in Hicks & Purvis (Reference Hicks and Purvis2011); future analyses could include compressibility as formulated in Mandre et al. (Reference Mandre, Mani and Brenner2009) and Hicks & Purvis (Reference Hicks and Purvis2013).
In (2.6), the term 
$-({h^3}/{12\mu _{g}})({\partial p_g}/{\partial r})$ is the flow rate for plane Poiseuille flow predicted by the Navier–Stokes equations with no-slip boundary conditions. The rarefied-gas adjustment to this term, derived by Chubynsky et al. (Reference Chubynsky, Belousov, Lockerby and Sprittles2020), is the gas-kinetic factor, plotted in figure 2
\begin{equation} \varDelta_P(\textit{Kn})=1+6.88\textit{Kn} +(6\textit{Kn}/{\rm \pi})\ln(1+2.76\textit{Kn}+0.127\textit{Kn}^2), \end{equation}
which is a function of the (local) Knudsen number, 
$\textit {Kn}=\lambda /h(r,t)$, where 
$\lambda$ is the mean free path of the gas.
As 
$\textit {Kn}\to 0$, (2.8) has the asymptotic form (to first order in 
$\textit {Kn}$) of the flow-rate correction found when, instead of no slip, the Navier slip boundary condition is applied: 
$\varDelta _P\sim 1+6.88\textit {Kn}$ (Lauga & Stone Reference Lauga and Stone2003). As 
$\textit {Kn}\to \infty$, it has the correct asymptotic behaviour approaching free molecular flow: 
$\varDelta _P\sim (6\textit {Kn}/{\rm \pi} )\ln (4\textit {Kn}^2/{\rm \pi} ^2)$ (Cercignani et al. Reference Cercignani, Lampis and Lorenzani2004). The asymptotic expressions are shown in figure 2. The combined expression is derived using asymptotic matching and is fitted to within 2 % of solutions of the Boltzmann equation from Cercignani et al. (Reference Cercignani, Lampis and Lorenzani2004) (also shown in figure 2), and so accurately describes the flow rate of Poiseuille flow across the full range of Knudsen numbers. Here, diffuse scattering of molecules from the boundaries is assumed, which is reasonable for most systems. This expression assumes constant temperature, and the simulations in Cercignani et al. (Reference Cercignani, Lampis and Lorenzani2004) use the linearised Bhatnagar–Gross–Krook operator, which is reasonable for monatomic gases or those with small molecules, such as nitrogen.
Li (Reference Li2016) used a similar but less well-matched expression than (2.8), formulated as an effective viscosity, 
$\mu _{g,P}(\textit {Kn})=\mu _g/\varDelta _P$; however, they used this to solve the Navier–Stokes equations in the gas film, rather than exploiting the lubrication approach, as in Chubynsky et al. (Reference Chubynsky, Belousov, Lockerby and Sprittles2020) and here.
For an ideal gas, the mean free path is related to the gas pressure 
$p$, collisional cross-sectional area 
$d^2$ and temperature 
$T$, by
\begin{equation} \lambda=\frac{k_BT}{\sqrt{2}{\rm \pi} d^2p}. \end{equation}
Additionally, assuming that the excess pressure in the gas film is small compared with the ambient pressure, in calculating the local 
$\textit {Kn}$ we can use the mean free path at the ambient pressure 
$p_0$, and assume constant temperature: i.e. we assume that the mean free path is constant along the film and inversely proportional to the ambient pressure. We denote the ratio of the ambient pressure to a reference pressure as 
$\mathcal {P}=p_0/p_{atm}$ (i.e. the ambient pressure relative to standard atmospheric pressure 
$p_{atm}=101$,325 Pa). Then for each gas we can use known values for the mean free path at atmospheric pressure 
$\lambda _{atm}^{gas}$, which differ due to the different collisional cross-sectional area 
$d^2$, to find
\begin{equation} \lambda=\lambda^{gas}_{atm}/\mathcal{P}. \end{equation}
This means that both the gas species and the ambient pressure are captured in our model through the choice of 
$\lambda$. Values used are 
$\lambda _{atm}^{air}=69$ nm, 
$\lambda _{atm}^{nitrogen}=70$ nm and 
$\lambda _{atm}^{helium}=180$ nm.
Our open-source code (available here, DOI: 10.5281/zenodo.10927009) is implemented in the open-source C++ finite element package oomph-lib (Heil et al. Reference Heil, Hazel and Matharu2022). We solve the Navier–Stokes equations (2.1)–(2.3) in the liquid drop on a triangulated mesh, with the gas manifesting itself through the drop's boundary conditions. An implicit, adaptive time stepper is used. We use oomph-lib's pseudo-solid arbitrary-Lagrangian–Eulerian mesh, where the boundary of the mesh is moved with the normal velocity of the liquid at the drop surface, and the interior elements are moved by treating them as an elastic solid. The time derivatives are adjusted to account for the motion of the mesh (Cairncross et al. Reference Cairncross, Schunk, Baer, Rao and Sackinger2000). Periodically the domain is remeshed using an adaptive method (Zienkiewicz & Zhu Reference Zienkiewicz and Zhu1987) that refines the size of the mesh elements as needed to minimise the inter-element stress jump. Here, our focus is purely on predicting the transitions between bouncing and merging. Capturing the post-contact dynamics is tricky but achievable in the current computational framework – it is essentially a routine computational fluid dynamics calculation, involving no additional physics, but is beyond the scope of this work.
A comparison of our computed drop profiles with experiments in Pan et al. (Reference Pan, Law and Zhou2008) are shown in figure 1, although the real power of our model is in being able to predict the transitions between bouncing and merging, rather than just capture the macroscopic motion of the drops.
Figure 1. Comparison of our simulations (red outlines) with the experiments of Pan, Law & Zhou (Reference Pan, Law and Zhou2008) (images have been rotated; snapshots in time from left to right, at intervals of approximately 100 
$\mathrm {\mu }$s, see Pan et al. (Reference Pan, Law and Zhou2008) for exact timings). Examples of drop bouncing near to the soft transition (decreasing impact speed results in coalescence; (a)) and the hard transition (increasing speed results in coalescence; (b)). Physical parameters are: 
$\rho =762$ kg m
$^{-3}$, 
$\mu _l=2.128$ mPa s, 
$\mu _g=18.27$ 
$\mathrm {\mu }$Pa s, 
$\gamma =26.5$ mN m
$^{-1}$, 
$A=5.0\times 10^{-20}$ J; above: 
$R=170.6$ 
$\mathrm {\mu }$m, 
$V=0.486$ m s
$^{-1}$; below: 
$R=167.6$ 
$\mathrm {\mu }$m, 
$V=0.992$ m s
$^{-1}$.
Figure 2. (a) Schematic showing the initial position of the drop. The dashed blue line shows the position of the second drop, which is included via symmetry. (b) The gas-kinetic factor 
$\varDelta _P$ as a function of 
$\textit {Kn}$, showing the simulation data from Cercignani, Lampis & Lorenzani (Reference Cercignani, Lampis and Lorenzani2004) used in Chubynsky et al. (Reference Chubynsky, Belousov, Lockerby and Sprittles2020) in fitting it. The dashed lines show the asymptotic expressions for large and small 
$\textit {Kn}$.
3. Comparison with experimental results of Huang & Pan (Reference Huang and Pan2021)
Table 1 presents quantitative predictions for 
$\textit {We}_S$ and 
$\textit {We}_H$ from our computational model and compares these with experiments from Huang & Pan (Reference Huang and Pan2021). The physical parameters used, taken from Huang & Pan (Reference Huang and Pan2021), are given in the table's caption. Note that we define 
$\textit {We}$ in terms of the drop radius, so that 
$\textit {We}=\rho _lV^2 R/\gamma$; most experimental articles including Huang & Pan (Reference Huang and Pan2021) use the drop diameter, so that the values reported in those articles are twice those referenced here. Huang & Pan (Reference Huang and Pan2021) reported 
$\textit {We}_S$ and 
$\textit {We}_H$ for both the lower-speed transition from merging to bouncing with increasing 
$\textit {We}$ (the ‘soft transition’ at 
$\textit {We}_S$) and the higher-speed transition from bouncing to merging (the ‘hard transition’ at 
$\textit {We}_H$). For comparison, also listed are values found for 
$\textit {We}_S$ and 
$\textit {We}_H$ when the GKE corrections are removed by setting 
$\varDelta _P=1$.
Table 1. Values of the critical Weber numbers, 
$\textit {We}_S$ and 
$\textit {We}_H$, found in our simulations with and without gas-kinetic effect (GKE) corrections, in comparison with the experimental values found in Huang & Pan (Reference Huang and Pan2021). Entries labelled NB (‘no bounce’) show where no transition between merging and bouncing was found, as the drops always merged. The starred entry shows where no soft transition was observed for the values of 
$\textit {We}$ used in the experiments. The simulations results are 
$\pm 5$ of the smallest digit shown (e.g. 0.265 means the value lies in the interval 0.26–0.27; 4.05 means the interval 4.0–4.1).
The physical parameters of the liquids and gas, shown below, are those used in Huang & Pan (Reference Huang and Pan2021). Decane: 
$\rho =728$ kg m
$^{-3}$, 
$\mu _l=0.82$ mPa s, 
$\gamma =23.8$ mN m
$^{-1}$, 
$A=5.0\times 10^{-20}$ J; dodecane: 
$\rho =750$ kg m
$^{-3}$, 
$\mu _l=1.33$ mPa s, 
$\gamma =24.9$ mN m
$^{-1}$, 
$A=5.0\times 10^{-20}$ J; tetradecane: 
$\rho =759$ kg m
$^{-3}$, 
$\mu _l=2.05$ mPa s, 
$\gamma =26.0$ mN m
$^{-1}$, 
$A=5.2\times 10^{-20}$ J; water: 
$\rho =998$ kg m
$^{-3}$, 
$\mu _l=0.92$ mPa s, 
$\gamma =72.0$ mN m
$^{-1}$, 
$A=3.7\times 10^{-20}$ J; air: 
$\mu _g=18.6$ 
$\mathrm {\mu }$Pa s, 
$\lambda =69.0$ nm.
What is immediately apparent is that the simulations which exclude the GKE predict transition values that are wildly different from the experimentally observed ones. Specifically, the bouncing window is massively extended, with values for 
$\textit {We}_S$ often three orders of magnitude too low and 
$\textit {We}_H$ often double what it should be. Without GKE, regimes in which no bouncing (NB) is seen experimentally cannot be reproduced, making this model particularly inept for small drops, where merging becomes the dominant outcome of a collision. In contrast, simulations with GKE routinely predict the correct magnitudes for the transitions and are able to reproduce NB regimes. We now analyse these results in further detail.
For the 
$R=300$ 
$\mathrm {\mu }$m hydrocarbon (decane, dodecane and tetradecane) drops, our simulation results agree very well with the experimental values for 
$\textit {We}_S$ and 
$\textit {We}_H$. Note that no soft transition was observed experimentally for 
$R=300$ 
$\mathrm {\mu }$m tetradecane; this makes sense, as our computed 
$\textit {We}_S$ is below the minimum 
$\textit {We}$ possible in those experiments. Our simulations also agree very well with the hard transition for 
$R=500$ 
$\mathrm {\mu }$m water and 
$R=150$ 
$\mathrm {\mu }$m dodecane drops. For all of the small drops, both our simulations and the experiments found that NB regime can be observed, in contrast to the simulations without the gas-kinetic factor when a bouncing regime was always found, highlighting again that the GKEs are essential to predicting these bounce–merge transitions.
The agreement between the simulations and the experiments is less good for the intermediately sized drops, and, unlike Huang & Pan (Reference Huang and Pan2021), for 
$R=150$ 
$\mathrm {\mu }$m decane and 
$R=115$ 
$\mathrm {\mu }$m tetradecane we observe the NB regime. However, at these radii we predict that 
$\textit {We}_S$ and 
$\textit {We}_H$ are particularly sensitive to radius, as shown in figure 6 and discussed in § 6. Our values for 
$R=115$ 
$\mathrm {\mu }$m tetradecane are very close to the experiments when plotted on the 
$R-\textit {We}$ regime diagram. The agreement is also less good for 
$R=150$ 
$\mathrm {\mu }$m tetradecane, and the soft transitions for 
$R=150$ 
$\mathrm {\mu }$m dodecane and 
$R=500$ 
$\mathrm {\mu }$m water, and we cannot currently explain the cause of these discrepancies. However, the agreement is still far superior to the simulations without GKE, see table 1.
Focusing on a particular set of results, figure 3 shows our simulations for the decane 
$R=300$ 
$\mathrm {\mu }$m drops for 
$\textit {We}$ either side of 
$\textit {We}_S$ and 
$\textit {We}_H$. The main plot shows the minimum film thickness on a log scale, with videos in the Supplemental Material showing the dynamics. Note that we regularly see film thicknesses less than 
$100$ nm, compared with the mean free path of air at atmospheric conditions of 
$69$ nm, justifying the need for gas-kinetic corrections.
Figure 3. Time evolution of minimum film height at 
$\textit {We}$ just above and below 
$\textit {We}_S=0.265$ (blue lines) and 
$\textit {We}_H=6.95$ (orange lines); for 
$R=300$ 
$\mathrm {\mu }$m decane drops. Drop profiles are shown at key stages of the collision for the Weber numbers (
$\textit {We}=0.3$ and 
$\textit {We}=6.5$) in the bouncing regime close to the transition (videos are available as Supplementary Material Movies 1 and 2 available at https://doi.org/10.1017/jfm.2024.722). Below the drop profiles, the lower boundary of the drop is shown, zoomed vertically by a factor of 500, so that the shape of the gas film can be seen. The black dashed line is the critical height (3.1) predicted by Chubynsky et al. (Reference Chubynsky, Belousov, Lockerby and Sprittles2020).
In blue (
$\textit {We}=0.3$, shown in the top panels of figure 3) are snapshots of the simulation for 
$\textit {We}$ just above 
$\textit {We}_S$, with a corresponding video in Supplementary Material Movie 1. For 
$\textit {We}$ just below 
$\textit {We}_S$, contact occurs once the bulk of the drops has starting receding (at the point labelled v in figure 3), as previously found in Bach et al. (Reference Bach, Koch and Gopinath2004). Contact happens at the ‘kink’ at the edge of the film, with a significant gas bubble trapped within.
In orange (
$\textit {We}=6.5$, shown in the bottom panels of figure 3) are snapshots of the simulation for 
$\textit {We}$ just below 
$\textit {We}_H$, with the video given in Supplementary Material Movie 2. Here, the minimum height is reached at the kink at the edge of the gas film (point III in figure 3), and it is here that contact happens for 
$\textit {We}$ just above 
$\textit {We}_H$. The gas film is noticeably affected by the macroscopic behaviour of the drop, such as at point IV in figure 3, where the top of the drop has become concave before rebounding, causing the oscillations in film height seen at point IV. These effects can mostly clearly seen in Movie 2. This film dynamics has not been seen experimentally, unlike for drop-on-solid impacts where interferometry has been able to measure it (de Ruiter, Van Den Ende & Mugele Reference de Ruiter, Van Den Ende and Mugele2015c).
The black line shows the predicted critical film height 
$h_c$, derived in Chubynsky et al. (Reference Chubynsky, Belousov, Lockerby and Sprittles2020), at which the disjoining pressure becomes dominant, overcoming the surface tension, and pulls the interfaces together (i.e. creates contact). This is based on a linear stability analysis of the lubrication equation, which is complicated by the height-dependent gas-kinetic factor 
$\varDelta _P$. Notably, the critical film height does not depend on the impact speed or the drop radius, and is calculated by solution to
\begin{equation} \frac{h_c^5}{\varDelta_P(\lambda/h_c)}=\frac{A^2}{96 {\rm \pi}^2\sigma_{max}\gamma \mu_g}, \end{equation}
where 
$\sigma _{max}$ is an estimated growth rate of the fastest growing mode, which we must impose. In an experiment, the instability will be initiated by a thermal fluctuation, whereas in a simulation it will be a small numerical error. We use an estimate for the size of a thermal fluctuation to approximate this, as described in Chubynsky et al. (Reference Chubynsky, Belousov, Lockerby and Sprittles2020). For decane at room temperature, the characteristic size of interfacial thermal fluctuations is 
$h_{fluc}\sim 0.1$ nm (Sides, Grest & Lacasse Reference Sides, Grest and Lacasse1999). This must grow to cover the gas film of thickness 
$h\sim 20$ nm in a characteristic time of 
$t_{c}\sim 0.1$ ms (these are approximations based on the simulations shown in figure 3). This means that, for an initial perturbation which grows exponentially to span the gap during the collision period, we require 
$\sigma _{max}=\log (h/h_{fluc})/t_c=5\times 10^4$ s
$^{-1}$. The resulting estimate for 
$h_c$ is plotted in figure 3. The theory predicts that, if the minimum film height falls below 
$h_c$, the film becomes linearly unstable and the disjoining pressure can drive the drops into contact, as shown for 
$\textit {We}=0.2$, 0.25, 7.0 and 7.5 in figure 3.
4. Influence of gas pressure
Having validated our approach at atmospheric pressures, we now consider whether our model is able to capture how transitions between bouncing and merging vary with gas type and ambient pressure, for fixed temperature – these effects have only been considered experimentally in Qian & Law (Reference Qian and Law1997). There, tetradecane drops with a radius of 
$100$–
$200$ 
$\mathrm {\mu }$m were considered for nitrogen and helium environments at different pressures. Given the small variation in viscosity for these two gases (differing only by approximately 10 %), we fix 
$\mu _g$ at an intermediate value of 
$18.4$ 
$\mathrm {\mu }$Pa s and run simulations with varied ambient gas pressure, which we recall only appears in the model due to the incorporation of GKEs (imposed via the mean free path) and varied 
$\textit {We}$ (by changing the impact speed, fixing 
$R=150$ 
$\mathrm {\mu }$m).
The inset to figure 4 shows our results for 
$\textit {We}_S$ and 
$\textit {We}_H$ as a curve plotted against the ambient gas pressure, with merging on one side of the curve and bouncing on the other. This shows excellent agreement with the experiments of Qian & Law (Reference Qian and Law1997) and Huang & Pan (Reference Huang and Pan2021), for both nitrogen and helium gases. For a large enough fixed pressure the curve becomes multivalued, showing the merge–bounce–merge transitions. At sufficiently low pressure, drops always merge and here this occurs very near atmospheric pressure for nitrogen (and hence air also). Surprisingly, given good agreement elsewhere, one higher-pressure nitrogen experiment from Qian & Law (Reference Qian and Law1997) does not agree with our simulations, and in fact shows 
$\textit {We}_H$ decreasing with gas pressure. We do not have an explanation for this, and think it is worthy of further experimental analysis: note there are relatively few data points shown in the relevant figure 7(c) in Qian & Law (Reference Qian and Law1997).
Figure 4. Regime diagram for 
$R=150$ 
$\mathrm {\mu }$m tetradecane drops. The Weber number is varied by changing the velocity. The white region is where the drops bounce. The contour shows the relative contact position, 
$r_c/r_{max}$, with labelled contact modes. The inset shows the experimental and computed transitions for the different gases at the gas pressure relative to the reference atmospheric pressure. Above are the drop profiles at the moment of contact for characteristic contact modes, at the numbered points in the regime diagram. Points 1, 2a, 3a and 3b are all for 
$\textit {We}=12$ and with 
$\lambda$ corresponding to 
$\mathcal {P}= 0.6, 1, 1.9, 2.0$ (for air), respectively. Point 2b is 
$\textit {We}=0.5$, 
$\mathcal {P}=1$. Below the drop profiles, the shape of the gas film is shown (zoomed vertically by a factor of 500) so that the position of contact can be seen. The regions of the contour corresponding to each contact mode have been separated by the dashed white lines, and are discussed in §5.
The main section of figure 4 shows 
$\textit {We}_S$ and 
$\textit {We}_H$ against the inverse mean free path. The white region is where the drops bounce, and the coloured region where they merge. The gas type and ambient pressure only appear in the model through the mean free path, calculated from the gas pressure using (2.10). The experimental values have also been plotted in terms of their mean free path (
$\lambda =70$ nm for nitrogen and 
$\lambda =180$ nm for helium), using (2.10) to account for both ambient pressure and gas type. A great boost for the model proposed, with GKEs at its core, is the collapse of experimental data onto a single transition curve showing that the dependence on both ambient pressure and gas species can be characterised purely by the mean free path. As expected, as 
$\mathcal {P}\to \infty$ (or equivalently, 
$\lambda \to 0$, where the mean free path will be much smaller than the gas film heights), 
$\textit {We}_S$ and 
$\textit {We}_H$ tend to those that are found when the gas-kinetic factors are removed.
Previous analyses of drop impact on solids and liquid films have attempted to build regime maps based on the modes of contact (Zhang et al. Reference Zhang, Soori, Rokoni, Kaminski and Sun2021a). Here, by covering all of the relevant parameter space with simulation data (totalling 798 simulations, shown in figure 10), we introduce a new quantitative method for presenting this information and identifying changes in contact mode. This is done by contour plotting relevant measures for the position and time at which contact occurs. For position, in figure 4 the relative contact position is given, i.e. the radial position of contact, 
$r_c$, divided by the maximum radial extent of the drop 
$r_{max}$ at the time of contact, see illustration in figure 4. For time, in figure 5, the normalised contact time indicates how far through the bouncing process contact occurs: the ratio of the contact time, 
$t_c$, to the equivalent time taken for the drops to bounce, 
$t_b$, when the disjoining pressure is removed (this the time until the drops are again at their initial separation distance). Notably, this denominator, 
$t_b$, does not depend on the gas pressure, and is plotted in the inset, along with theoretical predictions for large and small 
$\textit {We}$.
Figure 5. Regime diagram for 
$R=150$ 
$\mathrm {\mu }$m tetradecane drops, as in figure 4, but instead the contour shows normalised contact time, the ratio of the contact time to the equivalent time for the drops to bounce without the disjoining pressure (the time until the drops are again at their initial separation distance). The white region is where the drops bounce, and the contact modes shown in figure 4 are labelled. The bounce time is plotted in red in the inset, where the dashed lines showing asymptotic predictions: for large 
$\textit {We}$ the limit to the natural frequency of the drop 
$(t=2.2(\rho R^3/\gamma )^{1/2})$ (Rayleigh Reference Rayleigh1879), and for small 
$\textit {We}$ the prediction from Gopinath & Koch (Reference Gopinath and Koch2002).
5. Contact modes
The contours in figures 4 and 5 show that there are different contact modes: distinct ways that the drops can come into contact. Characteristic drop and film profiles at the moment of contact are shown at the top of figure 4. Where contact occurs will have a significant effect on how the drops merge beyond the initial contact, and the size of the entrapped gas bubble. There are 3 distinct modes that can be seen, which we call the ‘inner kink’ mode, the ‘outer kink’ mode and the ‘film’ mode. To understand these modes and the boundaries in phase space between them, we can look at the contours in figures 4 and 5 together.
The inner kink mode (with characteristic profile shown at point 1 in figure 4) occurs at high 
$\textit {We}$ and/or low pressure. Contact occurs at a ‘kink’ near the centre of the drop where the local curvature is high, where the entrapped gas bubble meets with an expanding film (as described by Mandre et al. Reference Mandre, Mani and Brenner2009), with this mode having the smallest relative contact position of all the modes. The normalised contact time is also small, as contact happens before the deformed surface of the drop extends beyond the initial radius of the drop. This is the high-speed mode seen for drop-on-solid impact, e.g. in de Ruiter et al. (Reference de Ruiter, Van Den Ende and Mugele2015c), and the ‘film’ mode, seen in Zhang et al. (Reference Zhang, Soori, Rokoni, Kaminski and Sun2021a).
The outer kink mode (with characteristic profile shown at point 2a in figure 4) occurs to the right of the inner kink mode on the regime diagrams. Here, the drop is still ‘spreading’, and contact now occurs at a kink at the edge of the film, where the film meets the more spherical part of the drop. This mode is equivalent to the lower-speed contact seen for drop-on-solid impact in de Ruiter et al. (Reference de Ruiter, Van Den Ende and Mugele2015c), the ‘kink’ mode seen by Zhang et al. (Reference Zhang, Soori, Rokoni, Kaminski and Sun2021a) and the ’kink’ mode found numerically by Chubynsky et al. (Reference Chubynsky, Belousov, Lockerby and Sprittles2020). It is also the mode of the simulations for 
$\textit {We}$ above 
$\textit {We}_H$ in figure 3.
The outer kink mode also extends below the soft transition, with characteristic profile shown at point 2b in figure 4. The normalised contact time is largest here, with contact happening once the drops start rebounding, as shown in figure 3.
There is a clear jump in the relative contact position across the boundary between the outer kink mode and the inner kink mode. In this part of the regime diagram, during the part of the collision when contact occurs there are two local minima in the film: the two ‘kinks’ where contact occurs for the inner and outer kink modes. This can be seen most clearly in the zoomed gas film profile at point 2a in figure 4. The boundary between the kink modes is formed due to a change in the kink which has minimum film height at the moment of contact. A flip in where this minima is can be seen in figure 3: the snapshots of the 
$\textit {We}=6.5$ simulation show that at point I the minima is relatively near the centre of the drop (at the inner kink), but that by point II, the minima is at the outer kink. For drop-on-solid impacts this two-kink behaviour was observed experimentally and discussed in detail by de Ruiter et al. (Reference de Ruiter, Van Den Ende and Mugele2015c), and contact occurring at both an inner and an outer kink is seen by Li et al. (Reference Li, Langley, Tian, Hicks and Thoroddsen2017).
The film mode has two characteristic contact profiles shown at points 3a and 3b in figure 4. Here, the drop has spread significantly, and contact can occur either in the middle of the gas film (3a) or at the edge (3b). This is the ‘film mode’ at high pressures seen in Chubynsky et al. (Reference Chubynsky, Belousov, Lockerby and Sprittles2020), and appears similar to the contacts seen in impacts of ultra-viscous drops in Langley et al. (Reference Langley, Li and Thoroddsen2017), where an extended air film formed before contact occurred. The exact position of contact depends on where numerical noise triggers the disjoining pressure driven instability, but in reality thermal fluctuations could trigger this film collapse (Fetzer et al. Reference Fetzer, Rauscher, Seemann, Jacobs and Mecke2007; Sprittles et al. Reference Sprittles, Liu, Lockerby and Grafke2023). This is why the film mode region of the contour in figure 4 is ‘speckled’, i.e. non-uniform in colour: each simulation has made contact at a different noise-derived point, which is sensitive to the precise numerical set-up. From a practical perspective, we may expect the position of contact seen in experiments to be stochastic in nature in this region, as is indeed the case for the similar drop–solid experiments in Langley et al. (Reference Langley, Li and Thoroddsen2017).
Across the boundary between the outer kink mode and the film mode, there is a jump in the normalised contact time, most clearly at the moderate 
$\textit {We}$ where the film mode first appears. This is because when the outer kink contact disappears, the drop spends time spreading before the film mode collapse can occur. This also highlights that having information from both temporal and spatial information for the contact helps define modes in a less ambiguous manner.
For drop-on-solid impact, Chubynsky et al. (Reference Chubynsky, Belousov, Lockerby and Sprittles2020) also observed contact at the very end of the bouncing process, during the collapse of the dimple where gas is trapped between the drops (as shown at points vii and VII in figure 3). This occurred for high-pressure, low-speed impacts, but we have not observed this in our drop-collision simulations. It is yet to be determined if this is a fundamental difference between drop collisions and drop-on-solid impacts, or is simply due to the different parameter regimes of the studies.
6. Dependence on drop size
The majority of studies have focused on how transitions depend on impact speed, and have thought of 
$\textit {We}$ as a dimensionless measure of this. Here, to systematically study the dependence of 
$\textit {We}_S$ and 
$\textit {We}_H$ on drop radius, considered experimentally in Huang & Pan (Reference Huang and Pan2021), we simulate the collision of tetradecane drops at atmospheric pressure (fixed 
$\lambda$) with varying drop radius and impact speed. Figure 6 shows the resulting regime map in terms of 
$\textit {We}$ and drop radius, with the contour showing relative contact position as in figure 4. It shows very good agreement with the experimental results of Huang & Pan (Reference Huang and Pan2021), and that for a fixed 
$\textit {We}$ the drops can bounce or merge depending on 
$R$.
Figure 6. Above: regime diagram for head-on collisions of tetradecane drops at atmospheric pressure (and therefore fixed 
$\lambda$). For a particular radius, 
$\textit {We}$ is varied by changing only the impact velocity. The contour shows the relative contact as in figure 4. The dashed red line shows 
$\textit {We}_S$ and 
$\textit {We}_H$ when GKEs are not considered. Below: the lower section of the above plot, but with a log scale on the y-axis to show the soft transitions more clearly. Also included is the merging–halting transition discussed in § 8.
In figure 6, the boundary between the inner kink (light region) and outer kink modes (dark region) can clearly be seen. The film mode does not appear here, and was only observed for higher ambient pressures than the atmospheric pressure used to create this figure.
The dashed red lines in figure 6 show what happens to the hard transition when GKEs are removed, by setting 
$\varDelta _P=1$. This curve is dramatically different – for 
$R=150$ 
$\mathrm {\mu }$m, it is well over twice the experimental value for 
$\textit {We}_H$, and increases as the radius decreases, despite the experiments and our primary simulations showing that the bouncing regime vanishes for small drops. Moreover, without GKEs, 
$\textit {We}_S$ is below 
$\textit {We}=0.001$, far beneath the soft transitions found experimentally or by the simulations with GKE, as shown in table 1 and the bottom plot of figure 6.
As 
$R$ becomes large, 
$\textit {We}_H$ tends to a constant value. For simulations up to 1 mm, 
$\textit {We}_H$ tends to 
$10.5\pm 0.5$, and without GKEs, 
$\textit {We}_H$ tends to 
$12.5\pm 0.5$. Here contact occurs in the film mode, once the drop has spread significantly. Contact is determined by the stability of the extended gas film. If the height of the film is below the critical height for the van der Waals driven instability, contact will occur. Our results suggest that (for large 
$R$ at least) the stability of the film depends on the radius only via the Weber number. The height (and stability) also depend on the mean free path, and hence the values of 
$\textit {We}_H$ with and without GKEs do not converge for large 
$R$.
The transition boundary from figure 6 is plotted again in figure 7, this time in terms of velocity and drop radius. This reveals that for some values of the velocity (here, around 
$V=1$ m s
$^{-1}$) there is also a merge–bounce–merge transition with increasing radius (see the arrow in figure 7). The experiments of Huang & Pan (Reference Huang and Pan2021) also demonstrate this, but they only have 3 data points (as shown in figure 7), and did not present this in terms of velocity. By systematically varying the radius and presenting the results in this way, we are able to identify this merge–bounce–merge transition with fixed velocity.
Figure 7. The regime diagram for tetradecane, as in figure 6, but in dimensional variables, showing critical velocities as a function of drop radius, at atmospheric pressure. The blue arrow shows that for fixed velocity, there is an additional merge–bounce–merge transition.
Both the transitions found with fixed radius (in figures 4 and 5) and with fixed 
$\lambda$ (figure 6) could be plotted against an inverse global Knudsen number 
$R/\lambda$. However, the curves are similar but not the same, demonstrating that 
$\textit {We}_S$ and 
$\textit {We}_H$ are not simply a function of a global Knudsen number of the ratio of 
$\lambda$ and 
$R$.
The theory of Chubynsky et al. (Reference Chubynsky, Belousov, Lockerby and Sprittles2020) (the solution of (3.1)) tells us whether or not the film is stable at a given height; however, without being able to predict this height for a given set of parameters it cannot tell us a priori whether or not a given pair of colliding drops will merge or bounce; clearly such a theory would be highly sought after. Some previous analytic models, such as those of those of Mandre et al. (Reference Mandre, Mani and Brenner2009), Huang & Pan (Reference Huang and Pan2021) and Sharma & Dixit (Reference Sharma and Dixit2021), have predicted the minimum film heights, but these are yet to incorporate full GKEs which we have seen is essential. The model of Zhang & Law (Reference Zhang and Law2011) includes a gas-kinetic factor and predicts both 
$\textit {We}_S$ and 
$\textit {We}_H$, but it uses a parameter fitted to experimental results and requires solving a nonlinear ordinary differential equation. Bach et al. (Reference Bach, Koch and Gopinath2004) considered 
$\textit {We}_S$ for small drops with small deformations, predicted a minimum height incorporating Knudsen effects, and used the results of their simulations to fit a critical height, comparing favourably with the experiments of Qian & Law (Reference Qian and Law1997).
An ideal model for predicting the bounce–merge transitions would characterise 
$\textit {We}_S$ and 
$\textit {We}_H$ as functions of a rarefaction parameter 
$\lambda /H$, where 
$H$ is a characteristic length scale of the gas film height. Even then, however, figures 4 and 5 show us how complex the situation is, with different modes of contact that would have to be predicted by any simplified theory hoping to understand the instability mechanism. Furthermore, such a model should predict both soft and hard transitions, meaning most simple expressions will fail.
García-Geijo, Riboux & Gordillo (Reference García-Geijo, Riboux and Gordillo2024) compared several models for the minimum height of the film during the initial approach of a drop towards a solid, but these theories could easily be adapted to drop–drop collisions. Combined with a theory for the critical height at which the disjoining pressure causes contact, this could be a step towards such a model. However, at present, they only predict the height of the ‘inner kink’, and we have shown that (at least in the parameter regime considered), the bounce–merge transitions are between contact at later times, such as at the ‘outer kink’.
7. Water drops
Similar to the results for tetradecane in § 4, we consider the bounce–merge transitions for water in figure 8. Again, we plot 
$\textit {We}_S$ and 
$\textit {We}_H$ against the inverse mean free path, and compare with the experiments of Qian & Law (Reference Qian and Law1997) for different gases and ambient pressures. These experiments varied 
$\textit {We}$ by varying both the drop radius (between 100 and 200 
$\mathrm {\mu }$m) and the impact speed, and the results are reported only in terms of 
$\textit {We}$. As with the hydrocarbons discussed above, for a fixed value of 
$\textit {We}$, the drops can bounce or merge depending on the drop radius. We therefore vary 
$\textit {We}$ by only changing the impact speed, and run two sets of simulations: for 
$R=100$ 
$\mathrm {\mu }$m and 
$R=150$ 
$\mathrm {\mu }$m.
Figure 8. Computed transitions for water drops collision of two different radii, compared with the experiments of Qian & Law (Reference Qian and Law1997).
At lower pressure, the experimentally reported transitions lie between the two bounce–merge transition curves found from our model, demonstrating the sensitivity to the drop radius, and the need for experiments to report the radius and velocity independently. The agreement is less good for the higher-pressure (larger 
$\lambda ^{-1}$) experiments, similar to the high-pressure tetradecane in nitrogen experiment discussed in § 4. Again, there is a large amount of experimental uncertainty in the relevant figures 6(c) and 8(c) of Qian & Law (Reference Qian and Law1997), and this regime is clearly worthy of further experimental study.
8. Very-low-speed impacts
For sufficiently low impact speeds, there is an additional transition that cannot be seen in the figures above. Here, the lubrication pressure is enough to slow the drops until they are stationary. This is shown for tetradecane in figure 9, with halting occurring between 
$\textit {We}=10^{-6}$ and 
$\textit {We}=10^{-5}$. A critical Weber number curve for this transition is also shown in the bottom part of figure 6. This was predicted theoretically by Bach et al. (Reference Bach, Koch and Gopinath2004), who suggest a value of 
$\textit {We}=9.8\times 10^{-5}$ based on a fit for the coefficient of restitution. Note that our simulations include only the lubrication pressure from the gas film, not the Stokes drag from the surrounding air which would cause halting to occur for higher 
$\textit {We}$ than seen here. The velocities here are far smaller than those investigated experimentally with free drops, and so this transition has not yet been observed. Altogether then, there can be three transitions in the outcome of drop collisions with increasing 
$\textit {We}$: starting from 
$\textit {We}=0$, with increasing 
$\textit {We}$ the drops first fail to merge (‘halting’), then merge, then bounce, then merge again!
Figure 9. Very-low-speed 
$R=150$ 
$\mathrm {\mu }$m tetradecane collisions at atmospheric pressure, showing the transition from merging to halting.
9. Conclusions
This article has presented a computational model for the head-on collision of identical liquid drops in gases of different ambient pressures, with the aim of predicting the conditions under which the drops will merge rather than bounce, due to the presence of a trapped gas film. We compare favourably with experiments, and have constructed regime diagrams showing how the outcome of a drop collision varies with drop size and ambient gas pressure.
We have seen that the conventional model of a continuum, incompressible lubricating gas film completely fails to describe the merge–bounce transitions in drop collisions, and nanoscale effects are essential to the accurate modelling of these collisions. Our new computational model, incorporating these nano-effects, is able to efficiently capture experimental observations with relatively good quantitative accuracy across a wide range of parameter space.
We have shown that the critical Weber numbers can be surprisingly sensitive to drop radius, and that for fixed velocity, there is a merge–bounce–merge with increasing drop radius. Finally, we have identified that for very slow drops, the lubrication pressure can cause the drops to halt, rather than merge.
10. Future directions
The results in this article suggest some regimes that would be of particular interest for future experimental investigation. In particular, experiments with controlled drop size at different ambient pressures, looking at: (i) the critical pressure below which the bouncing regime disappears, and (ii) the trend towards constant 
$\textit {We}_H$ for high ambient pressure.
In this article we consider only head-on collisions. Notably, for off-centre collision some models for 
$\textit {We}_S$ and 
$\textit {We}_H$ are fitted through the values for the head-on transitions (Hu et al. Reference Hu, Xia, Li and Wu2017; Al-Dirawi & Bayly Reference Al-Dirawi and Bayly2019). Clearly, however, to capture off-centre collisions correctly, the computational model must become three-dimensional, with the lubrication equations becoming two-dimensional. Here, it is the dramatic increase in computational cost that is the main limiting factor, but our lubrication approach should help to reduce this burden and the extension to three dimensions should be the focus of important future work.
Our computations are currently only valid up to the moment of contact, and it would be necessary to perform ‘numerical surgery’ on the mesh to continue the simulation beyond that point (Anthony et al. Reference Anthony2023). This is certainly possible and will be an avenue for future work. An alternative approach would be to use our results at the moment of contact as an initial condition to a solver better suited to simulations involving changes in topology, but that cannot predict when and where the initial contact will occur.
The model described in this paper can also easily be applied to impacts of drops onto smooth solids, as in Chubynsky et al. (Reference Chubynsky, Belousov, Lockerby and Sprittles2020). Recent experiments have been able to measure the thickness of the gas layer with remarkable accuracy (Kaviani & Kolinski Reference Kaviani and Kolinski2023), and comparisons could easily be made using this model.
Beyond drop–drop collisions and drop-on-solid impacts, there are good prospects of including additional physical processes. Drop impacts onto films could be considered with the addition of a second lubrication equation for the liquid film (Duchemin & Josserand Reference Duchemin and Josserand2020). Thin liquid films have been used in experiments as an analogue for smooth solids, but unexplained discrepancies between the behaviours of solid and films have been identified (Lo et al. Reference Lo, Liu and Xu2017; Pack et al. Reference Pack, Hu, Kim, Zheng, Stone and Sun2017; Zhang et al. Reference Zhang, Soori, Rokoni, Kaminski and Sun2021b), with different contact modes observed. By incorporating a model for a liquid film, our methodology could isolate the differences of a film as compared with a smooth solid.
The impact of a drop onto a liquid bath could also be modelled, using a second meshed domain for the bath solving the Navier–Stokes equations, along with the required adjustments to the lubrication equation. In Tang et al. (Reference Tang, Saha, Law and Sun2019), the dynamics of the trapped gas film was measured in an impact with a bath of thickness similar to the drop radius. Our methodology could investigate these dynamics, identify contact modes and predict the critical Weber numbers for merging to occur. By combining this approach with a liquid thin-film model, a regime diagram of Weber number vs bath thickness could be constructed.
As investigated by Pan et al. (Reference Pan, Law and Zhou2008), the critical Weber numbers are also sensitive to the value of the Hamaker constant. We have used Hamaker constants for each liquid as calculated by Huang & Pan (Reference Huang and Pan2021) using Lifshitz theory (Israelachvili Reference Israelachvili2011), and it is not fitted to any drop-collision results. This constant will be different for each pair of contacting materials, which must be considered carefully when comparing drop-on-solid and drop-on-film impacts.
If a drop is placed above a solid surface which is hotter than the liquid's boiling temperature, it can levitate on a layer of its own evaporating vapour, in the ‘Leidenfrost effect’ (Quéré Reference Quéré2013). The model described in this article has previously been extended to quasi-static drops in Leidenfrost conditions (without GKEs), and used to predict the equilibrium shapes of these drops (Chakraborty, Chubynsky & Sprittles Reference Chakraborty, Chubynsky and Sprittles2022). Transitions between bouncing and contact have been observed for impacts in Leidenfrost conditions (Tran et al. Reference Tran, Staat, Prosperetti, Sun and Lohse2012; Chantelot & Lohse Reference Chantelot and Lohse2021), though at higher Weber numbers than for impacts in isothermal conditions. Leidenfrost drops that are statically levitating have also been shown to start spontaneously ‘trampolining’ (Liu & Tran Reference Liu and Tran2020; Graeber et al. Reference Graeber, Regulagadda, Hodel, Küttel, Landolf, Schutzius and Poulikakos2021). To investigate these dynamic processes, the model could be coupled to heat flow inside the drop, along with a more complex evaporation model. This would allow the thermal Marangoni effect to be included, which has been predicted to have a significant effect on the stability of Leidenfrost drops (van Limbeek et al. Reference van Limbeek, Ramírez-Soto, Prosperetti and Lohse2021).
Supplementary movies
Supplementary movies are available at https://doi.org/10.1017/jfm.2024.722.
Funding
This work was supported by the EPSRC under grants EP/W031426/1, EP/S022848/1, EP/S029966/1, EP/P031684/1 and EP/V012002/1. P.L.-J. is supported by a studentship within the UK EPSRC-supported Centre for Doctoral Training in the Modelling of Heterogeneous Systems (HetSys), EP/S022848/1. For the purpose of open access, the author has applied a CC BY public copyright licence to any Author Accepted Manuscript version arising from this submission.
Declaration of interests
The authors report no conflict of interest.
Data availability statement
The data that support the findings of this study are openly available in figshare at http://doi.org/10.6084/m9.figshare.24894237.
Appendix A. Sampling points
In figure 10 we present the sampling points used in figures 4 and 5.
Figure 11. (a) The force on a sphere approaching a symmetry plane. The sphere has radius 
$300\,\mathrm {\mu }$m like the decane drops in figure 3, and the air has the properties of that of figure 3. Shown is the actual (‘full’) force on the spheres with a Navier slip condition as derived by Reed & Morrison (Reference Reed and Morrison1974), the force derived from the lubrication approximation for the air with slip due to Hocking (Reference Hocking1973), lubrication without slip and the limit for large separation. All forces are relative to Stokes drag 
$-6{\rm \pi} \mu _g R U$. The dashed black line shows the initial separation of the drops in our simulations. The dotted black line shows the separation when the drops start to deform for the 
$\textit {We}=0.3$ simulation. (b,c) The resulting motion of a sphere due to these forces, between 
$h=0.2$ and 
$h=0.01$.
Appendix B. Validity of assumptions
Two assumptions we make are in (i) neglecting the gas's force on the drop outside of the lubricating gas film and (ii) assuming our results are relatively insensitive to small changes in the initial condition. To investigate these assumptions we will consider known results for the forces on head-on approaching spheres in Stokes flow (applicable for our impact speeds), with Navier slip to account for GKEs in the ‘slip regime’ (
$Kn\ll 1$). For sake of argument, we take the case of 
$\textit {We}=0.3$ and 
$R=300$ 
$\mathrm {\mu }$m decane drops from figure 3. Specifically, we compare the forces on the spheres due only to the lubricating gas film (with Navier slip on the gas-sphere surface) as computed by Hocking (Reference Hocking1973) with the full force (again with Navier slip) as computed by Reed & Morrison (Reference Reed and Morrison1974). For drops much smaller than the ones considered in this paper, these assumptions will no longer hold (Poydenot & Andreotti Reference Poydenot and Andreotti2024).
The aforementioned forces are plotted in figure 11(a), along with the force due to lubrication without slip and the Stokes drag (with slip) on the spheres when they are far apart. In the inset, it is shown that the lubricating force is already over half of the total force at the initial condition of the simulation (shown at the dashed black line). The dotted black line shows the separation (
$h=0.01$) when the drops start to noticeably deform (for the 
$\textit {We}=0.3$, 
$R=300$ 
$\mathrm {\mu }$m decane case from figure 3), at which point the lubrication force is over 90 % of the total. Once the drops start to deform the lubrication force will become even larger. This confirms that assumption (i) is reasonable and that the lubrication force is capturing the majority of the force just before the collision takes place.
Justification for (ii) and further justification for (i) come from studying the effect on the spheres’ speed in this period. The resulting motion of the sphere between 
$h=0.2$ (initial condition) and 
$h=0.01$ (first deformation) due to the different forces is plotted in figure 11(b), with the resulting velocities shown in figure 11(c). Notably, whilst the drop decelerates slightly in this period, this is minimal for both the full force case (where the effective Weber number is reduced from 0.3 to 
$\textit {We}=0.2979$) and the lubrication-only case (
$\textit {We}=0.2985$). These differences are all so small that they will not affect the outcome of the collision, and we are justified in both assumptions (i) and (ii).
In future, it would be interesting to compare these differences in more detail across a range of parameter regimes (e.g. including liquid–liquid systems), comparing computations from our model with those of the two-phase Navier–Stokes system that incorporate drop deformation, inertial effects and slip. We could also then evaluate the range of applicability of recent simplified models for drop collisions, e.g. in the context of predicting precipitation (Poydenot & Andreotti Reference Poydenot and Andreotti2024).
Open access
