1. Introduction
Coalescence of fluid drops is a central process in many natural, industrial and medical settings, including raindrop formation (Bowen Reference Bowen1950), ink-jet printing (Stringer & Derby Reference Stringer and Derby2009) and drug delivery (Schwarz et al. Reference Schwarz, Mehnert, Lucks and Müller1994). Consequently, the fundamental problem of the coalescence of two isolated drops has been the the focus of a large body of work, with the many different dynamical regimes of coalescence being studied theoretically (e.g. Hopper Reference Hopper1984; Eggers, Lister & Stone Reference Eggers, Lister and Stone1999; Duchemin, Eggers & Josserand Reference Duchemin, Eggers and Josserand2003), numerically (e.g. Baroudi, Kawaji & Lee Reference Baroudi, Kawaji and Lee2014; Anthony, Harris & Basaran Reference Anthony, Harris and Basaran2020) and experimentally (e.g. Paulsen, Burton & Nagel Reference Paulsen, Burton and Nagel2011; Chireux, Tordjeman & Risso Reference Chireux, Tordjeman and Risso2021).
The early-time rate of drop coalescence is strongly influenced by the initial surface profile of the drops. For example, theoretical models of coalescence for viscous drops in an inviscid exterior fluid often take the initial drop geometry to be two spheres touching at a point (e.g. Hopper Reference Hopper1984; Eggers et al. Reference Eggers, Lister and Stone1999), which gives an early-time rate of coalescence 
$\sim - (\log t)/{\rm \pi}$ for time 
$t$. This result, which has been confirmed numerically with full Navier–Stokes simulations (e.g. Sprittles & Shikhmurzaev Reference Sprittles and Shikhmurzaev2014; Anthony et al. Reference Anthony, Harris and Basaran2020), is faster than, and has a different time dependence to, the apparently constant early-time rate of coalescence that has been observed in experiments (e.g. Paulsen et al. Reference Paulsen, Burton, Nagel, Appathurai, Harris and Basaran2012; Chireux et al. Reference Chireux, Tordjeman and Risso2021). One possible reason for this discrepancy might be a finite separation between the drops at the onset of coalescence, with some local surface deformation leading to (local) contact: Anthony et al. (Reference Anthony, Harris and Basaran2020) showed numerically that such initial conditions can give rise to a nearly constant early-time rate of coalescence, and Paulsen (Reference Paulsen2013) suggests that there may have been a finite initial separation in his experiments. Due to the micrometre scales of the initial contact region, it is difficult to observe directly the initial conditions in experiments.
The drop profile at the onset of coalescence is determined by any deformation of the drops prior to contact. At contact, there is a topological transition between disjoint and connected drops. Singular problems such as the approach to transition often exhibit self-similar behaviour (see Eggers & Fontelos Reference Eggers and Fontelos2005), in which case we would expect a generic surface profile prior to contact that gives the surface profile at the onset of coalescence. For example, the similar topological transition for capillary pinch-off of a viscous drop has self-similar dynamics prior to pinch-off that sets the conical initial shapes of the resulting drops at the onset of subsequent recoil (Lister & Stone Reference Lister and Stone1998).
Experiments (Crassous et al. Reference Crassous, Charlaix, Gayvallet and Loubet1993) have shown that for drops that are effectively stationary with any bulk motion occurring on a time scale much greater that the time scale of local deformation, the drop surfaces deform due to van der Waals intermolecular attraction across the gap between the drops. If the gap width is sufficiently small, then the van der Waals attraction overcomes stabilizing effects, such as surface tension, inertia and viscosity, to rapidly deform the drop surfaces into contact. This process is often known as jump-to-contact. When the drops are surrounded by a viscous exterior fluid, the dynamics is controlled by the thinning intervening sheet (Jones & Wilson Reference Jones and Wilson1978; Davis, Schonberg & Rallison Reference Davis, Schonberg and Rallison1989; Jiang & James Reference Jiang and James2007; Chan, Klaseboer & Manica Reference Chan, Klaseboer and Manica2011; Duchemin & Josserand Reference Duchemin and Josserand2020). When the sheet is much more viscous than the drops, the thinning of the sheet is self-similar, with van der Waals attraction, surface tension and viscous forces all acting at leading order (Zhang & Lister Reference Zhang and Lister1999; Vaynblat, Lister & Witelski Reference Vaynblat, Lister and Witelski2001; Moreno-Boza, Martínez-Calvo & Sevilla Reference Moreno-Boza, Martínez-Calvo and Sevilla2020). However, such analyses are not appropriate when the exterior fluid is effectively inviscid, including for jump-to-contact between liquid drops in air (e.g. Chireux et al. Reference Chireux, Protat, Risso, Ondarçuhu and Tordjeman2018), and in atomic force microscopy of a fluid, where understanding jump-to-contact is essential to prevent the wetting of the solid probe (e.g. Ledesma-Alonso, Legendre & Tordjeman Reference Ledesma-Alonso, Legendre and Tordjeman2012; Quinn, Feng & Stone Reference Quinn, Feng and Stone2013).
For two very viscous drops in an effectively inviscid exterior fluid Beaty & Lister (Reference Beaty and Lister2022) (referred to henceforth as BL22) found that jump-to-contact is self-similar. The dynamics is driven by the van der Waals attraction in a small neighbourhood of the tip, inducing a flow in the drop that is resisted by viscosity. Surface tension does not act at leading order, consistent with the description of jump-to-contact as occurring when the van der Waals attraction is too great to be constrained by surface tension. The dynamics remains viscously dominated throughout jump-to-contact for drops with Ohnesorge number 
$Oh \equiv \mu / (\rho \gamma a)^{1/2}\gtrsim 1$, where 
$\mu$ is the dynamic viscosity, 
$\rho$ is the density, 
$\gamma$ is the surface tension, and 
$a$ is the radius of the unperturbed spherical drop profile. However, for only moderately viscous fluid drops, 
$Oh \ll 1$, inertia becomes significant before contact is reached.
In this paper, we elucidate the various dynamical regimes for jump-to-contact between two drops in an effectively inviscid exterior fluid by including the inertia of the fluid within the drops. There are distinct Stokes and inertial–viscous regimes in which either viscous forces or both inertial and viscous forces are balanced with the driving van der Waals force at leading order. As contact is approached, the dynamics is then given, respectively, either by the viscous self-similar solution of BL22, or by a new self-similar solution in the inertial–viscous regime. We argue that there is no distinct inertial regime.
In § 2, we formulate the problem, obtaining the general governing equations for jump-to-contact. We also simplify the equations under a shallow-slope approximation that is shown to hold initially. In § 3, we discuss the possible leading-order force balances and scaling relations, and include a recap of numerical results for the Stokes regime studied in our previous paper (BL22). We identify a possible new balance between inertial, viscous and van der Waals forces, and predict its similarity scalings. In § 4, we perform time-dependent simulations of jump-to-contact with inertia, and demonstrate that the numerical results confirm the new predicted scalings. We conclude in § 5, where we give a full description of the different dynamical regimes for jump-to-contact and the dynamics in each regime. We also discuss how the results suggest new initial conditions for the coalescence of two drops that differ from the condition of touching spheres often used in models of coalescence.
2. Formulation of the problem
We consider the jump-to-contact dynamics between two identical drops of Newtonian fluid in an effectively inviscid surrounding fluid. The drops are assumed to be sufficiently close to each other for the van der Waals attraction between them to be significant. The van der Waals attraction between the opposing drops draws them together, deforming them in the process, until the separation between the drops is comparable to molecular length scales, at which point the drops have effectively reached contact. Deformation of the drop profile is resisted by surface tension. Throughout this paper, we assume that the unperturbed drop radius 
$a$ is much greater than the characteristic length scale of the van der Waals attraction. The deformation to the drop radius due to the van der Waals attraction is then small compared to the unperturbed radius and localized to some ‘near-contact’ region. We focus on the deformation to the drops in the near-contact region, and neglect any bulk motion due to the net attractive van der Waals force. In the experiments that motivate this study (e.g. Paulsen et al. Reference Paulsen, Burton, Nagel, Appathurai, Harris and Basaran2012; Chireux et al. Reference Chireux, Protat, Risso, Ondarçuhu and Tordjeman2018), the bulk motion is resisted by some apparatus holding the drops on the macroscopic scale.
As the exterior fluid is effectively inviscid, there is no resistance to the contact dynamics from the gap. In this case, contact is made at a single, central point determined by the initial point of closest approach (Jiang & James Reference Jiang and James2006). We thus consider axisymmetric point contact, as depicted in figure 1. We take a cylindrical coordinate system with the plane 
${z}=0$ equidistant from both drops, the 
${z}$ axis through the centres of the drops, and a radial coordinate 
${r}$. The system has both up–down symmetry about 
${z}=0$ and axisymmetry about 
${r}=0$. The total vertical gap width between the drops at some radius is 
$2 h( r , t) \ll a$, with the surfaces of the drops at 
${z} = \pm {h}$. We will henceforth consider the evolution of the surface 
${z}={h}$ by utilizing the symmetry of the system.
Figure 1. Two drops with surfaces 
$z = \pm h(r,t)$ attracted by the van der Waals force across an inviscid gap.
The van der Waals attraction between the two drops arises from pairwise interactions between molecules across the gap. The force experienced by one drop due to the presence of the other can be expressed as an effective pressure 
$\varPi$, known as the ‘disjoining pressure’, exerted by the gap on the surface of the drop. Thus the total forcing on the drop due to the van der Waals attraction, surface tension and pressure in the gap is given by an applied normal stress on the drop surface, expressed as a dynamic boundary condition
\begin{equation} \boldsymbol{\sigma}\boldsymbol{\cdot}\boldsymbol{n} ={-}(p_e+ \gamma \kappa + \varPi)\boldsymbol{n}\quad \mbox{on}\ z=h , \end{equation}
where 
$\boldsymbol {{\sigma }}$ is the stress tensor for the flow in the drop, 
$\boldsymbol {n}$ is the outward surface normal, 
${p}_e$ is the pressure in the gap, 
$\gamma$ is the surface tension, and 
${\kappa }$ is the surface curvature. As the gap is effectively inviscid, 
${p}_e$ is constant and so can be set to zero.
The van der Waals attraction drives the dynamic jump-to-contact process through gradients of the normal stress on the surface (2.1) which force a flow in the drop. The flow satisfies the Navier–Stokes equations
\begin{equation} \rho \frac{\partial \boldsymbol{u}}{\partial t} + \rho\boldsymbol{u} \boldsymbol{\cdot} \boldsymbol{\nabla} \boldsymbol{u} ={-}\boldsymbol{\nabla} p+\mu \nabla^2 \boldsymbol{u},\quad \boldsymbol{\nabla} \boldsymbol{\cdot}\boldsymbol{u}=0. \end{equation}The gap between the drops evolves according to the velocity on the surface of the drops, as given by the kinematic boundary condition
\begin{equation} \frac{\partial h}{\partial t} + u_r\,\frac{\partial h}{\partial r}= u_z \quad \mbox{on}\ z = h, \end{equation}
where 
$u_r$ and 
$u_z$ are the radial and vertical velocity components, respectively.
Conversely, when the normal stress on the surface is constant, there is no flow in the drop and thus no surface evolution. This is consistent with experiments (e.g. Chireux et al. Reference Chireux, Protat, Risso, Ondarçuhu and Tordjeman2018) that reveal that jump-to-contact between drops is preceded by a regime where the drop surfaces are locally in a static equilibrium. In equilibrium, the gap profile is determined by a balance between surface tension and van der Waals attraction in (2.1), with the constant value of the normal stress set on the macroscopic scale by the capillary pressure 
$2\gamma /a$ due to a spherical far-field shape. Analytical and numerical studies (Ledesma-Alonso et al. Reference Ledesma-Alonso, Legendre and Tordjeman2012; Quinn et al. Reference Quinn, Feng and Stone2013; Hilaire et al. Reference Hilaire, Siboulet, Ledesma-Alonso, Legendre, Tordjeman, Charton and Dufrêche2020) show that such an equilibrium exists when the separation between the drops is above a critical value. Below the critical separation, no equilibrium is possible as the van der Waals attraction is sufficiently large that surface tension is unable to restrain it. Motivated by these experimental and theoretical results, we consider jump-to-contact as an instability between two drops where initially the van der Waals attraction and surface tension are approximately in balance.
As jump-to-contact begins when the equilibrium balance fails, the initial length scales of the near-contact region are given by the equilibrium balance between the van der Waals attraction, surface tension and the constant capillary pressure set by the far field. The van der Waals disjoining pressure between two surfaces separated by a gap of width 
$2h$ scales as 
$A/h^3$ (Hamaker Reference Hamaker1937), where 
$A$ is the Hamaker coefficient. Here and later in this paper, we use 
$h$, 
$r$ and other terms in scaling arguments to denote the typical magnitudes of the respective variables. The disjoining pressure in equilibrium is in balance with the surface tension in the deformed region 
${\sim }\gamma h/r^2$ and the capillary pressure 
${\sim }\gamma /a$, giving
\begin{equation} h \sim a \mathcal{H}^{1/3}, \quad r \sim a \mathcal{H}^{1/6}, \end{equation}where the non-dimensional parameter
\begin{equation} \mathcal{H} = \frac{A}{3 {\rm \pi}\gamma a^2} \end{equation}
represents the square of the molecular length scale relative to the drop radius. The numerical factor 
$3{\rm \pi}$ is included for later convenience. Typically, 
$\mathcal {H}$ is very small, 
$\mathcal {H} \ll 1$. For example, the silicone-oil drops used in the post-contact coalescence experiments of Paulsen et al. (Reference Paulsen, Burton, Nagel, Appathurai, Harris and Basaran2012) have 
$A = 10^{-20}\ \mbox {J}$, 
$\gamma = 0.02\ \mbox {N} \mbox {m}^{-1}$ and 
$a = 0.001\ \mbox {m}$, which gives 
$\mathcal {H} \approx 5 \times 10^{-14}$.
For 
$\mathcal {H}\ll 1$, the gap-width scale is much shorter than the radial length scale, with 
$h/r \sim \mathcal {H}^{1/6}\ll 1$. This suggests that the deformed surface profile has a shallow gradient, 
$\partial h / \partial t \ll 1$, as is indeed confirmed by a full calculation of the equilibrium profile (e.g. see Quinn et al. Reference Quinn, Feng and Stone2013). We therefore make a shallow-slope approximation that considerably simplifies the disjoining pressure and the dynamics of jump-to-contact, as we will see in the next subsection.
2.1. Shallow-slope approximation
For a drop with a shallow surface gradient, 
$\partial h / \partial r \ll 1$, the van der Waals disjoining pressure is given at leading order by the Derjaguin approximation (Derjaguin Reference Derjaguin1955)
\begin{equation} {\varPi}={-} \frac{A}{6 {\rm \pi}(2{h})^3}. \end{equation}
In this approximation, the disjoining pressure at some point 
$r$ is given by the disjoining pressure between two infinite, flat, parallel surfaces with separation 
$2 h( r, t)$. The Derjaguin approximation has been shown to give excellent agreement with full body-force calculations in many problems similar to that considered in this paper, provided that the shallow slope is maintained (e.g. Jiang & James Reference Jiang and James2006). Furthermore, under the shallow-slope approximation, the curvature can be linearized to
\begin{equation} \kappa = \frac{1}{ r}\,\frac{\partial}{\partial r} r\,\frac{\partial h }{\partial r}. \end{equation} The shallow-slope approximation also allows the geometry of the flow problem to be simplified. As the flow in the drop is driven by gradients of the normal stress on the surface, the characteristic flow length scale is given by the scale on which the surface stress varies, that is, the radial scale 
$r\gg h$ of the near-contact region. Thus the gap width 
$2h$ is small compared to the flow length scale 
$r$, so the flow domain can be linearized to 
$z > 0$, with boundary 
$z = 0$ and outward surface normal 
$\boldsymbol {n} = - \boldsymbol {e}_z$. The flow problem reduces to the flow in a half-space due to an applied normal stress on the surface. The evolution of the gap width 
$2h$ influences the flow only though the dynamic boundary condition, with both the surface tension and van der Waals disjoining pressure depending on it through (2.6) and (2.7).
It can also be shown that the nonlinear terms in (2.2a,b) and (2.3) can be neglected. As 
$\partial h / \partial r \ll 1$, the consistent leading-order balance in the kinematic boundary condition (2.3) is 
$\partial h / \partial t = {u}_z$, so the time scale for evolution of the gap width and contact is 
$t \sim h / u$. The relative sizes of the inertial terms in the Navier–Stokes equations can then be compared:
\begin{equation} \frac{\boldsymbol{u} \boldsymbol{\cdot}\boldsymbol{\nabla} \boldsymbol{u}}{\partial \boldsymbol{u} / \partial t} \sim \frac{u^2/ r}{u/(h/u)} \sim \frac{ h}{ r }\ll 1. \end{equation}
Consequently, the nonlinear term can be neglected in (2.2a,b) to give the unsteady Stokes equations in 
$z>0$:
\begin{equation} \rho\,\frac{\partial \boldsymbol{u}}{\partial t} ={-}\boldsymbol{\nabla} p + \mu\,\nabla^2 \boldsymbol{u},\quad \boldsymbol{\nabla} \boldsymbol{\cdot} \boldsymbol{u}=0. \end{equation} We non-dimensionalize the linearized problem according to the macroscopic scale of the drop by scaling lengths with the drop radius 
$a$, velocities with a capillary speed 
$\gamma /2 \mu$, pressures with a capillary pressure 
$\gamma /2 a$, and time with 
$a \mu / \gamma$. The numerical constants are included for convenience. As the deformation due to van der Waals attraction is localized to the small near-contact region of radius 
$r \ll a$, we define perturbations from the far field by subtracting the spherical far-field profile (
${\sim }r^2/2a$ for 
$r \ll a$) and the associated uniform capillary pressure 
$2 \gamma / a$ to obtain a local perturbation gap width 
$d = 2 h - r^2/a$ and perturbation pressure 
$p-2 \gamma / a$. Thus, using tildes temporarily to denote dimensionless variables, we define
\begin{equation} \tilde{d} = \frac{1}{a }\left(2 h - \frac{ r^2}{a}\right),\quad (\tilde{r},\tilde{z}) = \frac{(r,z)}{a }, \quad \tilde{\boldsymbol{u}} = \frac{2 \mu\boldsymbol{u}}{\gamma }, \quad \tilde{p} = \frac{2 a }{\gamma}\left( p - \frac{2 \gamma}{a}\right), \quad \tilde{t} = \frac{\gamma t }{a \mu}. \end{equation}
The dimensionless governing equations (2.9a,b) in 
$\tilde {z}>0$ become
\begin{equation} \frac{1}{Oh^2}\,\frac{\partial {\tilde{\boldsymbol{u}}}}{\partial \tilde{t}} ={-}\tilde{\boldsymbol{\nabla}}\tilde{p} + \tilde{\nabla}^2 \tilde{\boldsymbol{u}},\quad \tilde{\boldsymbol{\nabla}}\boldsymbol{\cdot}\tilde{\boldsymbol{u}}=0, \end{equation}
where the Ohnesorge number 
$Oh \equiv \mu / (\rho a \gamma )^{1/2}$ relates viscous, inertial and surface tension forces on the scale of the drop. The dynamic boundary condition (2.1) becomes
\begin{equation} \tilde{\boldsymbol{\sigma}} \boldsymbol{\cdot}\boldsymbol{e}_z ={-} \left( \frac{1}{ \tilde{r}}\,\frac{\partial}{\partial \tilde{r}} \tilde{r}\,\frac{\partial \tilde{d}}{\partial \tilde{r}} -\frac{\mathcal{H}}{ (\tilde{d}+\tilde{r}^2)^3}\right)\boldsymbol {e_z} \quad \mbox{on} \ \tilde{z} =0 . \end{equation}From now on, we drop the tildes. The coupled evolution of the surface is given by the dimensionless kinematic condition
\begin{equation} \frac{\partial d}{\partial t} = u_z \quad \mbox{on}\ z =0. \end{equation}The flow in the drop given by (2.11a,b) includes both viscosity and inertia. The relative strength of the viscous and inertial terms, determined by the scales of the flow in the near-contact region, sets the leading-order force balance for the flow in the drop. In the next section, we investigate the different dynamical regimes for jump-to-contact that arise depending on the dominant force balance.
3. Dynamical regimes
To compare the relative strength of the viscous and inertial forces during jump-to-contact, we define a Reynolds number 
$\mathcal {R}$ from the ratio of the inertial and viscous terms in the unsteady Stokes equation (2.11a,b):
\begin{equation} \mathcal{R} \sim \frac{1}{Oh^2}\,\frac{ \partial \boldsymbol{u} / \partial t}{\nabla^2 \boldsymbol{u}} \sim \frac{1}{Oh^2}\,\frac{u/t}{ u/r^2} \sim \frac{1 }{Oh^2 }\,\frac{u r^2 }{d } , \end{equation}
where, as before, 
$t \sim d /u$ from the surface evolution, and 
$Oh^2$ has taken the place of the kinematic viscosity 
$\nu$ in the dimensionless equations. This Reynolds number is greater by a factor 
$r/d\gg 1$ than the usual Reynolds number, 
$u r / Oh^2$, as the time scale 
${\sim }d/u$ of the gap evolution is much less than the time scale 
${\sim } r/u$ of the flow in the drop. We term 
$\mathcal {R}$ the ‘increased Reynolds number’, by analogy with the ‘reduced Reynolds number’ in lubrication theory. To estimate the velocity scale in (3.1), we assume that it is given by a viscous balance in (2.11a,b) between 
$\nabla ^2 \boldsymbol {u}$ and the pressure gradient 
$\boldsymbol {\nabla } p$. As van der Waals attraction drives jump-to-contact, it must be included at leading order in any dynamical regime, thus the pressure scales with the disjoining pressure and 
$\boldsymbol {\nabla } p \sim \mathcal {H}/rd^3$. Therefore, the viscous balance gives 
$u \sim \mathcal {H} r / d^3$, and the increased Reynolds number scales as
\begin{equation} \mathcal{R}\sim \frac{\mathcal{H}}{Oh^2 }\,\frac{ r^3 }{d^4 }. \end{equation}
The characteristic gap width 
$d$ decreases during jump-to-contact. We therefore expect the other scales of the flow, and consequently 
$\mathcal {R}$, also to evolve.
From the increased Reynolds number, the dynamics can be classified into three distinct possible regimes based on the current gap width and radial scale of the near-contact region, and the fluid properties of the drop. When 
$\mathcal {R} \ll 1$, the viscous forces are much greater than the inertial forces, suggesting a Stokes regime where the flow in the drop is characterized by viscous resistance to the van der Waals attraction. Conversely, when 
$\mathcal {R}\gg 1$, there is a possible inertial regime where viscosity is negligible. Finally, when 
$\mathcal {R}\sim 1$, both inertia and viscosity act at leading order. In general, 
$\mathcal {R}\sim 1$ could occur either as a transition between the Stokes and inertial regimes, or as a distinct inertial–viscous regime with a stable balance between the viscous, inertial and van der Waals forces. We assess this possibility using scaling arguments in § 3.2.
At the onset of jump-to-contact, the gap width and radial scale are, respectively, 
$d \sim \mathcal {H}^{1/3}$ and 
$r \sim \mathcal {H}^{1/6}$, from (2.4a,b). Thus initially,
\begin{equation} \mathcal{R}\sim \frac{\mathcal{H}^{1/6}}{Oh^2}. \end{equation}
As discussed previously, 
$\mathcal {H} \ll 1$ for the millimetre-scale fluid drops used typically in jump-to-contact and drop-coalescence experiments. However, the Ohnesorge number 
$Oh$ and hence 
$\mathcal {R}$ can vary considerably. For example, the silicone-oil drops used in Paulsen et al. (Reference Paulsen, Burton, Nagel, Appathurai, Harris and Basaran2012) had a range of viscosities such that 
$0.1 < Oh^2< 1000$ and initially 
$10^{-5} < \mathcal {R}<0.05$; thus the jump-to-contact dynamics began in the Stokes regime. However, for less viscous drops, inertia is significant from the outset: for example, millimetre drops of water have 
$\mathcal {R}\sim 450$ initially.
In the following subsections, we first consider the Stokes regime for jump-to-contact, before returning to the effects of inertia.
3.1. Stokes regime
In our previous paper (BL22), we examined the Stokes regime with 
$\mathcal {R} \equiv 0$. As the van der Waals attraction diverges as contact is approached, jump-to-contact is an example of a flow with a finite-time singularity, hence we might expect self-similar behaviour (see Eggers & Fontelos Reference Eggers and Fontelos2005). We searched for possible self-similarity with a power-law dependence on the time until contact 
$t'$ by positing
\begin{equation} r,z \sim t{'^{\alpha}}, \quad d \sim t{'^{\beta}}, \quad u \sim t{'^{\gamma}}, \end{equation}
where the exponents 
$\alpha$, 
$\beta$ and 
$\gamma$ are to be found. The gap-width evolution equation (2.13) determines 
$\gamma = \beta - 1$. Then gradients of the van der Waals disjoining pressure scale with time as 
$\boldsymbol {\nabla } p \sim 1 / r d^3 \sim t{'^{-\alpha -3 \beta }}$, and viscous diffusion scales as 
$\nabla ^2 u \sim u/r^2 \sim t{'^{\beta - 1 - 2 \alpha }}$. Force balance for Stokes flow requires that the time exponents of these two terms match, determining 
$\beta$ in terms of 
$\alpha$ as 
$\beta = (1+\alpha )/4$. To fix both exponents, we would require another term in the leading-order force balance. The only other force in the Stokes regime (inertia being ruled out) is from surface tension. However, as jump-to-contact begins when van der Waals attraction overcomes surface tension, it seems unlikely that surface tension should act at leading order in any similarity solution. We thus seem to have a second-kind similarity problem where the force balance does not fix the time exponent, which instead must be selected by other means, for example by an existence condition (cf. capillary pinch-off of viscous threads; Papageorgiou Reference Papageorgiou1995).
In BL22, we performed full time-dependent simulations of jump-to-contact in the Stokes regime, including surface tension. Figure 2 shows the evolution of the gap-width profile 
$d(r,t)$ at fixed values of the minimum gap width 
$d_0 \equiv d(0,t)$. Contact is made in finite time, consistent with the definition of 
$t'$. Figure 3 shows how the maximum surface velocity 
$\dot {d\hskip 1.4pt}\hskip -1.4pt_0 \equiv \mathrm {d} d_0 / \mathrm {d} t$ and a radial scale 
$r^*$, defined by 
$\dot {d\hskip 1.4pt}\hskip -1.4pt(r^*,t) = \dot {d\hskip 1.4pt}\hskip -1.4pt_0(t)/2$, vary with the minimum gap width 
$d_0$ in log–log space. After an initial transient, 
$\ln |\dot {d\hskip 1.4pt}\hskip -1.4pt_0|$ varies linearly with 
$-{\ln d_0}$. The observed gradient is 
$2$, suggesting that 
$\dot {d\hskip 1.4pt}\hskip -1.4pt _0 \propto d_0^{-2}$, which integrates to give 
$d_0 \propto t{'^{1/3}}$ as contact is approached. Similarly, 
$-{\ln r^*}$ varies linearly with gradient 
$1$, therefore 
$r^* \propto d_0 \propto t{'^{1/3}}$. The clear power-law behaviour is consistent with the self-similar scalings (3.4a–c) for the value 
$\alpha = 1/3$.
Figure 2. Gap-width profile 
$d(r,{t})$ in the Stokes regime scaled by the initial gap-width and radial scales (2.4a,b), and shown for minimum gap widths 
$d_0(t)/\mathcal {H}^{1/3} = 0.9$, 
$0.7$, 
$0.5$, 
$0.3$, 
$0.1$ and 
$10^{-8}$ corresponding to times 
$t = 0$, 
$2.552$, 
$3.372$, 
$3.569$, 
$3.604$ and 
$3.605$, respectively.
Figure 3. The evolution in the Stokes regime of 
$\ln |\dot {d\hskip 1.4pt}\hskip -1.4pt_0|$ (solid line) and 
$- \ln r^*$ (dashed line) plotted against 
$-{\ln d_0({t})}$. After the initial transient, the solid line has gradient 
$2$, and the dashed line has gradient 
$1$.
With 
$\alpha = 1/3$, the curvature scales as 
$\kappa \sim t{'^{-1/3}}$. Thus the surface tension is indeed much weaker than the disjoining pressure 
$\varPi \sim t{'^{-1}}$ as contact is approached, and as expected, is therefore not part of the leading-order force balance. Neglecting surface tension, we found an analytic solution for the self-similar Stokes flow and gap-width profile with 
$\alpha = 1/3$. Remarkably, the similarity profile is an exact hyperbola:
\begin{equation} d(r,t') = (d_0(t')^2 + \theta^2 r^2)^{1/2}, \quad d_0(t') = \left(\frac{3 \mathcal{H}\,t'}{2 \theta}\right)^{1/3}, \end{equation}
where 
$\theta$ is a constant defined by 
$\lim _{r \to \infty }\partial d / \partial r \to \theta$. The value of 
$\theta$ in (3.5a,b) is arbitrary; it is determined by the specific choice of initial conditions. Whilst the exact processes that set 
$\theta$ were not investigated fully, the time-dependent simulations demonstrated that the initial curvature at the origin strongly influenced the observed value of 
$\theta$, and hence the aspect ratio of the similarity solution. As there is some memory of initial conditions in the similarity solution, viscous jump-to-contact is an example of non-universal self-similarity.
We argued in BL22 that the value 
$\alpha = 1/3$ is selected by the regularity condition of finite curvature at 
$r=0$, which physically is enforced by surface tension, even though surface tension is not part of the leading-order force balance. (A loose analogy might be drawn with finger-width selection for vanishing surface tension in the Saffman–Taylor problem.) Selection of 
$\alpha = 1/3$ by regularity at 
$r=0$ is supported first by the fact that solutions to the similarity equations for 
$\alpha \ne 1/3$ are non-smooth at the origin, for example as shown in figure 4 for 
$\alpha = 0.44$. Second, in time-dependent simulations with no surface tension and non-smooth initial conditions, the gap profile converges to one of these non-smooth similarity solutions. Finally, in recent unpublished time-dependent simulations with surface tension and a generalized disjoining pressure 
$\varPi _q \propto -d^{-q}$, we have found convergence to smooth self-similar solutions with irrational exponents 
$\alpha$ and 
$\beta = (1+\alpha )/(1+q)$, where 
$\alpha (q)$ is determined by regularity at 
$r=0$, and 
$\alpha \ne \beta$ for 
$q \ne 3$.
Figure 4. Plot of two solutions to the similarity equations (equation (9) in BL22) for the gap profile in the Stokes regime with different 
$\alpha$. The profile for 
$\alpha = 0.44$ is not smooth at 
$r=0$, but the profile for the selected value 
$\alpha = 1/3$ is.
We note that the similarity solution (3.5a,b) has the same time scaling for both the gap-width and radial scales, thus 
$d/r$ remains of constant magnitude. Thus as the shallow-slope approximation 
$d/r\ll 1$ is initially valid for 
$\mathcal {H} \ll 1$, it remains valid throughout jump-to-contact in the Stokes regime.
An important inference from BL22 for the purposes of this paper is that in the Stokes regime with self-similar behaviour (3.5a,b), the increased Reynolds number evolves as
\begin{equation} \mathcal{R}\sim \frac{\mathcal{H}^{1/6}}{Oh^2}\,t{'^{{-}1/3}}. \end{equation}
Thus even if viscosity is initially dominant (
$Oh^2 \gg \mathcal {H}^{1/6}$), there is a crossover time 
$t'_c \sim \mathcal {H}^{1/2}/Oh^6$ at which 
$\mathcal {R}$ is order 
$1$ and inertia becomes significant. At crossover, the minimum gap width scales as 
$d_0 \sim \mathcal {H}^{1/3} t_c{'^{1/3}} \sim \mathcal {H}^{1/2}/Oh^2$, which is potentially very small. From (2.5) and (2.10a–e), 
$\mathcal {H}^{1/2}$ is the non-dimensionalized molecular length scale for the drop, thus the continuum approximation fails when 
$d_0 \sim \mathcal {H}^{1/2}$. For the purposes of continuum-scale modelling, inertia is significant only if the crossover gap width is reached before the continuum approximation fails, 
$\mathcal {H}^{1/2}/Oh^2 \gg \mathcal {H}^{1/2}$. We conclude that for moderately viscous drops with 
$Oh^2 \gtrsim 1$, inertia is negligible throughout the jump-to-contact dynamics, but for less viscous drops with 
$Oh^2 \ll 1$, inertia becomes significant before contact is reached.
3.2. Inertial–viscous regime
We now consider the possibility of a stable balance between inertial, viscous and van der Waals forces for 
$\mathcal {R} \sim 1$. Setting 
$\mathcal {R} \sim 1$ in (3.2) gives a scaling relation for the radial length scale in terms of the gap width:
\begin{equation} \frac{ \mathcal{H}}{Oh^2}\,\frac{r^3}{ d^4} \sim 1 \implies r \sim \frac{Oh^{2/3} }{\mathcal{H}^{1/3}}\,d^{4/3}. \end{equation}
As before, the velocity scales as 
$u\sim \mathcal {H} r / d^3$. With the radial length scale (3.7), the velocity scales with the gap width as 
$u \sim Oh^{2/3} \, \mathcal {H}^{2/3} /{d^{5/3} }$. Then the kinematic condition (2.13) can be integrated,
\begin{equation} \frac{\partial d}{\partial t} \sim { Oh^{2/3} \, \mathcal{H}^{2/3} }{d^{{-}5/3} } \implies d \sim { Oh^{1/4} \, \mathcal{H}^{1/4} }t'^{3/8}, \end{equation}
to give the dependence of the gap width on the time until contact 
$t'$. We note that the three-way force balance in an inertial–viscous regime fixes both the time exponent and the parametric dependence to give a first-kind similarity problem. Substituting the scaling (3.8) for the gap width back into the radial and velocity scales gives
\begin{equation} d \sim Oh^{1/4} \, \mathcal{H}^{1/4} t'^{3/8}, \quad r\sim Oh \, t'^{1/2},\quad u \sim Oh^{1/4} \, \mathcal{H}^{1/4} t'^{{-}5/8}. \end{equation} When determining (3.9a–c), we did not include the surface tension in the leading-order force balance. From (3.9a–c), the curvature increases as 
$\kappa \sim d / r^2 \sim Oh^{-7/4} \, \mathcal {H}^{1/4} t'^{-5/8}$, and the disjoining pressure increases as 
$\varPi \sim \mathcal {H}/d^3 \sim Oh^{-3/4} \,\mathcal {H}^{1/4} t'^{-9/8}$. We see that the growth of the disjoining pressure is stronger than the growth of the surface tension as 
$t' \to 0$. Therefore, as surface tension is initially no greater than the disjoining pressure, it can be neglected.
We also consider the validity of the shallow-slope approximation in an inertial–viscous regime. With scalings (3.9a–c), the surface gradient is 
$\partial d/\partial r \sim (\mathcal {H}/Oh^{2}\, d)^{1/3}$. Thus the shallow-slope approximation 
$\partial d/\partial r \ll 1$ remains valid until the gap width 
$d$ shrinks to 
$d_s \sim \mathcal {H}/Oh^2$. We compare 
$d_s$ to the molecular length scale 
$\mathcal {H}^{1/2}$. For 
$Oh^2 \gtrsim \mathcal {H}^{1/2}$, the critical gap width 
$d_s \lesssim \mathcal {H}^{1/2}$ is therefore not reached on continuum scales, so the shallow-slope approximation remains valid throughout jump-to-contact. We note that 
$\mathcal {H}^{1/2} \ll 1$ is typically very small, for example 
$\mathcal {H}^{1/2}\approx 2 \times 10^{-7}$ for millimetre-scale drops, so the condition 
$Oh^2 \gtrsim \mathcal {H}^{1/2}$ is satisfied by drops with even fairly low viscosity, including millimetre drops of water that have 
$Oh^2 \approx 10^{-5}$.
The arguments within this subsection have identified the possible scaling behaviour of an inertial–viscous regime and demonstrated that it is consistent with the shallow-slope approximation. We have not yet demonstrated such a regime to exist or be stable. In § 4, we describe numerical simulations to investigate the inertial–viscous regime further.
3.3. Inertial regime
For 
$\mathcal {R} \gg 1$, viscosity is negligible at leading order and there is a two-way force balance between inertia and van der Waals attraction. We again posit time-dependent scalings of the form (3.4a–c) with 
$\gamma = \beta - 1$ from the kinematic boundary condition (2.13). The inertial and pressure terms in (2.11a,b) scale as 
$\partial \boldsymbol {u}/ \partial t \sim t{'^{\beta - 2}}$ and 
$\boldsymbol {\nabla }p \sim t{'^{-3 \beta -\alpha }}$, thus equating time exponents gives 
$\beta = (2-\alpha )/4$ and consequently
\begin{equation} r,z \sim t{'^{\alpha}}, \quad d \sim t{'^{(2-\alpha)/4}}, \quad u \sim t{'^{-(2+\alpha)/4}}. \end{equation}
As in the Stokes regime, the time-dependent scalings in the inertial regime are under-determined with 
$\alpha$ not being fixed by the force balance. We defer further discussion of a possible inertial regime to § 4.3.
4. Numerical simulations of jump-to-contact with inertia
4.1. Method
We now consider time-dependent simulations of the full unsteady Stokes equations (2.11a,b)–(2.13). Motivated by the possibility of an inertial–viscous balance, we scale the dimensionless parameters 
$\mathcal {H}$ and 
$Oh$ out of the problem according to (3.9a–c), and define rescaled variables
\begin{equation} \bar{d} = \frac{d}{Oh^{1/4} \, \mathcal{H}^{1/4}},\quad (\bar{r},\bar{z}) = \frac{(r,z)}{Oh},\quad \bar{\boldsymbol{u}} = \frac{u}{Oh^{1/4} \, \mathcal{H}^{1/4}}, \quad \bar{p} = \frac{Oh^{3/4}}{\mathcal{H}^{1/4}}\,p. \end{equation}
The unsteady Stokes equations (2.11a,b) in 
$\bar {z}>0$ become
\begin{equation} \frac{\partial \bar{\boldsymbol{u}}}{\partial \bar{t}} ={-} \bar{\boldsymbol{\nabla}} \bar{p} + \bar{\nabla}^2 \bar{\boldsymbol{u}},\quad \bar{\boldsymbol{\nabla}}\boldsymbol{\cdot}\bar{\boldsymbol{u}}=0. \end{equation}As surface tension is negligible in an inertial–viscous regime, the dynamic boundary condition (2.12) is given by the disjoining pressure at leading order:
\begin{equation} \bar{\boldsymbol{\sigma}}\boldsymbol{\cdot}\boldsymbol{e}_z =\frac{1}{ \bar{d}^3}\,\boldsymbol{e}_z \quad \mbox{on}\ \bar{z} =0 . \end{equation}From now on, we drop the bars on the rescaled variables.
We solve (4.2a,b) as an evolution equation for the velocity due to the instantaneous pressure gradient 
$\boldsymbol {\nabla } p$ and viscous diffusion 
$\nabla ^2 \boldsymbol {u}$. As the fluid is incompressible, 
$\boldsymbol {\nabla }\boldsymbol {\cdot }\boldsymbol {u} =0$, the velocity field can be expressed in terms of a Stokes streamfunction 
$\varPsi$ with 
$\boldsymbol {u} = \boldsymbol {\nabla }\times (\varPsi \boldsymbol {e}_{\theta }/r)$. Furthermore, taking the divergence of (4.2a,b) shows that 
$\nabla ^2 p =0$. Thus 
$\boldsymbol {\nabla } p$ is also divergence-free and so can also be expressed in terms of a vector potential. As the problem is axisymmetric, this vector potential can be expressed in terms of a scalar function 
$\varPhi$ with 
$\boldsymbol {\nabla } p = \boldsymbol {\nabla } \times (\varPhi \boldsymbol {e}_{\theta }/r)$. Using the functions 
$\varPsi$ and 
$\varPhi$, the unsteady Stokes momentum equation (4.2a,b) can be uncurled to give
\begin{equation} \frac{\partial \varPsi}{\partial t} ={-} \varPhi + E^2 \varPsi + c_1, \end{equation}
where 
$E^2 \equiv \partial ^2/\partial r^2 - r^{-1}\,\partial /\partial r + \partial ^2 / \partial z^2$, and 
$c_1$ is an arbitrary constant that gives no flow and so can be set to 
$0$ without loss of generality.
We view (4.4) as a forced diffusion equation for the Stokes streamfunction 
$\varPsi$, and hence the velocity. The viscous diffusion term 
$E^2 \varPsi$ is given instantaneously by the current streamfunction field. The instantaneous forcing 
$\varPhi$ can be determined from the current streamfunction field and the gap profile by the dynamic boundary condition as follows. By the substitution of 
$\boldsymbol {\sigma } = - p \boldsymbol{\mathsf{I}} + \boldsymbol {\nabla }\boldsymbol {u} +(\boldsymbol {\nabla } \boldsymbol {u})^\mathrm {T}$ into (4.3), where 
${\boldsymbol{\mathsf{I}}}$ is the identity matrix, we obtain
\begin{equation} p = 2\,\frac{\partial u_z}{\partial z} - \frac{1}{ d^3}\quad \mbox{on}\ z =0 . \end{equation}
The Stokes streamfunction 
$\varPsi$ determines the viscous normal stress 
$2\,\partial u _z / \partial z$, hence the pressure 
$p$ is known on 
$z=0$. As the pressure is harmonic (
$\nabla ^2 p = 0$), it is determined throughout 
$z>0$ by its value on 
$z=0$. Thus the pressure is known everywhere, thence the vector potential for 
$\boldsymbol {\nabla }p$, and hence 
$\varPhi$, can be found everywhere. We present further details and an explicit boundary integral representation of 
$\varPhi$ in Appendix A.
The normal component of the dynamic boundary condition (4.3) is enforced when determining 
$\varPhi$. However, the tangential component remains and is applied as a boundary condition to (4.4). In terms of the Stokes streamfunction, this yields
\begin{equation} \frac{\partial }{\partial r}\frac{1}{r}\,\frac{\partial \varPsi}{\partial r}-\frac{1}{r}\,\frac{\partial^2 \varPsi}{\partial z^2}= 0 \quad \mbox{on}\ {z} =0. \end{equation} The gap-width evolution, and hence the evolution of the disjoining pressure, is described by the kinematic equation (2.13), given in terms of the Stokes streamfunction 
$\varPsi$ as
\begin{equation} \frac{\partial d}{\partial t} = \frac{1}{r}\,\frac{\partial \varPsi}{\partial r} \quad \mbox{on}\ z =0. \end{equation}
Equations (4.4)–(4.7), together with those in Appendix A, are sufficient to determine the evolution of 
$p$, 
$\varPhi$, 
$\varPsi$ and 
$d$.
We performed numerical simulations of (4.4)–(4.7) on an adaptive grid with tighter spacing of grid points in the centre and increasing spacing with distance, chosen both to resolve the radial length scale of the near-contact region and to capture the far-field behaviour. At a given time, the forcing term 
$\varPhi$ in (4.4) is calculated at each grid point by the integrals (A8)–(A10) using the current values of the Stokes streamfunction 
$\varPsi$ and gap-width profile 
$d$. The viscous term 
$E^2 \varPsi$ is calculated using central differencing and the boundary condition (4.6). Given 
$\varPhi$ and 
$E^2 \varPsi$, the Stokes streamfunction 
$\varPsi$ is advanced on the grid by time stepping (4.4) using the forward Euler method. The gap width 
$d$ is advanced likewise by (4.7). Throughout the dynamics, the grid is rescaled adaptively to maintain a good resolution of the evolving length scale of the flow, with 
$\varPsi$ and 
$d$ interpolated at new grid points using a quintic spline. Further details of the numerical scheme are given in Appendix B. Also included in the appendix is a verification of the numerical parameters that we used in the simulations.
4.2. Results
We took initial conditions with no flow in the drops and prescribed gap-width profiles, with results from a wide range of initial profiles presented below. We note that the initial size of the increased Reynolds number 
$\mathcal {R}$ is determined by the gap width 
$d$ and radial length scale 
$r$ of the initial profile by (3.2), which becomes 
$\mathcal {R} \sim r^3/d^4$ in the rescaled variables defined by (4.1a–d).
We first considered evolution from an initial gap-width profile 
$d(r,0) = (1+0.4r^2)^{1/2}$. This choice has increased Reynolds number 
$\mathcal {R} \sim 4$ and thus dynamics that are initially in the inertial–viscous regime. The evolution of the gap-width profile is shown at selected times in figure 5, with contact made in finite time. To examine the predicted self-similar behaviour, we show in figure 6 the evolution of 
$\ln |\dot {d\hskip 1.4pt}\hskip -1.4pt_0|$ with 
$-{\ln d_0}$ (a proxy for time until contact 
$t'$). After the initial transient, there is a clear linear relationship between the two quantities, with gradient 
$5/3$. This relationship implies that 
$\dot {d\hskip 1.4pt}\hskip -1.4pt_0 \propto d_0^{-5/3}$, which can be integrated to give 
$d_0 \propto t{'^{3/8}}$ and 
$\dot {d\hskip 1.4pt}\hskip -1.4pt_0 \propto t{'^{-5/8}}$. Thus the observed time dependence of both the velocity and the gap width matches the scaling behaviour (3.9a–c) expected in an inertial–viscous regime. To examine the evolution of the radial length scale of the flow, we define 
$r^*$ to be the radius at which the surface velocity is half its maximum value, 
$\dot {d\hskip 1.4pt}\hskip -1.4pt(r^*,t) = \dot {d\hskip 1.4pt}\hskip -1.4pt_0(t)/2$. The observed surface velocity is peaked at 
$r=0$ (e.g. see figure 10), which suggests that 
$r^*$ gives a characteristic radial scale for the near-contact region. We also show in figure 6 the evolution of 
$-{\ln r^*}$ with 
$-{\ln d_0}$. After the initial transient, there is again a linear relationship with a different gradient, 
$4/3$, which implies 
$r^*\propto d_0^{4/3}\propto t{'^{1/2}}$. Thus the observed radial length scale also agrees with the expected scaling behaviour (3.9a–c) in the inertial–viscous regime.
Figure 5. Evolution with inertia under (4.4)–(4.7) of the initial gap-width profile 
$d(r,0)=(1+0.4r^2)^{1/2}$ with inertia, shown (from top to bottom) at minimum gap widths 
$d_0(t)=1$, 
$0.7$, 
$0.3$, 
$0.2$, 
$0.1$ and 
$10^{-8}$. The corresponding times are 
$t=0$, 
$1.085$, 
$1.464$, 
$1.490$, 
$1.502$ and 
$1.504$.
Figure 6. The evolution of 
$\ln |\dot {d\hskip 1.4pt}\hskip -1.4pt_0|$ (solid line) and 
$-{\ln r^*}$ (dashed line) with 
$-{\ln d_0(t)}$ for the initial gap-width profile 
$d(r,0)=(1+0.4r^2)^{1/2}$. The characteristic radial length scale 
$r^*$ is defined by 
$\dot {d\hskip 1.4pt}\hskip -1.4pt(r^*,t) = \dot {d\hskip 1.4pt}\hskip -1.4pt_0(t)/2$. After the initial transient, the solid line tends to a constant gradient 
$5/3$ and the dashed line tends to a constant gradient 
$4/3$, consistent with the predicted inertial–viscous regime.
Having established the adherence of the numerical simulations to the scaling relations (3.9a–c), we solved directly for the corresponding self-similar solutions in similarity variables 
$d/d_0$, 
$r/d_0^{4/3}$ and 
$u/d_0^{-5/3}$. We used an iterative procedure, starting from an initial guess, in which each iteration consisted of taking a fixed number of steps (
$1000$) with the time-dependent code on a fixed grid, followed by a rescaling of all variables with the relevant power of 
$d_0$. We continued to iterate until the gap width 
$d/d_0$ at radius 
$r/d_0^{4/3}=40$ did not change by more than 
$10^{-6}$ in one iteration, with 
$d_0$ more than halving in this time. The self-similar gap-width profile, shown in figure 7 (dashed line), is quadratic near 
$r=0$ with 
$d/d_0 \approx 1+0.076 r^2/d_0^{8/3}$, and has 
$d /d_0\approx 0.57 r^{3/4}/d_0$ in the far field. Figure 7 also shows the evolution of the gap-width profile in similarity variables, from time-dependent simulations of initial profiles 
$d(r,0) = (1+0.4 r^2)^{1/2}$ and 
$d(r,0) = (1+0.04 r^2)^{1/2}$. There is clear convergence to the same self-similar profile from both sides, confirming that the dynamics of jump-to-contact in the inertial–viscous regime is self-similar as contact is approached. Convergence to the same solution from both these initial conditions suggests that it is stable and universal.
Figure 7. The evolution of the rescaled gap-width profile 
$d(r,t)/d_0(t)$ plotted against the rescaled radial variable 
$r/d_0(t)^{4/3}$ for initial profiles (a) 
$d(r,0) = (1+0.4 r^2)^{1/2}$, and (b) 
$d(r,0) = (1+0.04 r^2)^{1/2}$. The profiles are shown (solid lines) at minimum separations (a) 
$d_0(t) = 1$, 
$0.03$, 
$0.003$, 
$10^{-4}$, 
$10^{-6}$, 
$10^{-9}$, and (b) 
$d_0(t) = 1$, 
$0.2$, 
$0.005$, 
$10^{-4}$, 
$10^{-6}$, 
$10^{-9}$. In both cases, the profiles converge to the same self-similar solution (dashed lines).
To investigate further the universality of this solution, we considered the evolution of the rescaled characteristic radius 
$r^*/d_0(t)^{4/3}$, which must converge for any similarity solution, for a range of initial profiles. Figure 8 shows that in each case, 
$r^*/d_0(t)^{4/3}$ converges to the same value 
${\approx }2.94$ observed in the self-similar solution, providing further strong evidence that the self-similar solution shown in figure 7 is universal. We note that the initial profiles considered have a wide range of radial length scales, corresponding to several orders of magnitude of initial reduced Reynolds numbers, 
$0.1 \lesssim \mathcal {R} \lesssim 4000$. Figure 8 also shows that the manner in which 
$r^*/d_0(t)^{4/3}$ converges is very similar in simulations with initial profiles 
$d(r,0)=(1+0.4r^2)^{1/2}$ and 
$d(r,0) = 1+0.2r^2$, which have the same quadratic behaviour, and hence radial scale, at 
$r=0$, but are linear and quadratic in the far field, respectively. This suggests that the dynamics is determined primarily by the surface profile near 
$r=0$.
Figure 8. Evolution of the parameter 
$r^{*}/d_0(t)^{4/3}$ with 
$- {\ln d_0(t)}$ for various initial gap profiles. In each case, there is convergence to 
$r^{*}/d_0(t)^{4/3} \approx 2.94$ (dashed line) corresponding to the similarity profile shown in figure 7. (a) Initial profiles 
$d(r,0) = (1+0.4 r^2)^{1/2}$ (solid line) and 
$d(r,0) = 1+0.2 r^2$ (dot-dashed line) have the same curvature near 
$r=0$ but are linear and quadratic in the far field, respectively. In these cases, 
$r^{*}/d_0(t)^{4/3}$ converges in a very similar manner. (b–d) Initial profiles 
$d(r,0) = (1+4 r^2)^{1/2}$, 
$(1+0.04 r^2)^{1/2}$ and 
$(1+0.004r^2)^{1/2}$, respectively, show the effect of increasing the initial increased Reynolds number 
$\mathcal {R}$.
Furthermore, we can show that the decay to the similarity solution shown in figure 8 is algebraic with a constant exponent by plotting the evolution of 
$\ln (r^*/d_0^{4/3}-2.94 )$ against 
$-{\ln d_0(t)}$ in figure 9. After the initial transient, there is a linear relationship with gradient approximately 
$-1/5$. Thus perturbations to the self-similar solution decay as 
$d_0^{1/5}\propto t{'^{3/40}}$ as contact is approached. We do not have a theoretical argument for 
$-1/5$, and it may be that 
${\approx }-0.2$ is just the least stable eigenvalue in a system of stable modes (for a similar problem, cf. Sierou & Lister Reference Sierou and Lister2003).
Figure 9. Evolution of 
$\ln ( r^{*}/d_0^{4/3}-2.94)$ with 
$- {\ln d_0(t)}$ for initial gap profiles (a) 
$d(r,0) = (1+0.4 r^2)^{1/2}$, and (b) 
$d(r,0) = (1+0.04 r^2 )^{1/2}$. After the initial transient, the lines tend to a constant gradient of approximately 
$-1/5$.
Finally, we examine the flow in the drop in the inertial–viscous regime using the self-similar solution for the Stokes streamfunction 
$\varPsi$ and gap-width profile 
$d$. Figure 10 shows the streamlines of 
$\varPsi$ in similarity variables, as well as the vertical velocity 
$u_z(r,0)$ and tangential velocity 
$u_r(r,0)$ on the surface, and the van der Waals disjoining pressure 
$\varPi = -1 / d(r)^3$. The van der Waals attraction drives a predominantly vertical flow with a small radial correction towards the origin. Vorticity 
$\omega \equiv \partial u_r/\partial z - \partial u_z / \partial r$, shown both on 
$z=0$ and in the drop, is generated along 
$z=0$ and is maximum at an order 
$1$ distance away from 
$r=0$ in the similarity frame, corresponding approximately to the inflection point in the gap-width profile. The vorticity diffuses into the drop and is also advected outwards in the similarity frame as the shrinking flow length scales leave the vorticity behind in the far field.
Figure 10. (a) Streamlines (solid lines) of the Stokes streamfunction 
$\varPsi$ and contours of the vorticity 
$\omega$ for the inertial–viscous self-similar flow inside the drop. (b) The inertial–viscous similarity solution for the normal velocity 
$u_z(r,0)$, tangential velocity 
$u_r(r,0)$, vorticity 
$\omega (r,0)$, and van der Waals disjoining pressure 
$\varPi = -1/d(r)^3$ on the surface of the drop.
4.3. Inertial regime
To investigate the possibility of an inertial regime, we attempted time-dependent simulations as in § 4.2 for purely inviscid drops, setting the viscous terms 
$E^2 \varPsi$ in the governing equation (4.4) and 
$2 \partial u_z / \partial z$ in the dynamic boundary condition (4.5) to zero. We found the inviscid dynamics to be unstable, with the radial length scale 
$r$ of the near-contact region quickly shrinking and becoming arbitrarily small relative to the gap width 
$d$. In these simulations, the observed radial length scale was regularized numerically only by the minimum grid spacing. In real drops, one would expect there to be some physical process that provides a short wavelength cut-off and regularizes the instability. In particular, for any small but non-zero viscosity, there is a radial length scale on which viscosity cannot be neglected, thus providing a cut-off scale at which the dynamics enters the inertial–viscous regime. We therefore conclude that there is no distinct inertial regime of jump-to-contact. Indeed, as shown in figure 8, in time-dependent simulation from 
$d(r,0)=(1+0.004r^2)^{1/2}$ with very large, but finite, initial increased Reynolds number 
$\mathcal {R} \sim 4000$, viscosity quickly became significant, with the increased Reynolds number suggesting that the inertial–viscous regime was reached before the minimum separation 
$d_0$ had decreased by one order of magnitude. Subsequently, the dynamics converged to the inertial–viscous similarity solution.
5. Discussion
We have shown that for jump-to-contact due to van der Waals attraction between fluid drops separated by an effectively inviscid exterior fluid, there are two distinct dynamical regimes distinguished by the forces – viscosity, or viscosity and inertia – that act at leading order alongside the strong van der Waals attraction. In each case, surface tension is negligible as contact is approached. The regimes are classified by an ‘increased’ Reynolds number 
$\mathcal {R} \sim \mathcal {H} a r^3/ Oh^2\, d^4$, which depends on the fluid properties through the ratio of the molecular and drop length scales 
$\mathcal {H}^{1/2}$ and the Ohnesorge number 
$Oh$, and on the current dimensional gap width 
$d$ and radial length scale 
$r$. As the length scales of the near-contact region evolve during jump-to-contact, 
$\mathcal {R}$ is not constant, thus the dynamical regime can change over time. Figure 11 shows the different dynamical regimes and how, depending on the Ohnesorge number, the relevant regime changes as the gap width decreases.
Figure 11. Diagram showing how the dynamical regime of jump-to-contact changes as the minimum dimensional gap width 
$d_0(t)$ decreases for different Ohnesorge numbers 
$Oh^2$ for fixed 
$\mathcal {H} = 5 \times 10^{-14}$ and 
$d_0(0) = a \mathcal {H}^{1/3}$. For 
$Oh^2 \gg \mathcal {H}^{1/6}\approx 0.006$, the dynamics is initially in the Stokes regime, with, for 
$Oh^2\ll 1$, a transition to the inertial–viscous regime when 
$d_0(t) \sim a \mathcal {H}^{1/2}/Oh^2$. For 
$Oh^2 \ll \mathcal {H}^{1/6}$, the dynamics is in the inertial–viscous regime throughout. At the dotted line, the minimum gap width reaches molecular length scales and the drops effectively reach contact.
Viscous drops with 
$Oh^2 \gg \mathcal {H}^{1/6}$ are initially in the Stokes regime, with 
$\mathcal {R} \ll 1$ and a leading-order force balance between viscosity and van der Waals attraction. BL22 showed that in this regime, the dynamics is given by a non-universal self-similar solution as contact is approached, with the gap-width and radial length scales both decreasing with the same 
$t{'^{1/3}}$ time dependence for time until contact 
$t'$. As both the gap-width and radial scales have the same time dependence, there is no interfacial steeping as contact is approached. Thus the gradient of the gap-width profile, which is initially shallow for 
$\mathcal {H} \ll 1$, remains shallow throughout, in spite of the diverging van der Waals attraction. However, the divergence of the velocity outpaces the decrease of the flow length scale such that the inertia in the drop becomes increasingly significant. For moderately viscous drops with 
$\mathcal {H}^{1/6} \ll Oh^2 \ll 1$, inertia becomes significant before contact, and a new, inertial–viscous, dynamical regime then takes over. Drops with 
$Oh^2 \gtrsim 1$ remain in the Stokes regime until contact.
Conversely, for less viscous drops with 
$Oh ^2\lesssim \mathcal {H}^{1/6}$, the dynamics is in an inertial–viscous regime throughout jump-to-contact, with 
$\mathcal {R} \sim 1$ and a three-way leading-order force balance between inertia, viscosity and van der Waals attraction. For initial increased Reynolds number 
$R \gg 1$, where viscous forces in the drop are initially much weaker than inertia, the radial length scale of the near-contact region decreases rapidly until 
$\mathcal {R}\sim 1$ and viscosity becomes significant. In the inertial–viscous regime, we have shown that there is a unique self-similar solution for the jump-to-contact dynamics, with the dimensional gap width 
$d$ and radial length scale 
$r$ decreasing as
\begin{equation} d \sim a \, Oh^{1/4}\,\mathcal{H}^{1/4} \left(\frac{\gamma t'}{a \mu}\right)^{3/8}, \quad r \sim a \, Oh\left(\frac{\gamma t'}{a \mu}\right)^{1/2}, \end{equation}
with time until contact 
$t'$. The similarity solution is seen to be stable in time-dependent simulations, and gives the generic behaviour as contact is approached. The three-way force balance does not admit any scaling symmetry, other than a rescaling of 
$t'$, which provides some justification for the uniqueness of the self-similar behaviour, and contrasts with the Stokes regime, where a time-independent radial rescaling gives rise to non-universal behaviour.
In the inertial–viscous regime, there is some interfacial steeping as contact is approached, with 
$d/r \propto t{'^{-1/8}}$. However, the growth rate is sufficiently weak that the gradient of the gap-width profile remains shallow, 
$d/r\ll 1$, throughout jump-to-contact for all but the least viscous drops, failing only if 
$Oh^2 \!\ll\! \mathcal {H}^{1/2}$, where 
$\mathcal {H}^{1/2}\approx 2 \times 10^{-7}\!\ll\! 1$ for millimetre-scale drops.
The jump-to-contact dynamics sets the surface profile of the drops at the point of contact and gives an initial profile on scales smaller than (2.4a,b) for the subsequent coalescence of the drops. The radial length scale of the near-contact region vanishes at contact, leaving behind the far field of the similarity solution as the surface profile. For less viscous drops with 
$Oh^2\ll 1$, the dynamics preceding contact is in the inertial–viscous regime, thus in dimensional variables, the gap-width profile at contact is 
$d \approx 0.57 ( \mathcal {H}/Oh)^{1/4} a^{1/4} r^{3/4}$. Conversely, for more viscous drops with 
$Oh ^2\gtrsim 1$, BL22 found that the Stokes regime gives a conical profile at contact, with the 
$O(\mathcal {H}^{1/6})$ angle set by the initial conditions of jump-to-contact. The surface profile in each regime is shown in figure 12.
Figure 12. The surface profiles after (a) inertial–viscous and (b) viscous jump-to-contact on the scale of the initial gap-width and radial length scales. The profile in each case is given by the outer behaviour of the self-similar solution, with (a) 
$h \propto r^{3/4}$, and (b) 
$h \propto r$. Note that in (b), the angle of the cone is arbitrary due to the non-universality of the self-similar behaviour in the Stokes regime.
There are no images of the drop profile on the length scale of this near-contact region during the very early stages of coalescence to confirm either of these initial profiles. However, both of these profiles would lead to slower rates of coalescence than would be expected for touching spheres: the driving surface tension in both cases is less effective as the curvature of the meniscus at the edge of the contact patch is smaller when matching to a far-field wedge from the Stokes regime than to a cusp, and smaller again when matching to the 
$h \propto r^{3/4}$ profile from the inertial–viscous regime. The early-time rate of coalescence on scales smaller than (2.4a,b) would still be proportional to 
$- {\log t}$ with these initial profiles, but with a smaller coefficient, which may partly help to explain the slower-than-expected early-time rates of coalescence observed in experiments (e.g. Paulsen et al. Reference Paulsen, Burton, Nagel, Appathurai, Harris and Basaran2012, Reference Paulsen2013).
We began the modelling in § 2 by considering two identical spheres with some small deformation. However, since the radial length scale decreases rapidly as contact is approached, the original gap-width profile is asymptotically negligible. Thus the analysis in this paper applies equally to jump-to-contact between fluids of any smooth shape across an inviscid gap, and jump-to-contact between a viscous fluid and a solid object. For example, the dynamics of jump-to-contact in atomic force microscopy will also follow the physical regimes and self-similar behaviour given in this paper.
Funding
E.B. gratefully acknowledges a UK Engineering and Physical Sciences Research Council studentship (grant no. 2267833).
Declaration of interests
The authors report no conflict of interest.
Data statement
All data accompanying this paper are available directly within the paper.
Appendix A. Calculation of 
$\varPhi$
To calculate the evolution of the streamfunction 
$\varPsi$ from (4.4), we need to know the vector potential 
$\varPhi \boldsymbol {e}_{\theta }/r$ for 
$\boldsymbol {\nabla }p$ in 
$z \ge 0$. In this appendix, we show how to determine 
$\varPhi$ from the pressure 
$p$ on 
$z=0$, which is known from (4.5).
The pressure 
$p$ satisfies Laplace's equation 
$\nabla ^2 p = 0$ in 
$z>0$, with 
$p$ regular on 
$z=0$ and 
$p =o(|\boldsymbol {x}|)$ as 
$|\boldsymbol {x}| \to \infty$. Thus 
$p$ in 
$z>0$ can be expressed as a boundary integral over 
$z'=0$,
\begin{equation} p(\boldsymbol{x})={-}\int_{z'=0}p(\boldsymbol{x}')\,\frac{\partial H(\boldsymbol{x}',\boldsymbol{x}) }{\partial z'}\,\mathrm{d}S', \end{equation}
where 
$H$ is the Green's function for Laplace's equation that satisfies 
$H=0$ on 
$z=0$, given by
\begin{equation} H(\boldsymbol{x}',\boldsymbol{x}) ={-} \frac{1}{4 {\rm \pi}} \left( \frac{1}{|\boldsymbol{x}'-\boldsymbol{x}|}- \frac{1}{|\boldsymbol{x}'-\boldsymbol{x}+2z\boldsymbol{e}_z|} \right) . \end{equation}
Here, the primes denote dummy variables. As the pressure is axisymmetric about the 
$z$ axis, the azimuthal integral in (A1) can be calculated analytically to give
\begin{equation} p(r,z) = \frac{2z}{\rm \pi} \int_0^\infty \frac{p(r',0)\, r'\,{E}(k)}{[(r-r')^2 + z^2][(r+r')^2 + z^2]^{1/2}} \, \mathrm{d}r' \quad \mbox{for}\ z>0, \end{equation}
where 
$k^2 = 4rr'/(z^2+(r+r')^2)$, and 
${E}(k)$ is the complete elliptic integral of the second kind.
The integral (A3) has an integrable singularity at 
$r'= r$ as 
$z \to 0$, which we can remove analytically as follows. We note that a constant pressure satisfies both Laplace's equation and the far-field and regularity conditions required in the derivation of (A3). We can thus set 
$p = 1$ in (A3) to obtain the identity
\begin{equation} \frac{2z}{\rm \pi} \int_0^\infty \frac{ r'\,{E}(k)}{[(r-r')^2 + z^2][(r+r')^2 + z^2]^{1/2}} \, \mathrm{d}r' = 1. \end{equation}The integrable singularity in (A3) is then removed by writing
\begin{equation} p(r,z)=p(r,0) + \frac{2z}{\rm \pi} \int_0^\infty \frac{[p(r',0)-p(r,0)] r'\,{E}(k)}{[(r-r')^2 + z^2][(r+r')^2 + z^2]^{1/2}} \, \mathrm{d}r' . \end{equation}
The pressure in 
$z>0$ is thus specified by its value on the boundary 
$z=0$, as given by the normal component of the dynamic boundary condition (4.5). We now wish to determine the scalar field 
$\varPhi$, which is defined, like a streamfunction for the solenoidal 
$\boldsymbol {\nabla }p$, by 
$\boldsymbol {\nabla }\times ( \varPhi \boldsymbol {e}_{\theta } /r)= \boldsymbol {\nabla } p$.
As 
$\boldsymbol {\nabla } p$ is also irrotational, 
$\boldsymbol {\nabla } \times \boldsymbol {\nabla } \times ({\varPhi } \boldsymbol {e}_{\theta }/r )= \boldsymbol {0}$, which for axisymmetric 
$\varPhi$ yields the azimuthal equation
\begin{equation} \left(\nabla^2-\frac{1}{r^2}\right)\frac{\varPhi}{r} \equiv \frac{1}{r}\,\nabla^2 \varPhi - \frac{2}{r^2}\,\frac{\partial \varPhi}{\partial r} = 0. \end{equation}
The definition of 
$\varPhi$ gives 
$\partial \varPhi / \partial r = r\,\partial p/ \partial z$, so (A6) implies that 
$\nabla ^2 \varPhi =2\,{\partial p}/{\partial z}$. We note also that 
$\nabla ^2 (zp) = 2\,\partial p / \partial z$ as 
$p$ and 
$z$ are harmonic. Combining these results, we obtain finally
\begin{equation} \nabla^2(\varPhi - z p) = 0 . \end{equation}
Having shown that 
$\varPhi - zp$ is harmonic as well as 
$p$, we can re-use (A5) to write down a boundary integral expression for 
$\varPhi$:
\begin{equation} \varPhi(r,z) - z\,p(r,z)=\varPhi(r,0) + \frac{2z}{\rm \pi} \int_0^\infty \frac{[\varPhi(r',0)-\varPhi(r,0)] r'\,{E}(k)}{[(r-r')^2 + z^2][(r+r')^2 + z^2]^{1/2}} \, \mathrm{d}r' . \end{equation} Thus 
$\varPhi (r,z)$ can be determined from 
$p(r,z)$ and the values of 
$\varPhi (r,0)$ on the boundary. As noted previously, 
$\partial \varPhi / \partial r = r\,\partial p/ \partial z$, which we integrate along the boundary to obtain
\begin{equation} \varPhi(r,0) =\int_0^r r'\left. \frac{\partial p}{\partial z}\right|_{(r',0)} \mathrm{d}r' + \varPhi(0,0). \end{equation}
Setting 
$\varPhi (0,0)=0$ is equivalent to setting the arbitrary constant 
$c_1=0$ in (4.4), which, as mentioned previously, makes no difference to the flow.
Finally, we note that the value of 
$\partial p/ \partial z$ on 
$z=0$ can be obtained from (A5). The integral in (A5) converges absolutely, so the normal derivative can be taken inside the integral to give
\begin{equation} \frac{\partial p}{\partial z} = \frac{2}{\rm \pi} \lim_{z \to 0} \int_0^\infty \frac{\partial}{\partial z}\frac{ z [p(r',0)-p(r,0)] r'\,{E}(k)}{[(r-r')^2 + z^2][(r+r')^2 + z^2]^{1/2}} \, \mathrm{d}r' \quad \mbox{on}\ z = 0. \end{equation}Given the pressure on the boundary from the normal component of the dynamic boundary condition (4.5), 
$\partial p/ \partial z$ on 
$z=0$ is determined by (A10). This, in turn, determines 
$\varPhi$ on 
$z=0$ through (A9). Furthermore, the boundary pressure determines 
$p$ in 
$z>0$ through (A5). Finally, 
$\varPhi$ is determined throughout 
$z > 0$ from its boundary value and the pressure field by (A8). Thus 
$\varPhi$ is determined completely by the normal component of stress on the boundary.
Appendix B. The numerical scheme
In this appendix, we provide further details of the numerical methods and parameters used in the simulations of § 4.
At each time step, the current values of the Stokes streamfunction 
$\varPsi$ are known on a grid, with the gap-width profile also known at the grid points along 
$z=0$. We calculate the viscous diffusion 
$E^2 \varPsi$ at each grid point using five-point central finite differencing. Symmetry on 
$r=0$ and the tangential dynamic boundary condition (4.6) on 
$z=0$ are used to calculate the normal derivatives on the respective boundaries. The derivatives required for the normal surface velocity 
$u_z=-r^{-1}\,\partial \varPsi / \partial r$ on 
$z=0$ and the normal viscous surface stress 
$2\,\partial u_z / \partial z= -2r^{-1}\,\partial ^2 \varPsi / \partial r\,\partial z$ on 
$z=0$ are calculated likewise.
The forcing term 
$\varPhi$ in (4.4) is calculated at each grid point using the integrals (A8)–(A10), with the pressure on the boundary known from (4.5), as follows. We use Simpson's method to integrate between consecutive grid points, evaluating the integrand at an adaptive number of subdivision points where the values of 
$p(r,0)$ and 
$\varPhi (r,0)$ are interpolated using a quintic spline. The number of subdivision points is doubled successively until the relative change to the result when doubling is less than 
$10^{-7}$. We restrict the integration range to the finite grid, and note that the contribution to each integral from 
$r>r_{{max}}$, for maximum grid extent 
$r_{{max}}$, is less than 
$r^*/r_{{max}}$ times the leading-order contribution from near 
$r=0$, where 
$r^*$ is the radial length scale of the near-contact region.
We time march the Stokes streamfunction 
$\varPsi$ and gap width 
$d$ using the forward Euler method. We take the time step to be the minimum of 
$d_0/200 \dot {d\hskip 1.4pt}\hskip -1.4pt_0$ and 
$\varDelta ^2/8$, where 
$\varDelta$ is the minimum separation between grid points. The first condition ensures that the surface does not evolve by more than 
$0.005\,\%$ of its current value in any time step, and the second ensures that the forward Euler scheme is numerically stable.
We continue to time step the solution on the same grid until the minimum gap width 
$d_0$ reaches 
$90\,\%$ of its starting value on the current grid. We then update the grid by rescaling both the vertical and radial directions with a characteristic radius 
$r^*$, defined as in § 4.2 as the radius at which the vertical velocity on 
$z=0$, which is peaked around 
$r=0$, is half of its maximum value: 
$\dot {d\hskip 1.4pt}\hskip -1.4pt(r^*,t) = \dot {d\hskip 1.4pt}\hskip -1.4pt_0(t)/2$. We use a quintic spline to interpolate the Stokes streamfunction 
$\varPsi$ and the gap width 
$d$ on the new grid.
For the simulations presented in this paper, we used a grid obtained from a fixed grid 
$\{ (\eta _i, \eta _j)\}_{ 0 \le i, j \le n}$ with uniform spacing 
$\delta \eta$ by the map
\begin{equation} r_i= R(t) \sinh(\eta_i),\quad z_j = R(t)\, \eta_j, \end{equation}
where 
$R(0)$ is an initial parameter, and 
$R(t)$ is determined by the rescaling procedure given above. We took 
$\delta \eta = 0.075$, 
$n = 135$ and 
$R(0) = 0.5$. To verify these choices of grid parameters, we individually varied 
$\delta \eta$, 
$n$ and 
$R(0)$, and considered the effect on the evolution of 
$\dot {d\hskip 1.4pt}\hskip -1.4pt_0$ as 
$d_0$ evolved over one order of magnitude. We found that the relative cumulative error between simulations with 
$\delta \eta = 0.075$, 
$n=135$ and 
$R(0) = 0.5$, and simulations with 
$\delta \eta = 0.075$, 
$n=125$ and 
$R(0) = 1$ (
$50\,\%$ of the resolution, same maximum extent), 
$\delta \eta = 0.1$, 
$n=100$ and 
$R(0) = 0.5$ (
$75\,\%$ of the resolution, same maximum extent), or 
$\delta \eta = 0.075$, 
$n=151$ and 
$R(0) = 0.5$ (same resolution, 
$10\,\%$ of the maximum extent), was less than 
$10^{-4}$. We similarly varied the choice of grid in the 
$r$ and 
$z$ directions independently. We note that with this choice of grid, 
$r^*/r_{{max}} = 1/R(0)\sinh (n \, \delta \eta ) \approx 10^{-4}$, so neglecting the far-field contribution to the boundary integrals used to calculate 
$\varPhi$ is justified.
Open access
