Capillary Breakup of a Liquid Bridge: Identifying Regimes and Transitions

Computations of the breakup of a liquid bridge are used to establish the limits of applicability of similarity solutions derived for different breakup regimes. These regimes are based on particular viscous-inertial balances, that is different limits of the Ohnesorge number $Oh$. To accurately establish the transitions between regimes, the minimum bridge radius is resolved through four orders of magnitude using a purpose-built multiscale finite element method. This allows us to construct a quantitative phase diagram for the breakup phenomenon which includes the appearance of a recently discovered low-$Oh$ viscous regime. The method used to quantify the accuracy of the similarity solutions allows us to identify a number of previously unobserved features of the breakup, most notably an oscillatory convergence towards the viscous-inertial similarity solution. Finally, we discuss how the new findings open up a number of challenges for both theoretical and experimental analysis.


Introduction
The breakup of liquid volumes is a process that is ubiquitous throughout industry and nature. Recently, this process has attracted significant attention due to its importance for the functioning of a range of microfluidic technologies where one would like to be able to control the generation of uniform-sized droplets which then become building blocks, as in three-dimensional (3-D) printers (Derby 2010), or modes of transport for reagents, as in lab-on-a-chip devices (Stone, Stroock & Adjari 2004). In order to optimise the functioning of these technologies it is key to be able to understand the physical mechanisms controlling the breakup process and, ideally, to have a tool capable of predicting the response of the system to alterations in operating conditions. From a theoretical perspective, the breakup phenomenon has attracted interest due to both its technological importance and its status as a free-surface flow which 30 Y. Li and J. E. Sprittles exhibits finite-time singularities (Eggers 1997). In the former case, motivated by a desire to describe engineering-scale systems, the focus has been on capturing the global dynamics of the process using multiphysics computational fluid dynamics (CFD) codes, with the actual breakup usually under-resolved due to the inherently multiscale nature of the problem. In contrast, in the latter case most interest has revolved around using asymptotic methods to study the micromechanics of the breakup process, with the larger-scale flow often not considered. Potentially, the advances made in understanding the micromechanics could be fed into large-scale CFD codes in order to exploit the advantages of both approaches, but this is yet to be achieved.
The failure of CFD packages to reliably capture breakup phenomena has been highlighted in Fawehinmi et al. (2005) for drop formation, where the appearance of satellite drops depended on grid resolution, and Hysing et al. (2009) for a two-dimensional (2-D) rising bubble in a regime where breakup was anticipated. In Hysing et al. (2009), six different codes were tested for the same problem and all gave different results ranging from no observed breakup through to multiple bubble detachments. In other words, the codes were unreliable for this process. The conclusion of the authors provides motivation for the current work: 'Although the obtained benchmark quantities were in the same ranges, they did not agree on the point of break up or even what the bubble should look like afterwards, rendering these results rather inconclusive. To establish reference benchmark solutions including break up and coalescence will clearly require much more intensive efforts by the research community'.
Computational approaches cannot, in principle, resolve the thinning of a liquid volume down to the point at which breakup occurs, when the thread is infinitesimally small. Instead, at some point the pre-breakup process must be terminated, by 'cutting' the liquid thread joining two volumes, and a post-breakup state must be initiated based on the final solution in the pre-breakup phase. Whilst some computational codes do this process 'automatically', notably those using Eulerian meshes such as volume-of-fluid, there is no guarantee that this approach accurately represents the physical reality. Theoretical approaches to this problem have been considered in different regimes and are described in Eggers (2014). To implement an appropriate post-breakup solution, one must identify which pre-breakup regime the process is in, based on the scale for the cutoff. Therefore, knowing what regime a given breakup process is in and determining the transitions between these regimes are critical for the topological change to be correctly handled. The complexity of this procedure has been revealed in Castrejón-Pita et al. (2015); see § 2.2.3 for details, where multiple regimes were discovered for each breakup event.
Here, by computing the pre-breakup solution to the smallest possible scales that our simulation will allow, we will measure the accuracy of various similarity solutions proposed in the literature and, importantly, establish their limits of applicability. To do so, rather than just identifying power-law behaviour, methods will be developed to determine the regions of phase space where similarity solutions proposed for the different 'regimes' accurately approximate the full solution. The procedure is intended to leave no ambiguity as to when a certain regime is encountered. This attempt at a more rigorous approach will lead us to discover previously unobserved features of the breakup process that open up several new avenues of enquiry.

Regimes and transitions in the breakup phenomenon
This article focuses on the breakup of incompressible Newtonian liquids with constant viscosity µ and density ρ, which have a constant surface tension σ at their liquid-gas interface and are surrounded by a dynamically passive gas. The breakup is considered axisymmetric with the axis of the thread along the z-axis of a cylindrical coordinate system (r, z), as shown in figure 1 for a liquid bridge geometry, with characteristic length scale R. In general, the flow is governed by the Navier-Stokes equations, but when inertial forces are weak the dominant balance of viscous and capillary forces gives a characteristic speed of breakup of U = σ /µ so that the capillary number Ca = µU/σ = 1. Then the appropriate dimensionless number characterising breakup is the Ohnesorge number Oh, which is obtained by substituting U into the Reynolds number Re = ρUR/µ = ρσ R/µ 2 = Oh −2 . The Ohnesorge number is the square root of the ratio of a viscous length scale µ = µ 2 /(ρσ ) to the characteristic scale of the system R. Then Oh −1 = 0 is Stokes flow and Oh = 0 gives inviscid flow. The scales used so far are based only on global quantities. To analyse the breakup process locally it is useful to introduce a time-dependent local Reynolds number Re local (t) = ρU z L z /µ based on instantaneous axial scales U z (t) and L z (t) (which will be larger than, or comparable to, radial scales). The practicalities of calculating this quantity are discussed in § 5.1. This local parameter will indicate the dominant forces during a breakup event and hence should identify the flow regime.

Regimes of breakup
Near breakup much research has focused on finding similarity solutions for the different regimes, all of which are categorised in § 4 of Eggers & Villermaux (2008). A summary of the key results is provided here, with particular attention paid to the dependence of the minimum thread radius r min , the maximum axial velocity w max and the local Reynolds number Re local on the time from breakup τ = t b − t, where t b is the breakup time. Quantities have been made dimensionless with characteristic scales for lengths, velocities and time of R, σ /µ and µR/σ .

The inertial regime
In the 'inertial regime', henceforth the I-regime, where there is a balance between inertial and capillary forces, dimensional analysis (Keller & Miksis 1983;Brenner et al. 1997) of inviscid flow (Oh = 0) gives as τ → 0 that

Y. Li and J. E. Sprittles
where A I is a constant of proportionality, given in Eggers & Villermaux (2008) as A I ≈ 0.7. Notably, computations performed in the inertial regime have never listed A I and so it will be of interest to determine this value so that (2.1) is precisely determined and becomes a predictive tool. It will be useful to note, in terms of r min , that Re local ∼ Oh −1 r 1/2 min and w max ∼ Oh r −1/2 min . A feature of the inertial regime is the 'overturning' of the free surface, observed in both experiments (Chen, Notz & Basaran 2002) and simulations in this regime (Schulkes 1994), which invalidates attempts at a slender description of the thread. In Day, Hinch & Lister (1998) it was predicted that at breakup the free surface forms a double-cone shape with angles of 18.1 • and 112.8 • with the z-axis. This phenomenon has been observed experimentally in Castrejón-Pita et al. (2012) and computationally in Wilkes, Phillips & Basaran (1999).

The viscous regime
The solutions derived for both the viscous and viscous-inertial regimes, henceforth referred to as the V-regime and the VI-regime, were derived in a one-dimensional approximation of the Navier-Stokes equations which relies on the thread remaining 'slender' as the breakup is approached. For a free surface represented by r = r(z, t) this means that the gradient of the free surface must remain small, ∂r/∂z 1.
For the V-regime, the similarity solution for Stokes flow (Oh −1 = 0) was derived in Papageorgiou (1995a,b) and is given by where β ≈ 0.175. Recently, Eggers (2012) showed that this is the only similarity solution of infinitely many possible ones which remains stable to small perturbations. The free-surface profile associated with this similarity solution remains symmetric about the pinch point and the local profile depends on an external length scale (i.e. it is not 'universal'), so that it is a self-similar problem of the second kind (Barenblatt 1996). The persistence of the V-regime for high-viscosity liquids has been observed experimentally in McKinley & Tripathi (2000).

The viscous-inertial regime
As τ → 0, from (2.1) and (2.2) we find in the I-regime that Re local → 0 and in the V-regime that Re local → ∞. This contradicts the initial assumptions behind their derivations (Lister & Stone 1998). Therefore, neither the I-nor the V-regime is valid right up to breakup, and a regime where viscous and inertial effects are in balance, so that Re local ∼ 1, must be considered. In this VI-regime, the Navier-Stokes equations are required. It was shown in Eggers (1993) that for the VI-regime a 'universal' set of exponents exist for which with a highly asymmetric free-surface shape joining a thin thread to a much steeper 'drop-like' profile. It was later demonstrated in Brenner, Lister & Stone (1996) that this solution is the most favourable of a countably infinite set of similarity solutions. Notably w max = max(w) is the maximum velocity out of the thin thread (w > 0 in what follows), rather than max|w|, which scales in the same way but has a pre-factor 3.1 (Eggers 1993). In terms of r min this gives w max = 0.3Oh r −1/2 min .
Capillary breakup 33 2.2. Transitions between regimes It has been established that whichever regime the breakup process starts in, it will eventually end up in the VI-regime. However, it is possible that the transition from either the V-or the I-regime will occur at such small scales that this regime is irrelevant, at least from a practical perspective. For example, experiments in Burton, Rutledge & Taborek (2004) for millimetre-sized mercury drops show that the similarity solution in the I-regime adequately describes the breakup down to the nanoscale. Therefore, determining the transitions between the regimes is an important part of understanding the breakup process.
Phase diagrams which identify the regions of (Oh, r min ) space associated with different regimes have been constructed for the related problem of drop coalescence in Paulsen et al. (2012), Paulsen (2013) and Sprittles & Shikhmurzaev (2014b) though, intriguingly, these publications disagree on the number of regimes -whereas for the breakup phenomenon an equivalent diagram is currently lacking. Most progress in this direction has been made in Castrejón-Pita et al. (2015), where sketches of two characteristic phase space trajectories were presented (in their figure 1E). Here, we provide a systematic exploration of (Oh, r min ) space in order to construct a phase diagram for breakup which is equivalent to those obtained for coalescence. To do so, the scalings proposed for the transitions between the various regimes as r min , τ → 0 are considered.

Viscous to viscous-inertial transition (V→VI)
When Oh 1, the V-regime adequately describes the initial stages of breakup. However, the inertial term in the Navier-Stokes equations does not remain negligible as r min → 0 and so a transition to the VI-regime occurs at a bridge radius r V→VI min . In Basaran (2002) and Eggers (2005), by balancing the solutions in the V-and VI-regimes, i.e. by taking Re local ∼ 1 in (2.2), it was shown that this should occur when However, experimental evidence in Rothert, Richter & Rehberg (2003) appears to contradict this result, finding instead that r V→VI min is constant, i.e. independent of Oh. Investigating these regimes computationally should provide new insight into the transition.

Inertial to viscous-inertial transition (I→VI)
When Oh 1 the Euler equations accurately capture the initial stages of breakup until the local Reynolds number in (2.1) drops to Re local ∼ 1, which occurs when τ ∼ Oh 2 so that r I→VI min ∼ Oh 2 . (2.5) This cross-over was first confirmed computationally in Notz, Chen & Basaran (2001) and experimentally in  for the case of a water-glycerol mixture with Oh = 0.16, but has never been studied systematically across a range of Oh. The paths of typical breakup events are given by arrowed lines 1 and 2, which show how different regimes are encountered as r min decreases. At higher Oh (path 2) a V→VI transition is both expected and computed, whilst at lower Oh (path 1), a single I→VI transition was expected but a more complex behaviour was computed, with a low-Oh V-regime encountered before the VI-regime, so that I→V→VI transitions are found, as discovered in Castrejón-Pita et al. (2015). The dashed lines show the scaling of the transitions between the different regimes as r min → 0, with I→V and V→VI transitions scaling as r min ∼ Oh 2 and V→VI occurring when r min ∼ Oh −3.1 .
given Oh) from the top axis downwards (through decreasing r min ), e.g. see paths 1 and 2. For small Oh (path 1) one has an I-regime crossing into a VI-regime when r min ∼ Oh 2 , whilst for large Oh (path 2) one has a V-regime crossing into the VI-regime at r min ∼ Oh −3.1 . However, a recent publication by Castrejón-Pita et al. (2015) shows that the picture can be more complex: the breakup can pass transiently through multiple different regimes. For example, Castrejón-Pita et al. (2015) show that, counter-intuitively, at Oh = 0.23 there is no I→VI transition but rather an I→V→VI transition, so that a new unexpected low-Oh V-regime is encountered. Similarly, a high-Oh I-regime is observed. This behaviour is sketched out qualitatively in a phase diagram (their figure  1E) showing how Re local varies as breakup is approached for different values of Oh.
In this work, we will use computations to build the first fully computed and quantitatively constructed phase diagram for the breakup so that the nature of the transitions can be established. To orientate the reader over the forthcoming sections, figure 2(b) gives a preview of the computed phase diagram, in a form that is only asymptotically accurate as r min → 0 but serves to neatly illustrate the computed flow transitions. The result shows clearly how the new phase diagram differs from the expected behaviour (figure 2a) due to the appearance of the low-Oh V-regime sandwiched between the I-and VI-regimes, so that path 1 now first encounters an I→V transition. In contrast to Castrejón-Pita et al. (2015), no evidence for the high-Oh I-regime can be seen. These features will be considered in further detail over the forthcoming sections.
In § § 7-9 we will show how the phase diagram in figure 2(b) was actually constructed, but before doing so we must specify the problem considered.

Problem formulation
A liquid bridge geometry (figure 1) is used with the liquid trapped between two stationary solid discs of (dimensional) radius R a distance 2R apart and surrounded Capillary breakup 35 by a dynamically passive gas. This set-up allows us to isolate the breakup dynamics, limit the elongation of the domain (e.g. in contrast to dripping phenomena) and retain an experimentally realisable set-up. A future work will consider the effects of geometry through elongation by allowing the plates to move apart. The contact line where the liquid-gas free surface meets the solid remains pinned at the disc's edge throughout. Assuming gravitational effects are negligible allows us to consider a plane of symmetry at z = 0 of a cylindrical polar coordinate system (r, z, θ ) so that, using R as a characteristic scale for lengths, the solid is located at dimensionless position z = 1 and the free surface is pinned at (r, z) = (1, 1). The extension to include gravity is not a difficult one, but it does not add significant value to a study focused on the small-scale breakup.
Using U = σ /µ as a scale for velocity, T = R/U as a scale for time and µU/R for pressure, the (dimensionless) Navier-Stokes equations are where the stress tensor is Here, u and p are, respectively, the velocity and pressure in the liquid, and the Ohnesorge number is Oh = µ/ √ ρσ R. On the free surface, whose location f (r, z, t) = 0 must be obtained as part of the solution, the kinematic equation is applied alongside the usual balance of fluid stresses with capillarity in the directions tangential and normal to the free surface: n · P · (I − nn) = 0, n · P · n = ∇ · n, where n is the inward normal, I is the metric tensor of the coordinate system and p is taken relative to the (constant) gas pressure. At the liquid-solid boundary, no-slip and impermeability conditions are applied, u = 0 and the free surface is pinned at the contact line f (1, 1, t) = 0.

Initial conditions and initiation of breakup
The dynamics of breakup is studied by generating a thread shape that is taken to the edge of its Rayleigh-Plateau stability limit (Slobozhanin & Perales 1993;Lowry & Steen 1995), so that small perturbations will trigger breakup. The movie 'FullBreakup' in the supplementary material available at http://dx.doi.org/10.1017/jfm.2016.276 shows a typical breakup event resulting from this approach. In this way, the breakup dynamics can be computed independently of a driving force such as the flux from a nozzle, as in jetting, or the pull of gravity, as in dripping (Rubio-Rubio, Sevilla & Gordillo 2013). This reduces the dimensionality of the system to just one controlling parameter Oh so that constructing a phase diagram becomes a tractable task.
The required initial thread shape can be generated by either solving the Young-Laplace equation for the free-surface position (Fordham 1948;Brown & Scriven 1980; 36 Y. Li and J. E. Sprittles Thoroddsen, Takehara & Etoh 2005) or using the finite element code, described below, in a quasi-static mode, to generate this shape by gradually removing fluid from the thread in a manner similar to the experimental set-up described in Meseguer & Sanz (1985). The latter method is used, as this also enables the shape on the stable-unstable boundary to be established. The result is the first profile in figure 3(a). The instability can be triggered either by continuing to gradually remove fluid from the thread or by very slowly moving the plates apart, and both approaches were shown to produce indistinguishable results, thus confirming that the only control parameter is Oh.

Multiscale finite element computations
The mathematical model formulated in § 3 is for an unsteady free-surface flow influenced by the forces of inertia, viscosity and capillarity. To study this system over the entire range of Oh requires computational methods.

Previous computational studies
The breakup of liquid volumes has been studied computationally using the Stokes (Gaudet, McKinley & Stone 1996;Pozrikidis 1999), Euler (Schulkes 1994;Day et al. 1998) and Navier-Stokes equations for both Newtonian (Ashgriz & Mashayek 1995;Popinet 2009) and non-Newtonian (Li & Fontelos 2003;Bhat et al. 2010) liquids. Considerable progress has been made by the group at Purdue University using finite element methods for axisymmetric flows such as dripping (Ambravaneswaran, Phillips & Basaran 2000), jetting (Ambravaneswaran et al. 2004), drop-on-demand ) and liquid bridge breakup (Suryo & Basaran 2006) with an array of different physical effects such as non-Newtonian fluids (Yildirim & Basaran 2001), surfactant dynamics (McGough & Basaran 2006) and electric field effects (Collins et al. 2008). These algorithms were the first to demonstrate a number of experimental features such as overturning of the free surface for a viscous fluid (Wilkes et al. 1999), transitions between scaling regimes  and the identification of multiple regime transitions (Castrejón-Pita et al. 2015).
Notably, it has been shown that with sufficient care, the results of sharp interface methods, such as those implemented at Purdue, can sometimes be recovered by the diffuse interface approach, in which topological changes are easier to handle computationally. In particular, the diffuse interface method developed in Yue et al. (2004) was shown in Zhou, Yue & Feng (2006) to accurately recover the sharp interface results of Wilkes et al. (1999). Furthermore, in Zhou et al. (2006), the flexibility of the diffuse interface method was demonstrated by computing compound drop formation, where additional complexity arises from the need to track two free surfaces. Despite these advances, we will implement a sharp interface approach (see § 4.2) whose reliability and accuracy has been repeatedly confirmed.
Despite the numerous successes of computational schemes in accurately predicting many of the global features of breakup, there have been fewer investigations of the local dynamics and its comparison with the similarity solutions. One of the most impressive works in this direction is Suryo & Basaran (2006), where scales comparable to those resolved in the current work (r min < 10 −4 ) were captured and the scaling behaviour of both Newtonian and power-law liquids was investigated. Other progress includes results in Notz et al. (2001) and , where I→VI transitions were shown, and Castrejón-Pita et al. (2015), in which multiple transitions were observed (see § 2.2.3). Our work will build on these previous investigations by systematically quantifying the accuracy of the similarity solutions across the Downloaded from https://www.cambridge.org/core. IP address: 52.11.211.149, on entire parameter space (of Newtonian liquids) in order to build the first quantitatively constructed phase diagram for breakup.
In principle, computations are the ideal tool with which to map a phase diagram for the breakup process and to investigate in which regions each of the similarity solutions is accurate; a similar procedure was used for the coalescence phenomenon in Sprittles & Shikhmurzaev (2014b). The difficulty of this procedure is that many decades of r min are required in order to reliably determine the scaling behaviour of different regime boundaries (if they exist). Although one may hope to compute r min until its behaviour falls into the VI-regime, in practice this may occur at scales below the possible computational resolution, e.g. r min < 10 −5 (Burton et al. 2004). Therefore, establishing the range of applicability of the similarity solutions using computations and then using these solutions to carry the dynamics to scales below realisable computational resolution seems to be a promising strategy for capturing breakup.

Computational approach
A focus of this work is to resolve the spatial and temporal dynamics of the pinch-off process to the smallest scales possible in order to compare with the similarity solutions proposed in the literature for the final stages of breakup. The approach used is based on the finite element framework originally developed in Sprittles & Shikhmurzaev (2012b to capture dynamic wetting problems and subsequently used to study the coalescence of liquid drops (Sprittles & Shikhmurzaev 2012a, 2014a, including two-phase calculations (Sprittles & Shikhmurzaev 2014b); drop impact phenomena (Sprittles & Shikhmurzaev 2012a); the detachment of bubbles from an orifice (Simmons, Sprittles & Shikhmurzaev 2015); and dynamic wetting in a Knudsen gas (Sprittles 2015). These flow configurations are all inherently multiscale, either due to the disparity of length scales in the problem formulation or because of the dynamics of the process itself, which generates small scales during its evolution, as is the case for breakup phenomena. The framework allows the global flow in these problems to be resolved alongside localised smaller scales, without relying on particular limits of Oh.
As a step-by-step user-friendly guide to the implementation of this computational framework has been provided in Sprittles & Shikhmurzaev (2012b), only the main details will be recapitulated here alongside some aspects which are specific to the current work. The code uses the arbitrary Lagrangian Eulerian scheme, based on the method of spines (Ruschak 1980;Kistler & Scriven 1983), to capture the evolution of the free surface in 2-D or 3-D axisymmetric flows. In order to keep the problem computationally tractable, the mesh is graded so that small elements can be used near the pinch-off region whilst larger elements are used where scales associated with the global flow are present.
The mesh starts with ≈4000 triangular elements (≈500 surface nodes) and then adaptively refines as the breakup progresses by adding spines whenever elements become too deformed. At the end of the computation, the number of elements is in the range 7000-50 000 (800-7000 surface nodes), depending on the particular breakup requirements. For example, for a computation at Oh = 10 −3 , the mesh is refined due to the formation of a 'corner' in the free surface (see figure 15), resulting in 869 surface nodes with a minimum surface element length of 3 × 10 −5 . In contrast, at Oh = 0.16, refinement is required due to the formation of a thin thread of liquid (figure 9) which requires 5263 surface nodes with surface elements in this region all having length ≈10 −4 . Overturning of the free surface is permitted by using angled spines (i.e. not only horizontal) whose slope varies along the thread.
Refinement of the mesh is restricted by controlling the smallest permitted element size h min , which affects the final number of elements. Reducing h min allows the code to converge to smaller r min so that more of the breakup is recovered. By reducing h min from 10 −2 down to 10 −5 , convergence of the scheme under spatial refinement has been established, with each decrease in h min revealing more of the solution (smaller r min ) but producing curves which are graphically indistinguishable from those obtained on cruder meshes, where both solutions exist.
The result of the spatial discretisation is a system of nonlinear differential algebraic equations of index two (Lötstedt & Petzold 1986), which are solved using the second-order backward differentiation formula, whose application to the Navier-Stokes equations is described in detail in Gresho & Sani (1999). The time step automatically adapts during a simulation to capture the appropriate temporal scale at each instant.
In the problem considered, the smallest scale that needs resolving is the minimum neck radius r min , and eventually this becomes so small that accurate converged solutions can no longer be obtained. For the entire range of Oh considered, computations remain accurate down to r min = 10 −4 and in certain cases they can go as far as r min = 10 −5 .

Quantitatively identifying regimes
Computations are unable to resolve down to r min = 0, i.e. the point at which the topological change occurs, so that the precise time t = t b of pinch-off is unknown. This can be problematic when comparing to scaling laws of the form r min = Aτ B + C, where τ = t b − t requires the time of breakup, as small changes in t b can have a large effect on a curve of r min versus τ , which is then used to determine which regime a breakup process is in. Similar conclusions were reached for coalescence in Thoroddsen et al. (2005). In the final regime one has C = 0, otherwise for 'transitional' or 'transient' regimes C must be determined as well. In general, there is no easy method to overcome these issues without assuming a particular functional form for the breakup.
However, in the V-and VI-regimes it is possible to exploit r min being a linear function of τ (B = 1) so that the speed at which the minimum thread radius evolves is constant, i.e. independent of both t b and C. Therefore, to identify these regimes the 'speed of breakup' given byṙ min = dr min /dt can be compared to the predictions of (2.2) thatṙ min = −0.071 and (2.3) thatṙ min = −0.030 as a function of r min without needing to know t b . This has the further advantage that C does not have to be fitted when a regime is observed transiently, as in Castrejón-Pita et al. (2015).
The same approach cannot be used to identify the inertial regime where B = 2/3, so that the speed of breakup is no longer constant. Therefore, we try a similar idea to that above and determine an Oh-independent quantity which should be linear in the inertial regime. This relies on C ≈ 0, so that this regime must be close to breakup. Then, given that r min = A I (Oh τ ) 2/3 , we introduce l min = Oh −1 r 3/2 min , which will satisfy l min = A 3/2 I τ if the expected scaling holds. Differentiating this quantity with respect to time givesl min = dl min /dt = −A 3/2 I , i.e. a time-independent constant. A key aspect of this work is to develop techniques to quantify when the breakup dynamics is in a certain regime (Sprittles & Shikhmurzaev 2014b). To do so, we must define what it means to be 'in a regime'. We choose to define a breakup process to be in the V-or VI-regime when the speed of breakupṙ min is within 0.015 of Capillary breakup 39 that predicted by (2.2) or (2.3), respectively. In other words, speeds in the rangė r min = −0.071 ± 0.015 indicate V-regime dynamics whilstṙ min = −0.030 ± 0.015 for the VI-regime. In a similar way, the inertial regime is defined by points at whichl min is within 0.15 of A 3/2 I , where we will later see that A 3/2 I = 0.5. Choosing smaller (larger) margins to define each regime will shrink (expand) the area of a given regime in the phase diagram but should not alter the qualitative picture. The particular values chosen avoided regimes overlapping and enabled us to satisfactorily map the phase diagram. The method for identifying regimes is illustrated in the Appendix.
Rather than assuming that the regimes fill all of phase space, we will be careful to split the entrance to and exit from a particular regime. For example, the exit from the V-regime will be labelled as r V→ min and the entrance to the VI-regime as r →VI min , rather than assuming a priori that these coincide at r V→VI min .
5.1. Analysis of the regimes To understand the dynamics in the different regimes and determine when the assumptions behind the similarity solutions of § 2.1 are valid, local to the breakup point we analyse the assumptions made about (i) the relative magnitudes of viscosity and inertia and (ii) the slenderness of the thread. In practice, this means calculating (i) the local Reynolds number Re local and (ii) a measure for the slenderness of the free-surface profile.
To measure these characteristics, one needs to define an appropriate axial length scale L z and velocity scale U z near the pinch point. The obvious choice for U z is the maximum axial velocity in the thread, U z = w max , which will occur near, but not at, the position where r = r min . Therefore, a characteristic axial scale can be defined as the vertical distance between the positions of minimum radius and maximum velocity, so that L z = |z(r = r min ) − z(w = w max )|. This defines (i) the local Reynolds number as Re local = Oh −2 U z L z and (ii) the slenderness as the ratio of the characteristic radial length scale L r = r min to L z , so that the geometry is slender when L r /L z 1. Importantly, w max = max(w) is used rather than max |w| (see figure 9c) to avoid discontinuities in L z , and hence Re local , which occur when the position of max |w| jumps from being above the pinch point to below it (as in the V→VI transition).
Other possible definitions of U z and L z exist (see Castrejón-Pita et al. 2015), so it is not a priori clear we have made the correct choices. However, computations will show that our definitions allow us to reproduce the expected scalings for Re local in both the I-regime (2.1) and the V-regime (2.2), which serves to validate our approach.

Overview
In § § 7-9, respectively, the dynamics of breakup for large (Oh > 1), intermediate (0.1 < Oh < 1) and small (Oh < 0.1) Ohnesorge numbers are analysed separately and used to establish the limits of applicability of the different similarity solutions described in § 2.1. This will allow us to define the regimes of breakup on a phase diagram.
In each section, the initial focus will be on a breakup event which highlights the features of the prominent regime, namely the V-regime (Oh = 10), VI-regime (Oh = 0.16) and I-regime (Oh = 10 −3 ), with a movie of each breakup in the supplementary material. In each case, snapshots of the free-surface shape, axial velocity and pressure distributions in the breakup region will be presented. Evidence that the breakup is in a given regime will be provided by comparing the computational results with the similarity solution's predictions for the evolution of the free-surface shape, minimum (1, blue), r min = 10 −2 (2, red), r min = 10 −3 (3, green) and r min = 10 −4 (4, brown).
bridge radius r min , maximum axial velocity w max and scaling of the local Reynolds number Re local . Having confirmed the existence of a regime, we then identify its boundaries on the phase diagram to determine the regime transitions. As explained in § 5, these boundaries are specified by calculating the values ofṙ min andl min as a function of r min . Having discovered the transitional behaviour, the local measures Re local and slenderness L r /L z (introduced in § 5.1) are studied in an attempt to rationalise these findings.
The analysis in § § 7-9 will not only allow us to construct a phase diagram (cf. figure 18) for the process but also open up new questions about the breakup process; to avoid distracting from the main focus in § § 7-9, these aspects will be considered in more detail in § 11.

Viscous-dominated flow (large Oh)
For Oh = 10 free-surface profiles in figure 3(a,b) show the development of a thin thread of liquid whose smallest minimum radius remains at z = 0 throughout, so that no satellite drops form. The axial velocity (figure 3c) at the free surface w fs , which is approximately r-independent, shows a rapid drainage from the thread driven by a pressure gradient from z = 0 (figure 3d).
Evidence that this breakup is in the V-regime can be seen in figure 4(b), where the free-surface profile obtained at r min = 10 −4 is almost indistinguishable from the  (1993), giving (r, z) = (τ φ VI , z 0 + Oh τ 1/2 ξ VI ) with no free parameters, and in (b) with (ξ V , φ V ) from Papageorgiou (1995b) and Eggers (1997), so that (r, z) = (τ φ V , Oh 2−2β τ β z ξ V ) with z = 1.85 × 10 −3 . In (a), the height of the pinch point z 0 = 0.268 for the VI-regime is calculated using the expression in Eggers (1993) for the drift of this point z(r = r min ) = z 0 − 1.6Oh τ −1/2 (shown in figure 11 as the dashed line). Assuming that the breakup remains in the V-regime, the data can be extrapolated to obtain the breakup time t b and, as a result, plot r min against τ in figure 5 where there is excellent agreement with the similarity solution of (2.2) for r min < 10 −2 .

Identification of the V-regime
When Oh 1 the breakup is expected to start in the V-regime before transitioning into the VI-regime. Figure 6(a) shows that when Oh = 10 (curve 5), the breakup speed r min reaches the value predicted in the V-regime, −0.071, at approximately r min = 10 −2 and remains there until r min = 10 −4 , i.e. until the end of the computation. For higher Oh, curves are graphically indistinguishable from the case of Oh = 10, so this is the Stokes flow solution. Therefore, for r min 10 −4 the VI-regime is not observed at large Oh. Figure 7 shows the region of Oh-r min phase space where V-regime dynamics is encountered, with the circles marking the computationally determined boundaries of this regime. The procedure for identifying this regime is shown in the Appendix and, as discussed in § 5, it is based on the speed of breakupṙ min remaining between −0.071 ± 0.015. Next, the changes in behaviour triggering the entrance and exit from this V-regime are considered.

Entrance into the V-regime
The V-regime is entered (r →V min ) at a constant Oh-independent value r min ≈ 10 −2 , which suggests that this transition may occur when the thread can be considered slender, so that the assumptions behind the similarity solution (2.2) become valid. To determine whether the geometry near the breakup point is 'slender' or not, we plot L r /L z in figure 8. For Oh = 10 (curve 5), and generally for larger Oh, initially L r /L z ∼ 1 (no scale separation), but as r min shrinks the breakup region becomes slender. Defining the region to be slender when L r /L z < 0.1, it is found that for all Oh > 1 slenderness occurs at r min ≈ 10 −2 , in agreement with our result for r →V min . Therefore, the transition into the V-regime appears to be dictated by the geometry of the thread.

Exit from the V-regime
The V-regime is left when the speed of breakupṙ min diverges from −0.071 ± 0.015 (e.g. curve 1 in figure 6a). Figure 7 shows that this transition follows the proposed scaling for V→VI that r min = A Oh −3.1 , with computations finding A = 5.5 × 10 −4 . . (Colour online) Evolution of the slenderness of the free surface in the breakup region for 1: Oh = 10 −3 , 2: Oh = 10 −2 , 3: Oh = 0.14, 4: Oh = 0.16, 5: Oh = 10. This quantity, which must be small for the thread to be slender (as is required by the V-and VI-regimes), is defined using L r = r min and L z = |z(r = r min ) − z(w = w max )|.
The transition out of the V-regime (r V→ min ) occurs when inertial effects become nonnegligible, so that (2.2) is no longer valid. Figure 6(b) shows that the increase in Re local follows remarkably well the scaling predicted by the similarity solution for the V-regime that Re local ∼ r −0.65 min . This also serves to validate our definition of Re local . Circles placed at r V→ min on curves in figure 6(b) show that the transition out of the V-regime occurs when Re local ≈ 0.85. One may expect that at this point, where Re local ∼ 1, the VI-regime will be encountered, and this will be the focus of the next section.

Viscous-inertial balance (intermediate Oh)
At Oh = 0.16 in figure 9(a,b) one can see the appearance of a satellite drop centred at z = 0. The close-up images show how this is driven by symmetry breaking about the pinch point after r min ≈ 10 −2 (curve 2), as observed experimentally in Rothert, Richter & Rehberg (2001). Consequently, two pinch points occur at a finite distance either side of the symmetry plane (z = 0) with a large pressure acting to push fluid out of the thinnest region (figure 9d). In contrast to profiles in the I-regime (cf. § 9), the geometry near the pinch point remains slender; see curve 4 in figure 9(b) and the slenderness data in figure 8 (curve 4).
Features that indicate the breakup is occurring in the VI-regime are as follows. Figure 4(a) shows that the free-surface profile at r min = 10 −4 compares well with the similarity solution (dashed line) derived in Eggers (1993) with no free parameters. Below the pinch point (z < 0.25) in figure 4, deviations of the computed free-surface shape from the similarity solution occur as there is no 'far field' in our set-up (as the plane of symmetry is at z = 0). This also affects the velocity profiles (see figure 9c), which agree with the similarity profiles in Eggers (1993) above the pinch point but not below -a point that will be discussed further in § 11.2.1. However, curve 2 in figure 5(a) shows that the maximum axial velocity follows w max = 0.3Oh r −0.5 min from (2.3), again with no adjustable parameters.
In figure 10(a), the breakup speed is shown for Oh > 0.15, values whose significance will become apparent. In this range, all curves tend towards the value of −0.030, indicating that breakup enters the VI-regime.
What is immediately striking about the curves in figure 10(a) is that they appear to oscillate around the value of −0.030, converging towards it as r min → 0. To highlight this behaviour, figure 10(b) compares curves for Oh = 10, which remains in the V-regime, and Oh = 0.16, where the VI-regime is most prominent. Whilst the value 10 -1 10 -2 10 -3 10 -4 Eggers (1993) Eggers (1993) Papageorgiou ( of −0.071 is approached monotonically from below by curve 2, the approach of curve 1 towards −0.030 is entirely different. The motion in the VI-regime is clearly oscillatory in the radial direction, and on closer inspection of figure 5(a) and figure 11 oscillations in axial quantities w max (and hence Re local in figure 6b) and z(r = r min ) can also be seen. Having identified this behaviour, it can also be recognised in the plot of r min against τ in curve 2 of figure 5(b), with the bridge radius oscillating around and converging towards the dashed line predicted by the similarity solution (2.3). As far as we are aware, this is the first time this oscillatory behaviour has been observed and it will be discussed further in § 11.2.  Figure 12(a) focuses on the narrow range Oh = 0.14 − 0.16 to highlight a sharp transition in the dynamics of breakup which occurs at a critical Oh c ≈ 0.15. For Oh = 0.16 > Oh c (curve 2), the transition into the VI-regime occurs at r min ≈ 10 −2 , whilst for Oh = 0.14 < Oh c (curve 1) the transition into this regime is delayed until r min ≈ 4 × 10 −4 . Remarkably, we will show that this is due to the existence of a 'low-Oh V-regime', first discovered in Castrejón-Pita et al. (2015), which precedes entry into the VI-regime for Oh < Oh c . In figure 12(a) for Oh = 0.14 (curve 1), it can be seen that this recently discovered regime is entered when r →V min = 3 × 10 −3 and exited again when r V→ min = 6 × 10 −4 , during which period the speed of breakup is in the V-regime range of −0.071 ± 0.015.
The transition in flow behaviour at Oh c is also seen from the free-surface shapes in the vicinity of the breakup region. In figure 13(b) significant differences in the breakup geometry are shown for Ohnesorge numbers just below (Oh = 0.14) and above (Oh = 0.16) the critical value. At r min = 10 −4 , for Oh = 0.16 (curve 1a) a long slender thread of length ≈0.22 connects the hemispherical volume above (z > 0.27) to a small satellite drop below (z < 0.05), whilst at Oh = 0.14 (curve 1) the thread is much shorter at ≈0.07. Consequently, the satellite drops formed in each case differ considerably, with Oh = 0.16 forming a much fatter drop connected to a longer thin thread.
The free-surface profile in the low-Oh V-regime (curve 2) does not have the slender symmetric shape in the breakup region that is characteristic of the conventional Vregime, which suggests that this regime is encountered in a weaker sense than the high-Oh V-regime which satisfies all expected characteristics (see § 7). Using more stringent criteria for how the regimes are defined could result in this becoming a region of phase space where none of the similarity solutions are accurate; however, our definition is based on the breakup speed alone and this identifies it as a low-Oh V-regime.

Identification of regimes
In figure 14 the VI-and low-Oh V-regimes are shown on a phase diagram. The extreme change in behaviour at Oh c = 0.15 is clearly visible.  ) for Oh > Oh c = 0.15. For Oh < Oh c the low-Oh V-regime appears and is bounded by curves 2: r min = 0.2Oh 2 and 3: r min = 10Oh 2 as r min → 0.

Entrance into the VI-regime for Oh > Oh c
For Oh > Oh c , the transition into the VI-regime follows the expected scaling for V→VI of r min = B Oh −3.1 , with B = 2.3 × 10 −4 . Curves 1-3 in figure 6(b) show that as expected Re local increases until it reaches ≈1, a value characteristic of the VI-regime. The sharp changes in the gradients of these curves at the point where Re local is a maximum marks the beginning of the formation of satellite drops, i.e. the axial shift of the minimum bridge radius to z(r = r min ) > 0, as can be seen from figure 11. This signals the appearance of satellite drops and entry into the VI-regime. Although L r /L z increases at this point (figure 8), the thread remains slender. Once in the VI-regime, we do not see any evidence of an exit from this regime.

Entrance into the low-Oh V-regime for Oh < Oh c
When Oh < Oh c , appearance of the VI-regime is delayed by the presence of a low-Oh V-regime. The entrance to this new regime follows a r min ∼ Oh 2 scaling as r min → 0, which is characteristic of the I→VI transition, although here, as shown in § 9, this is an I→V transition. Figure 12(b) corroborates the argument for a low-Oh V-regime by showing, counterintuitively, that once Oh < Oh c lower values of Re local are recovered. For example, for Oh = 0.16 there is a dip to Re local = 0.59, whilst for Oh = 0.14 it falls as low as Re local = 0.18. For Oh = 0.14, there is then a substantial period (4 × 10 −4 < r min < 7 × 10 −3 ) during which Re local remains at a value characteristic of the V-regime (Re local < 0.85). Identifying the V-regime from the breakup speed in figure 12(a) results in a slightly smaller period (6 × 10 −4 < r min < 3 × 10 −3 ), a mismatch which appears to be caused by an initial lack of slenderness in the thread (figure 8), with L r /L z < 0.1 only once r min < 9 × 10 −4 . This may explain why althoughṙ min ≈ −0.071 in the new regime, there is no asymptotic approach to this value but rather a transient passing through speeds associated with a V-regime.
Entrance into the low-Oh V-regime is triggered by the minimum radius moving away from z = 0 as a satellite drop is formed. This is accompanied by other features that are characteristic of an I-regime, namely a corner-like free-surface shape (curve 3 in figure 13a) and a decrease in the axial length scale L z due to the position of maximum velocity approaching the pinch point. However, in the newly created breakup region the axial velocity U z is still relatively small, so that U z L z (i.e. Re local ) is not large enough to trigger I-regime dynamics and instead a V-regime is encountered, a feature discovered in Castrejón-Pita et al. (2015). What remains unclear is why this mechanism persists at smaller Oh and prevents transitions from the I-regime directly into the VI-regime.

Entrance into the VI-regime for Oh < Oh c
As can be seen from figure 14, the transition out of the low-Oh V-regime is quickly followed by a transition into the VI-regime for Oh < Oh c . This boundary appears to follow a r min ∼ Oh 2 scaling as r min → 0, although for Oh > 4 × 10 −2 it is approximately constant. The ∼Oh 2 scaling was predicted for the I→VI transition and observed above for the entrance into the low-Oh V-regime. However, why this scaling is followed for V→VI is unclear.

Inertia-dominated flow (small Oh)
In figure 15 flow profiles for Oh = 10 −3 are shown. As seen in § 8, a satellite drop is formed, but now the free surface forms a corner at the pinch point (figure 15b). In a narrow region nearby a rapid increase in w fs and p fs (figure 15c,d) can be seen. Notably, w is no longer r-independent. Below the pinch point the free surface forms an angle close to 18.1 • (predicted for the I-regime in Day et al. (1998)) with the zaxis (dashed line), but there is no agreement with the predicted angle above (112.8 • ), which is found to be 78 • . Therefore, no 'overturning' of the free surface is observed, as discussed further in § 11.3.
The maximum axial velocity is shown in figure 5(a) to scale as ∼r −1/2 min as predicted by the similarity solution (2.1). In contrast to the other two regimes, the radial velocity −ṙ min scales in the same way and is close to w max when multiplied by a factor of 5, showing that in the I-regime there is no separation of scales in the r and z directions. This is supported by observations that the thread is not slender as pinch-off is approached, as can be seen from curve 1 in figure 8.  Day et al. (1998)), (c) the axial velocity at the free surface, w fs , and (d) the pressure at the free surface, p fs , at 1: r min = 10 −1 (blue), 2: r min = 10 −2 (red), 3: r min = 10 −3 (green) and 4: r min = 10 −4 (brown).
Recalling that regions wherel min is approximately constant result in a r min ∼ τ 2/3 scaling indicative of the I-regime, curve 1 in figure 16(a) shows the appearance of this regime for Oh = 10 −3 around 10 −3 < r min < 10 −2 . Computations show thatl min = −A 3/2 I ≈ −0.5 so that A I = 0.63, which can be used to extract the breakup time and plot r min against τ in curve 1 of figure 5(b). The value of A I is close to that given in Eggers & Villermaux (2008) of 0.7.

Identification of the I-regime
Usingl min to define the I-regime, and starting from the point at which the satellite drop is formed, i.e. r min < 3 × 10 −2 (figure 11), where the geometry of the I-regime starts to take shape, the boundaries of the I-regime are shown in figure 17.

Entrance into the I-regime
Similar to the V-regime, transitions into the I-regime r →I min are due to geometry, with a certain time required for the thread to form the corners which are characteristic of this regime. As one can see from curves 1 and 2 in figure 8, the breakup region is not slender and L r /L z ≈ 0.6 remains approximately constant as r min → 0.

Exit from the I-regime
The exit from the I-regime r I→ min occurs when viscous effects become significant and the local Reynolds number drops below a critical value. Figure 16(b) shows that the reductions in Re local follow the ∼r 0.5 min scaling predicted for the I-regime (2.1), which confirms that our definition of Re local is a good one. It is found that the I-regime is exited when Re local ≈ 14. From § 8 we know that at this point the breakup will transition into the low-Oh V-regime, and figure 17 shows that this entire boundary follows an ∼Oh 2 scaling.

A phase diagram for breakup
Having computed transitions into and out of the three regimes across parameter space, the results can be stitched together to produce the phase diagram shown in figure 18. Notably, once a regime is reached, which is around r min ≈ 10 −2 across Oh, the transitions between regimes are relatively sharp, so that there are not vast regions of parameter space where none of the scaling laws are applicable. There are also no overlaps between the regimes, partially due to our choices of how to define being 'in' each regime. Consequently, rather than worry about transitions into and out of the different regimes, it is now sensible to talk about transitions between different regimes. This is how figure 18 was converted to produce the neat phase diagram figure 2(b) in § 2.2.3, where (i) transitions between the regimes are provided and (ii) only the asymptotic form of these transitions as r min → 0 are kept (i.e. using only the dashed lines in figure 18). It is of course possible that there is further complexity in the phase diagram below r min = 10 −4 . Here, analytical work is required as computational studies are inherently confined to a finite resolution below which the breakup behaviour can only be implied. The most surprising results are the appearance of a low-Oh V-regime which prevents I→VI transitions and the dramatic change in behaviour around a critical value Oh c = 0.15. For Oh > Oh c the scaling for the V→VI transition follows that expected from theory (r V→VI min ∼ Oh −3.1 ), whilst both the transitions from the I-regime into the low-Oh V-regime and the exit into the VI-regime follow the r min ∼ Oh 2 scaling predicted for the I→VI transition, but only as r min → 0. Notably, it was shown that many of the transitions between regimes can be predicted from our measures Re local and L r /L z .
An experimentally verifiable prediction for the existence of a sharp transition at Oh c = 0.15 is that going from Oh = 0.16 down to Oh = 0.14 decreases the length of the thin liquid thread which connects the satellite drop to the main volume by a factor of three; see figure 13(b).

Discussion
Computations using a multiscale finite element method have allowed us to accurately construct a phase diagram for breakup and determine the transitions between different regimes. During the analysis, a number of features have been encountered which are worthy of further attention, and these are now considered.

Unexpected transitional behaviour
The possibility of multiple flow transitions was suggested in Eggers (2005) and first observed in Castrejón-Pita et al. (2015). The low-Oh V-regime can be seen in figure 2 52 Y. Li and J. E. Sprittles of Castrejón-Pita et al. (2015) at similar values of Oh to those found here. However, in contrast to Castrejón-Pita et al. (2015), no evidence of a transient high-Oh I-regime (their figure 4 for the case of Oh = 1.81) has been observed. This could be because our method of defining the I-regime does not pick up transient regimes far from breakup; however, as can be seen from figure 18, our phase diagram is full for r min < 10 −2 so that any I-regime that has been missed will not appear in the later stages of breakup. Moreover, for all cases where a high-Oh I-regime could exist, Re local remained at values well below that seen for the I-regime, where Re local > 14, giving strong evidence for the absence of a high-Oh I-regime.
To make unambiguous conclusions about the dominant forces as breakup is approached, one could directly extract the size of the inertial and viscous forces from the computed terms in the Navier-Stokes equations and compare their magnitudes in the breakup region. This would be a more advanced method for calculating Re local that may provide additional insight into the transitions between different regimes. However, we have found such measurements difficult to obtain in the breakup region due to the tiny elements found there and the singular nature of the dynamics. Such issues are amplified when determining second derivatives of the velocity, which itself is only approximated quadratically across each element, in order to calculate the viscous forces. This will be the focus of some of our future work in this area.
An intriguing possibility is that the low-Oh V-regime is responsible for the transitional behaviour observed in Rothert et al. (2003), where water-glycerol mixtures in the approximate range 0.06 < Oh < 3 were shown to transition V→VI at a constant r min rather than following the predicted scaling r min ∼ Oh −3.1 . Given the relatively low Oh considered in these experiments, it seems possible that it was the low-Oh V-regime that was being encountered rather than the usual one. Furthermore, figure 18 shows that in the range 0.04 < Oh < 0.14 the transition from the low-Oh V-regime into the VI-regime (r V→VI min ) is approximately constant. Unfortunately, no direct comparison to Rothert et al. (2003) is possible due to the different flow configurations used, so at present we can only speculate about the possible observation of a low-Oh V-regime in their experiments. 11.2. Oscillatory convergence to the VI-regime Figure 10(b) clearly illustrated that convergence towards the V-and VI-regimes is fundamentally different: at Oh = 10 the bridge speed increases monotonically towards the value predicted by the similarity solution in the V-regime of −0.071, whilst for Oh = 0.16 the bridge speed converges in an oscillatory manner towards the speed predicted by the similarity solution in the VI-regime of −0.030 (preliminary results for the case where the plates confining the liquid bridge move apart are similar except the VI-regime is entered earlier, i.e. at larger r min ). The results are reinforced by figure 19, which shows the same oscillations whenṙ min is plotted against the time from breakup τ . Notably, the wavelength of the oscillations is approximately constant in logarithmic time ln τ , a natural choice of variable, used e.g. in Eggers (1997), in the derivation of similarity solutions for breakup. Our results suggest that the similarity solutions in both the V-and VI-regimes are stable to perturbations, but that in the former case the associated eigenvalues of the perturbation are negative and purely real whilst in the latter case they are complex with negative real part.
For the V-regime our findings agree with those in Eggers (2012), where it is shown that out of an infinite hierarchy of similarity solutions, only (2.2) is stable with purely real eigenvalues. Our findings for the VI-regime are entirely new and have not been predicted computationally, analytically or experimentally. The likely reasons for this are as follows: (1) Computations. Very few simulations have captured the small scales resolved in our work for the axisymmetric Navier-Stokes system. However, in cases where they have, e.g. in Castrejón-Pita et al. (2015), the focus has been on r min against τ , whereas the oscillations are most clear when plottingṙ min against τ or r min . What is surprising is that this behaviour has not been observed in slender jet codes which recover the similarity solution in Eggers (1993). Our own preliminary simulations indicate that the same oscillations can be observed in this setting too.
(2) Analysis. The current situation is summarised in p. 39 of Eggers & Villermaux (2008): thus far no one has analytically calculated the eigenvalues of the Navier-Stokes system for breakup. Consequently, it is unknown (i) whether the similarity solution is stable and (ii) if the associated eigenvalues are purely real or complex. The framework for the stability analysis is provided in appendix B of Brenner et al. (1996) and progress on similar flows has been achieved in Bernoff, Bertozzi & Witelski (1998), but extending these results to the Navier-Stokes pinch-off remains an open problem. (3) Experiments. As with computations, so far the focus has been on r min against τ , where the oscillations can be missed, so that converting some of this data intoṙ min against τ or r min would be useful. However, here one faces the additional complication of noisy data, which could make the accurate extraction of derivatives a tricky procedure. Despite this, the oscillations observed are quite significant; at Oh = 0.16, between r min = 10 −3 -10 −2 the bridge speed doubles, from −0.02 to −0.04 ( figure 10b). If the length scale is R = 1 mm this occurs when the dimensional bridge radius is around 1-10 µm, i.e. near the limits of optical resolution. Therefore, there is hope that experimental data can identify these changes.
11.2.1. Bumps on the VI similarity solution Interestingly, continuing some of our simulations in the VI-regime beyond our self-imposed minimum radius of r min = 10 −4 resulted in the growth of a number of interesting features which, as far as we can see, appear unrelated to the aforementioned oscillations. These should be treated with caution, as the accuracy of the numerical scheme is close to its limit at such small radii, but the features appear to be robust. As can be seen from figure 20(a) for Oh = 0.16, at around r min = 5 × 10 −5 the free surface has developed waves close to the pinch point. These feature in plots of axial velocity in figure 20(b), where one can see that the maximum absolute velocity (which will have w < 0) shifts from around z = 0.08 (in curve 1) to z = 0.23 (curve 4). The velocity profile in curve 4, with the maximum absolute velocity near the pinch point, is closer to that predicted by the similarity solution in Eggers (1993).
It is possible these bumps could be related to the iterated instabilities observed experimentally in Brenner, Shi & Nagel (1994) and interpreted as successive instabilities of the similarity solution of Eggers (1993). In this work, the instability was attributed to thermal fluctuations which were added to the lubrication formulation. However, it was pointed out that because this formulation only approximates the full Navier-Stokes system, such external noise may be unnecessary. The iterated instabilities have been seen computationally in McGough & Basaran (2006) for a surfactant-covered jet, but our results suggest that the instabilities can be generated without requiring additional effects such as noise or surfactants.
There could also be a relation to the destabilisation of similarity solutions due to finite-amplitude perturbations observed in Brenner et al. (1996), with intriguing similarities between our figure 20(a) and their figure 5 for the evolution of an unstable similarity solution.
11.3. Overturning of the free surface in the I-regime Figure 15 showed that at Oh = 10 −3 no overturning of the free surface was observed and, as far as we are aware, this does not contradict any previous findings in the liquid bridge geometry. Although simulations in Suryo & Basaran (2006) show overturning in the same geometry, this phenomenon only occurs for non-Newtonian cases. In contrast, for drop formation from an orifice, experimental studies show cases that do (Castrejón-Pita et al. 2012) and do not (Brenner et al. 1997 whilst simulations using the inviscid equations do show this feature (Schulkes 1994) and those using slender jet theory are unable to (Brenner et al. 1997).
The most comprehensive computational study of overturning is given in Wilkes et al. (1999), where the finite element method is used to study the effect of Oh and the Bond number Bo (i.e. gravitational forces) on this phenomenon. It is shown that the critical Oh at which overturning is observed and the angles which the free-surface shape makes with the z-axis at the pinch point depend on Bo, which alters the geometry of the breakup. Furthermore, the angles recovered vary, with a maximum angle of 95 • observed, in contrast to the value of 112.8 • predicted from the inviscid theory (Day et al. 1998).
In the drop formation geometry, threads of arbitrary large lengths can be formed, whereas in our set-up the liquid is confined between two stationary plates. This appears to be the reason why no overturning has been observed even in cases of small Oh which fall below the overturning limit in Wilkes et al. (1999). In particular, rather than a thread connected to a free drop, our geometry has a short thread (the satellite drop) connected to a hemispherical volume pinned to a plate.

Implications for breakup in CFD
Computational approaches that 'go through' the topological change of breakup all rely on cutting the thread, either manually (as in finite element methods) or automatically (as in volume-of-fluid approaches -VoF) once it reaches a specific r min , which may be implicitly related to the grid size (as in VoF). Therefore, it is of interest to know how accurate this approach is and what could be lost from a premature truncation. To do so, for the sake of argument consider a scheme with a fixed resolution (mesh size) of h = 5 × 10 −3 , giving a (relatively large) grid of 200 × 200 for our problem, which cuts the thread wherever r min < h.
The most obvious feature that can be lost from the cutoff is the satellite drop, as seen in Fawehinmi et al. (2005). Looking at figure 11, satellite drop formation, indicated by z(r = r min ) becoming non-zero, can clearly be seen for both Oh = 0.16 (curve 2) and Oh = 0.5 (curve 3). This occurs when the symmetric V-regime transitions into the asymmetric VI one, as can be seen clearly from free-surface profiles at Oh = 0.16 in figure 9. For Oh = 0.16 the pinch point moves when r min ≈ 10 −2 , whilst for Oh = 0.5 this does not happen until r min ≈ 10 −3 . Therefore, using our cutoff, at Oh = 0.16 the thread would be cut near to the correct pinch point, whilst for Oh = 0.5 the thread would be severed at z = 0 so that the formation of a satellite drop would be missed. Most worryingly, in the latter case the post-breakup state would be entirely wrong, as instead of having two breakup points a distance z = 0.35 apart (giving three distinct volumes) there would be just one at the centre of the drops (and hence two volumes). Notably, for R = 1 mm the minimum radius r min = 10 −3 corresponds to a dimensional radius of 1 µm, so that these are truly macroscopic quantities.
To capture the features of breakup, such as satellite drops, one has to either (i) develop codes with huge levels of resolution/refinement, which will be extremely costly (particularly in 3-D) due to the multiscale nature of this phenomenon, or (ii) intelligently utilise the similarity solutions to take the codes through the topological change and up to a suitable post-breakup state. Such a scheme would be complex, but could lead to increased accuracy at a lower computational cost. By identifying in which regions of parameter space the different similarity solutions are accurate and devising methods to analyse in which flow regime a given breakup is occurring 56 Y. Li and J. E. Sprittles we have taken just some of the first steps in the development of such a scheme. Forthcoming works will extend these results to cases in which there is a strong externally driven elongation of the thread and/or interface formation physics which creates a singularity-free breakup; see Shikhmurzaev (2007).
In figure 21, the calculation of the boundaries of the viscous regime is shown for the case of Oh = 1. The V-regime is defined by speeds −0.086 <ṙ min < −0.056 (dashed lines). The transition into this regime is found to occur when r →V min = 1.4 × 10 −2 and the exit from this regime is when r →V min = 6.1 × 10 −4 . Once these values have been calculated they are placed on the phase diagram, as shown in figure 7 (as circles), to define the V-regime.