On the physics of transient ejection from bubble bursting

Abstract Using a dynamical scaling analysis of the flow variables and their evolution due to bubble bursting, here we predict the size and speed of ejected droplets for the whole range of experimental Ohnesorge and Bond numbers where ejection occurs. The transient ejection, which requires the backfire of a vortex ring inside the liquid to preserve physical symmetry, shows a delicate balance between inertia, surface tension and viscous forces around a critical Ohnesorge number, akin to an apparent singularity. Like in other natural phenomena, this balance makes the process extremely sensitive to initial conditions. Our model generalizes or displaces other recently proposed ones, impacting on, for instance, the statistical description of sea spray.


Introduction
Everyday experience teaches that radially convergent flows close to a liquid surface produce vigorous transient liquid ejections in the form of a jet perpendicular to the surface, as those seen after bubble bursting (Kientzler et al. 1954), droplet impact on a liquid pool (Yarin 2006) or cavity collapse (Ismail et al. 2018). The mechanical energy of the flow comes from the free surface energy of the initial cavity. Among these processes, bubble bursting can be considered the parametrically simplest and most ubiquitous one: just two non-dimensional parameters, namely the Ohnesorge (Oh) and Bond (Bo) numbers, determine the physics and the outcome (Oh and Bo, weighting the viscous and gravity forces with surface tension, respectively). It has received special attention in the scientific literature due to its impact at global planetary scales (bubble bursting at the sea surface, e.g. Kientzler et al. (1954), Blanchard & Woodcock (1957), Veron (2015) and Sampath et al. (2019); disease and pandemic transmission, e.g. Bourouiba (2021)) and in more 'indulging' applications (e.g. Séon & Liger-Belair 2017).
Current works state that the subject is already deeply understood (Berny et al. 2020;Sanjay, Lohse & Jalaal 2021), yet despite optimistic statements, even in the simplest case where the liquid properties are constant and the gas-to-liquid density and viscosity ratios are very small, unsettling but crucial questions remain open. Within an ample range of Oh and Bo, bubble bursting exhibits a strong focusing effect due to the nearly cylindrical collapse of a main wave at the bottom of the parent bubble. In these cases, a tiny bubble gets trapped and an extremely thin and rapid spout emerges. That initial spout grows into a much taller, thicker and slower transient jet that may eventually expel droplets. Their eventual size and speed is determined by Oh and Bo. However, for an apparently critical Oh, (Oh * ), the initial ballistic spout keeps its large speed and extremely thin radius until it ejects a nearly invisible droplet, such that the spout momentum (proportional to its velocity times its volume) seems to vanish. Strikingly, a nearly symmetric behaviour is observed around that Oh * (Séon & Liger-Belair (2017), figure 16). Around that special parametrical region, experiments and numerical simulations exhibit a strong dispersion attributed to local viscous effects and interaction with the gas environment (Brasz et al. 2018;Dasouqi, Yeom & Murphy 2021).
In general, the global streamlines of the flow resemble a dipole with its axis perpendicular to the interface (figure 1), which is responsible for the vigorous liquid-air transfer process characteristic of this phenomenon (Lee et al. 2011). The energy excess of the convergent flow leads to both a fast transient capillary jet and a vortex ring in opposite directions but equivalent effective momenta by virtue of Newton's third law of motion. Their radically different kinematics, due to a large density disparity across the interface, is not an obstacle to their profound symmetry, as we will show here. We claim that this symmetry and its scalings are the keys to determine the eventual ejected droplet size and speed, and the appearance of singular dynamics. In this regard, some of the above cited authors and many others (e.g. MacIntyre 1972; Duchemin et al. 2002;Thoroddsen et al. 2018) have noticed the trapping of tiny bubbles after bursting. In particular, the one trapped at the bottom of the formed cavity is the relic of a singularity produced by the collapse of a pilot capillary wave from the rim formed after the upper film rupture (Duchemin et al. 2002;Thoroddsen et al. 2018). This singularity would explain why the range of Oh numbers where this occurs exhibits the largest ejection speeds and smaller drops.
In this work, we present: (i) a universal non-trivial scaling of the evolution of the flow kinematic variables (characteristic lengths and velocity); (ii) a prediction of the eventual ejected droplet size and speed, for the complete range of Oh and Bo for which ejection occurs (which generalizes Gañán-Calvo (2017)); and (iii) an explanation why an apparently singular value of Oh * appears. The existing data from prior works (e.g. those collected in Gañán-Calvo (2017); Brasz et al. (2018), Deike et al. (2018), Gañán-Calvo (2018) and Berny et al. (2020)) and novel numerical simulations are compared with our theoretical results. When the droplet size and ejection velocity are properly scaled according to our model, the data show a remarkable collapse around our model except in the vicinity of Oh * . A careful inspection reveals that all data from numerical simulations exhibit an enormous degree of dispersion that experimental data do not display. In fact, while numerical series should be independently fitted, the available experimental series appear to obey a single universal fitting to our model. Given that it relies on just two non-dimensional parameters alone, Oh and Bo, this has a twofold consequence: (i) it highlights the critical importance of the artificial initial conditions after film bursting in the numerical simulations among the different authors, which precludes universality; and (ii) it provides a proof that the natural 929 A12-2  Figure 1. General overview of the flow development around the critical time t 0 of collapse of the main pilot wave at the bottom of the cavity, for Oh = 0.032, Bo = 0. The three instants here illustrated are t 0 − t = 3.54 × 10 −3 t c , t 0 − t = 8.4 × 10 −6 t c and t − t 0 = 4 × 10 −4 t c (with t c = (ρR 3 o /σ ) 1/2 ). (a) Global flow streamlines (similar to dipole contours at a liquid-gas surface) showing a particular one (line A-B) ending at the point (B) just above where the collapse eventually occurs. The blue sphere of the three-dimensional rendering is the initial bubble immediately after bursting. (b) Local details of the same instants. The stream function levels shown are closer around the tiny trapped bubble to exhibit the vortex ring. Note that the upper axial point of the ellipsoid defining the vortex ring coincides with the point where the surface collapsed at t = t 0 and remains at that position: observe the horizontal line connecting the three subpanels of panel (b). The main flow velocities W (radial) and V (axial, ejection) are indicated. The dashed lines represent the instantaneous stream tube which feeds the collapsing point or the issuing jet. In these simulations using the volume of fluid method (VOF) ( Basilisk, Popinet 2015), the density and viscosity of the liquid is 1000 and 100 times that of the gas, respectively, to reflect similar relations to the air-salt water ones. process should be fundamentally universal in terms of Oh and Bo alone, as expected from the universality of initial film drainage and rupture (Eggers & Fontelos (2015) and references therein), which supports the universal validity of our model once the three proposed physically relevant fitting parameters are found.

Formulation and dynamical scaling analysis
Consider a gas bubble of initial radius R o very close or tangent to the free surface of a liquid with density, surface tension and viscosity ρ, σ and μ, respectively, in static equilibrium under a relatively small action of gravitational acceleration g normal to the free surface far from the bubble. The flow properties are considered constant. At a certain instant, the thin film at the point of tangency breaks and the process of bubble bursting starts (figure 1a). The problem is thus determined by the Ohnesorge (Oh = μ/(ρσ R o ) 1/2 ) and Bond (Bo = ρgRo 2 /σ ) numbers alone (Montanero & Gañán-Calvo 2020). Alternatively, instead of Oh, some authors prefer to use the Laplace number (La), which is simply La = Oh −2 (e.g. Duchemin et al. 2002;Lai, Eggers & Deike 2018). Also, the Morton number, Mo = Bo Oh 4 , has been used as a relative measure of the gravitational forces (e.g. Berny et al. 2020). The liquid rim that is initially formed pilots two main capillary wavefronts, one that advances along the bubble surface towards its bottom, and the other propagating away from the cavity. In some parametrical regions of the domain (Oh, Bo) of interest, a significant fraction of the energy content in the nonlinear wave spectrum of the initial surface configuration (after the bubble film burst) is in the small-wavelength domain (wavelets). For small Oh, these wavelets may arrive first to the bottom, but they do not produce any significant effect compared with the collapsing main wave (Gañán-Calvo 2018) in the parametrical realm considered here. When the main wavefront approaches the bottom, it becomes steep. If the wavefront leads to a nearly cylindrical collapsing neck in a certain region, a tiny bubble gets trapped below the surface once the ejection process commences. The process of radial collapse and bubble trapping is locally described by the theory of Eggers et al. (2007), Fontelos, Snoeijer & Eggers (2011 and Eggers & Fontelos (2015) when a strong asymmetry in the axial direction occurs. However, the asymmetry of bubble trapping here observed leads to an extremely rapid and thin initial ejection in the opposite direction from the trapped bubble. Figure 1(b) shows three illustrating instants of the flow development around the critical instant t 0 at which the axial speed of the liquid surface on the axis reaches its maximum (immediately after the tiny gas neck collapses). When t < t 0 , the flow is predominantly radial (spherical) with a characteristic speed W at the interface. After collapse (t > t 0 ), while the main flow keeps running radially with speed W, the axial asymmetry of the flow and the kinetic energy excess of the liquid is axially diverted in the two opposite directions, with radically different results: (i) towards the open gas volume, producing the liquid spout with speed V and characteristic radial length R; and (ii) towards the liquid bulk, as a backward reaction vortex ring (figure 1b, in colours) with characteristic length L, around the trapped microbubble. Eventually, the advancing front of the resulting capillary jet expels a droplet or droplet train with characteristic radius R.
Obviously, all scales (W, V, R, L) are time-dependent along their evolution. Some relevant prior works (e.g. Zeff et al. 2000;Lai et al. 2018) have analysed the time self-similarity of the flow variables. Assuming that those time self-similarities exist, our focus is here to obtain closed relationships among those scales and to predict the size and speed of ejections reflected by the eventual values of R and V at the end of the process. This approach, though, demands the identification of the key symmetries appearing in the problem that may lead to an effective problem closure (Gañán-Calvo, Rebollo-Muñoz & Montanero 2013). To this end, a set of relations among the radial and axial characteristic lengths and velocities was formulated in Gañán-Calvo (2017, 2018 on the basis that inertia, surface tension and viscous forces should be comparable very close to the instant of collapse of the free surface. Under the light of detailed numerical simulations, the previously proposed simplified models are put into perspective after a rigorous derivation of the new model. Following a similar rationale to that in Gañán-Calvo (2017), this derivation is made in two steps: (i) we obtain explicit relationships among the four time-dependent variables along the process, verifying their validity numerically; and (ii) we close the problem identifying the critical key factors on the energy balance of the process. The final model is compared with the ample set of experimental and numerical available data, and its validity is extended to the whole range of Oh numbers where ejection takes place.

Instantaneous scalings of kinematic variables along the process
The general momentum equation of the liquid can be written as where v is the velocity vector, subscript t denotes a partial derivative with time, p is the liquid pressure, z the axial coordinate and P a the static gas pressure. Equation (2.1) can be multiplied by the unit vector l tangent to any instantaneous streamline, in particular the streamline flowing through a point A where it meets the free surface to a point B at the vicinity of the point of collapse (see figure 1a). Integrating with respect to the streamline coordinate s from A to B yields since the velocity is negligible and pressure is P a at the point A. Here, n is the unit normal on the liquid surface, and z is the depth of point B respect to A. As a general consideration, the liquid velocities are very small everywhere compared with those at distances L to the collapsing region, which may exhibit a self-similar flow structure (Zeff et al. 2000;Duchemin et al. 2002;Lai et al. 2018). The length scale L also characterizes the inverse of the mean local curvature of the liquid surface around the region of collapse, for any given time t. Thus, L obviously changes with time around the instant of collapse. To look for symmetries around t 0 , let us consider two situations, one for t < t 0 and the other for t > t 0 such that their characteristic length scales L are the same (see figure 1). Then, one may estimate the characteristic values of each term of (2.2) for both t < t 0 and t > t 0 .
In any case, the first two terms of (2.2) are always of the same order for this unsteady problem.
For t < t 0 , both the left-hand integral and the kinetic energy term at B in (2.2) should scale as ρW 2 . The surface tension term at B should be proportional to σ/L, and the gravity term to ρgR o . Finally, the curvature of streamlines suggests that the viscous stresses should scale as μW/L, which fixes the scale of the right-hand side term of (2.2). This is so because boundary layers cannot be thinner than the scale L of the collapsing main wave, as shown by Gañán-Calvo (2018): due to the near-zero tangential viscous stress at the free surface, boundary layers cannot be maintained by other means than the internal inertia. In effect, at the vicinity of the surface, there are no other sources of momentum and no other concentrated sources for viscous diffusion than the internal flow, whose fate is dictated by the free surface geometry and, consequently, its energy content. Since that characteristic length L is fixed by the collapsing main wave, this should be the minimum (if any) characteristic length of viscous diffusion as well. Thus, the overall scaling balance of (2.2) can be expressed as where prefactors α 1 and β 1 should be universal constants reflecting the dominance or subdominance of the corresponding terms. Using the natural length and velocity, μ = μ 2 /(ρσ ) and v μ = σ/μ, respectively, and defining ζ = L/ μ , ω = W/v μ = and 1 = β 1 Oh 2 Bo, this equation can be written as When t > t 0 , the vigorous ejection in the axial direction generates the new characteristic velocity V and radial size R scales of the liquid jet. The physical symmetry around t o should keep the scale of the radial velocity W in the bulk and the spatial scale L (including the axial scale of the jet) the same: see the central subpanel in figure 1(b) for Oh = 0.032, Bo = 0. Even before collapse, the flow develops a stagnant region in the liquid just below the collapsing point. From this point, a vortex ring of characteristic (growing) size L and speed W like those beautifully described by Maxworthy (1972) develops; see figure 4 of his work (dissecting absolute and relative speeds is not trivial in these flows, as Maxworthy showed). By virtue of total energy preservation and assuming that the trapped microbubble is much smaller than the vortex, the mechanical energy of the vortex, ρW 2 L 3 , should mirror that of the vigorous microjet, ρV 2 R 2 L. Furthermore, the streamlines below the jet form a feeding stream tube coming from deep in the liquid domain (surrounding the vortex ring, see figure 1b, right-hand subpanel), which should ultimately originate far away from the bubble at the free surface. In contrast, Gordillo & Rodriguez-Rodriguez (2019) ignores such a vortex and proposes a simplified radial (cylindrical) flow model, which hampers the correct identification of the local flow scales. In summary, one can write where χ = R/ μ and υ = V/v μ . As a secondary consequence, the incoming conical flow raises the conical surface at characteristic speed W.
One can now estimate the scaling of the different terms of (2.2) for a streamline ending at a point B at the surface of the jet (see figure 1a), and write the following balance: where, again, prefactors α 2 and β 2 should be universal constants at the end of the process. Using natural scales and defining 2 = β 2 Oh 2 Bo, (2.6) reads As in Gañán-Calvo (2017), (2.4), (2.5) and (2.7) define three relationships between the time dependent scales (χ, υ, ζ, ω). In effect, an explicit relationship between ω and ζ is easily derived after rewriting (2.4): Also, introducing (2.5) and (2.8) in (2.7), the following relationship between υ and ζ results: (2.9) shows the essential proportionality between both spatial variables along the whole process as predicted, consistently with prior analyses (Zeff et al. 2000). (b) Plot of υ versus ω. Black dashed line (υ ∼ ω) also shows the proportionality between both variables for υ below 20; the deviation above that value reflects the very high initial velocity acquired by the spout front at its inception for Oh around Oh * (focusing effect).
Observe that χ ∼ ζ and υ ∼ ω always hold. Testing the validity of (2.10a-c) involves showing that the two limits (i) and (ii) and their corresponding scaling relationships are indeed reached as the time-dependent characteristic scales (χ, ζ, υ, ω) evolve along the process. To do so, we analyse their evolution from numerical simulation (momentum-conserving VOF, Basilisk code). We have recorded the time-dependent values of the maximum radius of curvature in the meridional plane at the point of maximum radial velocity of the interface, and the radius of curvature and velocity of the jet front on the jet axis to represent the time-dependent characteristic values of L, R and V, respectively, for t > t 0 (no scales of R and V appear before ejection, i.e. for t < t 0 ). Since (2.10a-c) does not provide the time dependencies of the variables, in figure 2(a-c) we plot the non-dimensional evolving values of χ versus ζ , ω versus υ and χ versus υ, respectively, for a detailed range of Oh numbers.
Firstly, figure 2(a) demonstrates the fundamental proportionality between ζ and χ along the whole process, as anticipated. In this regard, the small constant of proportionality between both scales, χ 0.03ζ , should be noticed. The importance of this observation will be apparent after the global assessment of the energy balance of the process in § 2.2. Secondly, the remarkable collapse of the curves of χ versus υ in figure 2(c) reveals the two clear regimes as anticipated: viscosity-dominated (α 1,2 1) when χ 1 (R μ ) and surface tension-dominated (α 1,2 → 0) when χ 1 (R μ ), corresponding to χ ∼ υ −1 and χ ∼ υ −2 , respectively. The robustness of the results is thus confirmed: the formation and pinch-off of a droplet is expected to take place once surface tension becomes important, resulting in a dependency as observed. The relatively small deviations from the theoretical prediction (2.10a-c) for the smallest and largest values of υ in figure 2(c) reflect the vibrations of the droplet about pinch-off and ejection (not considered here) and the very initial evolution just after t = t 0 , respectively.

Energy balance
For a given Bo that fixes the initial static geometry of the bubble, droplet ejection is prevented by the loss of mechanical energy above a certain Oh threshold value. In the absence of gravitational forces (Bo 0), a reported value of that threshold is Oh th = 0.052 (Lee et al. 2011). Below that threshold, the finest and fastest jets and droplets are observed for a certain value Oh * < Oh th . Therefore, to close the problem, a key question which is not resolved by (2.10a-c) is the mechanical energy available to perform the ballistic ejection close to t o . One can express the mechanical energy equation in a sufficiently ample fluid volume Ω(t) around the initial bubble from the instant of bubble bursting (t = 0) up to the point of collapse t o . The fluid volume Ω(t) can be formed by the free surface and a hemisphere around the cavity with a radius approximately twice or three times larger than R o (see figure 3): where τ is the viscous stress and I the identity matrix. The surface of the integral domain is S(t) = D 1 (t) + D 2 , where D 2 is a hemisphere with a fixed radius. We define the capillary velocity as V o = (σ/ρR o ) 1/2 , which is the global expected value of the average velocities along the process.
Observing the geometric configuration of streamlines in figure 1, the most relevant feature is a stream tube which ends up at the collapse point, coming from below (figure 1b -subsequently also shown in figure 6). This stream tube, with velocity W and cross-sectional area proportional to L 2 , would ultimately feed the ejection. The global view of the streamlines' configuration in figure 1(a) reveals that the characteristic length of that stream tube is indeed R o , since its configuration remains practically unaltered along the process sufficiently far from the collapse region. Thus, the kinetic energy term should be proportional to (2.12) Note that the initial kinetic energy is zero. If one neglects viscous and gravity works, expression (2.11) suggests the scaling  Deike et al. 2018). There is not a general consensus around the value of Oh * and the origin of this behaviour. Naturally, Oh * depends on the Bond number through the initial shape of the bubble before bursting, since that shape determines the detailed dynamic path of the process. Reducing that geometry to a sphere, i.e. Bo 0, the simulations presented here and those of Berny et al. (2020) and Deike et al. (2018) agree on Oh * 0.03; in contrast, Gañán-Calvo (2017) suggests Oh * 0.043, in the range that would roughly point to Duchemin's data. Therefore, a detailed analysis of the process of collapse is needed, in particular the effects of microbubble trapping and the onset of ejection at the time singularity represented by t 0 (if such a singularity occurs).
2.2.1. The nonlinear behaviour around Oh * : focusing, bubble trapping and the start of ejection Previous works (Duchemin et al. (2002), Ghabache et al. (2014) and Deike et al. (2018), among others) have postulated that the origin of these ballistic ultrafast jets is a finite-time singularity event at Oh * . In reality, a true singularity followed by an instantaneous ultrafast jet occurs in all instances where a capillary wave collapses nearly cylindrically at the bottom of the cavity, locally resembling the bubble collapse phenomenon described by Eggers & Fontelos (2015). However, given the geometrical asymmetry of that collapse in the axial direction (see figure 4a), its consequences are twofold.
(i) A microbubble gets trapped. The dependence of its size on Oh is summarized in figure 4 for an ample range of Oh numbers. (ii) More importantly, an extremely thin ligament of vanishing radial size and very large speed is ballistically emitted in the opposite axial direction from the tiny trapped bubble, as illustrated in figure 5. As discussed in § 2, a vortex ring (usually much larger than the tiny trapped bubble) with a global mechanical energy equivalent to that of the ejected liquid spout is emitted towards the liquid bulk (see figure 1) to keep the overall neutrality of the initial axial momentum (which can be termed a bazooka effect). The described phenomenon of microbubble entrapment and initial ejection can even occur with smaller collapsing waves before the main wave generates the more energetic dominant ejection that eventually overwhelms the former ones (Duchemin et al. 2002;Deike et al. 2018;Gañán-Calvo 2018). However, the total mechanical energy of that initial singular spout is quite small, and surface tension together with axial dissipation immediately blunts it (figure 5) in a self-similar fashion described by Eggers's theory (Eggers 1993). This blunting is equivalent to the ultrafast recoil of a liquid ligament after pinch-off, as observed from a framework moving axially with a large ballistic convective velocity from the pinching point (see figure 5 for Oh = 0.03). As far as sufficient mechanical energy is left after collapse, the spout will grow in size (scaling as χ ) and decrease in speed (scaling as υ), while keeping the scaling relationships (2.10a-c) previously analysed. This growth continues until a droplet is formed by surface tension at the front of the spout. What happens at Oh * has been termed a focusing effect ) of the global evolution around collapse. However, its physical understanding requires further consideration. Figure 6 shows the cases Oh = 0.026, 0.033 and 0.047 just at the onset of ejection (figure 6a-c) and when the length of the liquid ligament is a fixed fraction of R o (approximately 0.27R o , figure 6d-f ). The stream tube that meets the free surface with zero radial velocity is highlighted as a red thick line. This figure shows that the size of the vortex ring just after collapse is much larger when Oh is around Oh * . A detailed observation of the streamlines reveals that the presence of the trapped cavity is associated with the formation of the vortex ring even before collapse when Oh is around Oh * : note that the vortex streamlines originate at the front of the advancing trapped cavity and flow backward, terminating at the rear of the cavity (see figure 6b). As a consequence of the early formation of that backward vortex ring of maximum size, the energy left for ejection immediately after collapse when Oh = Oh * is the minimum on the whole Oh range. The detailed phenomenon just described does not hamper the overall maintenance of the axial momentum balance. In fact, a vortex ring is always formed (see figure 6d-f ). However, as explained, the total mechanical energy left for raising the ejected liquid column becomes minimized for Oh around Oh * . This is consistent with the fact that the cross-section of the stream tube which feeds the liquid column keeps smaller along the process when Oh is around Oh * (see figure 6e). However, the fact that the total remaining mechanical energy is minimal does not imply that the local mechanical energy per unit volume in the ejected ligament must also be minimal. Quite on the contrary, that energy per unit volume is maximized, since the cross-sectional size and thus the total volume of liquid to be pushed is markedly minimized. Consequently, the initial pilot spout can remain ballistic for a longer time just after collapse. The radial scale of that pilot spout, comparable to that of the front of the stream tube responsible for the ejection feed, should be μ , which is the natural spatial scale selected when inertial, viscous and surface tension forces dominate locally (Eggers 1993). All this is confirmed by our numerical simulations, where we use a spatial precision an order of magnitude smaller than μ to adequately resolve the finest details of the ejection.
When Oh > Oh * , though, viscous dissipation takes over, which causes the widening of the emitted spout. Besides, for very small Oh numbers when no bubble trapping occurs, the sheer curvature reversal of the bubble bottom and the appearance of a stagnation point below the surface guarantee the ejection. This occurs even in the absence of bubble trapping and the initial vigorous pilot emission just described, due to the large overall excess energy remaining in the absence of any significant viscous dissipation.

Final energy balance including viscous dissipation
The observed deviation of VR from its expected (inviscid) value V o R o can be quantified expressing the counterpart of (2.11) after collapse as (2.14) The control volume remains the same as that used for (2.11), considering now the highly deformed free surface due to the ejection and possible bubble trapping. While gravity is simply volumetric, proportional to (ρgR o )R 3 o and reflecting the effect of gravity on the initial bubble shape, i.e.
the viscous terms should be proportional to the dominant one, as where constant α 3 should be a universal constant at the end of the process, like α 1,2 . Note that the average total displacement of fluid particles that this term retains should be proportional to vt c R o , being v ∼ V o and t c = (ρR 3 o /σ ) 1/2 . This is so because the main contribution to the integral comes from the larger side D 2 of the control volume (see figure 3) since the normal viscous stresses on the free surface are at most comparable to the surface tension, which is maximum at the front of the growing spout. Moreover, viscosity is, on the one hand, an energy sink that lowers the energy available for ejection, and on the other, it modifies the interface shape at the instant of collapse, favouring or preventing bubble trapping and overall focusing of the flow as described previously in § 2.2.1.
Finally, the surface tension term can be estimated as (2.17) Given that the integral in (2.17) is in reality the difference between the available surface energy between the instant of collapse and any time after that (in particular, at the end of ejection), one should expect the constant k to be quite small, since the main energy conversion from the initially available surface energy to kinetic energy has already been performed from t = 0 to t 0 , implying To this end, one should bear two important considerations in mind: (i) figure 2(b) shows that υ/ω O(1) at the points where these variables have been evaluated, suggesting υ to be comparable in size to ω as expected from basic conservation of momentum at the collapse; and (ii) χ/ζ 0.03. Interestingly, this number is rather comparable to the experimentally observed Oh * values. Therefore, from (2.12) and (2.17) one may expect υ 2 χ 2 /(ω 2 ζ 2 ) k ∼ Oh * 2 , which is justified by the relative size between the initial bubble and the ejected droplet.

Problem closure
While μV/R stays comparable to ρV 2 at the start of ejection, viscous forces decelerate the spout until σ/R becomes strong enough to form a droplet at its tip. Hence, the size of that minimum drop should reflect the ultimate balance between inertia, viscosity and surface tension forces: Thus, a minimum surface tension energy proportional to (σ/ μ ) 2 μ × R o should necessarily remain available along the process for the emission and formation of that droplet, consistent with experimental observations. In conclusion, keeping that in mind together with (2.17), one can formulate the following scaling of (2.11) at the end of the ejection, when the droplet of radius R forms: where prefactors Oh * , α 3 and β 3 should be universal positive constants for very small gas-to-liquid density and viscosity ratios. Equation (2.19) is equivalent to the energy equation in Gañán-Calvo (2017), with the exception of the minimum contribution of the surface tension at the natural limit scale μ . Dividing (2.19) by ρR o μ v μ = μ 2 R o /ρ one gets Oh . (2.20)

Scalings of the ejected droplet size and speed
The right-hand side of (2.20) is a quadratic expression for Oh with a minimum at Oh = 2Oh * 2 /α 3 . Observe that this minimum is located exactly at Oh = Oh * if α 3 = 2Oh * . Thus, if one defines this parameter would measure both the deviation of the location of the minimum from Oh = Oh * and the value of the energy left for ejection. In effect, using definition (2.21), (2.20) reads (2.22) From (2.22) the value of Oh where R is minimum and V maximum is given by around which the dependence is nearly symmetric. When Bo = 0, if α 3 is sufficiently small, Oh c Oh * . Therefore, α 3 constitutes not only a fitting parameter but also a key indicator of the degree of focusing that a certain set of experiments (either physical or numerical) achieve around Oh * , as shown in the next section. As we will see next, experiments reveal that α 3 Oh * , confirming our hypothesis. This fundamental finding explains and justifies the experimental observations by many authors (Duchemin et al. (2002), Walls et al. (2015), , Séon & Liger-Belair (2017), Brasz et al. (2018), Deike et al. (2018) and Berny et al. (2020), among others), but not only this: experiments confirm the continuous validity of (2.22) for the whole Oh range where droplet ejection occurs.
In addition, one can assume α 1,2 = β 1,2 = 0 for simplicity since (i) α 1,2 → 0 is a very good approximation, consistent with the fact that surface tension should become dominant when droplet formation and ejection occurs, as previously demonstrated, and (ii) the role of gravity in the expressions (2.4) and (2.7) is expected to be very small compared with inertia and surface tension forces. This is subsequently confirmed by experiments. With these justified simplifications, (2.9) simply reduces to φ = 1. This finally reduces the scalings (2.8)-(2.10a-c) together with (2.22) to explicit expressions in terms of Oh and Bo as follows: For convenience, one may write (2.24) in terms of the initial bubble radius R o and the capillary velocity V o as and k r,v are prefactors that will be obtained from the experimental fittings.

Experimental verification
To verify our model and obtain the relevant scaling constants for the whole Oh range experimentally explored, 600 experimental and numerical measurements of first ejected droplets and approximately 100 of their corresponding initial velocities have been collected from the literature (Garner, Ellis & Lacey 1954;Hayami & Toba 1958;MacIntyre 1972;Tedesco & Blanchard 1979;Blanchard 1989;Sakai 1989;Spiel 1995;Duchemin et al. 2002;Séon & Liger-Belair 2017;Brasz et al. 2018;Deike et al. 2018;Berny et al. 2020). In addition, to appraise the predictive value of (2.25a,b) in the range Oh ∈ (0.026, 0.052), where the maximum variation among the published results is found, additional extensive numerical simulations using a momentum-conserving VOF scheme (Basilisk) have been performed.
2.4.1. General verification: the constants (Oh * , α 3 , β 3 ) Numerical simulations have been made with gas-to-liquid viscosity and density ratios equal to 0.01 and 0.001, respectively, up to a minimum cell size 1.5 × 10 −4 R o (approximately 0.17 μ ) around the point of collapse and ejection, equivalent to approximately 2 nm for a gas bubble in water for Oh = Oh * = 0.03.
For a general verification of our model, the best first option is to obtain the values (Oh * , α 3 , β 3 ) that minimize the normalized standard deviation of the non-dimensional first ejected droplet radius R/R o from all available experimental and numerical data sets, divided by Ψ , according to the scalings (2.25a,b). This shows a very robust minimum of 24 % and a prefactor k r = 0.18Oh * −2 for Oh * = 0.0288, α 3 = 0.00308 and β 3 = −0.000205 (best fitting), with α 1 = α 2 = β 1 = β 2 = 0. The sensitivity of these results to small variations around α 1 = α 2 = β 1 = β 2 = 0 is rightly inappreciable, as expected, which confirms the negligible role of gravity and viscous forces compared with inertia and surface tension forces when the droplet is ejected. We can use these optimal values of (Oh * , α 3 , β 3 ) in a representable two-dimensional plot of the data against the model to test its goodness of fit.
Since the data show a stronger variability with Oh, this is the parameter usually employed as the abscissa in this kind of study. Thus, to produce a two-dimensional plot, Spiel 1995 Blanchard 1989 Garner 1954 Hayami 1958 Tedesco 1979 Oh * = 0.03, α 3 = 10 -2.9 Oh * = 0.03, α 3 = 10 -2.4 we can define a variable ξ −1 R/R o and represent it as a function of Oh, where ξ is a function derived from (2.26) which tends to 1 for very small Oh and Bo, given by That plot is given in figure 8, where we represent the variable ξ −1 R/R o (calculated with the fixed values (Oh * , α 3 , β 3 ) of the best fitting) as a function of Oh. This first representation shows a strong variability of data in the vicinity of Oh * . However, this variability does not seem entirely stochastic: one may observe that despite that the data have been made non-dimensional with fixed values of (Oh * , α 3 , β 3 ), the closeness of some data sets to curves corresponding to different values of these parameters suggests that those data sets could be independently fitted to our model, as the dashed curves reveal. Intriguingly, the data sets showing largest variability are the ones corresponding to numerical simulations.

Experimental data sets: universality in terms of Oh and Bo
The above observations lead to the imperative need to analyse the origin of the deviations in the numerical data sets. In fact, a careful examination of the nature of the different data sources reveal an apparently paradoxical finding: if one separately fits all the experimental data, i.e. excluding all data from numerical simulations, one obtains a robust best fit around our model for Oh * = 0.0296, α 3 = 0.00405 and β 3 = −0.000158. This yields a minimum normalized standard deviation of 9 % (significantly smaller than the 24 % obtained including numerical simulations) and a prefactor k r = 0.186 Oh * 2 . Figure 9 shows the resulting optimum fitting, and the distribution of data around the model compared with a normal distribution. Observe that the data of Brasz et al. (2018), and of  and Séon & Liger-Belair (2017), provide very valuable resources for comparison around Oh * despite the experimental errors and the particular difficulty of that region due to the droplet smallness.  Oh Spiel 1995Blanchard 1989Garner 1954Hayami 1958Tedesco 1979Sakai 1989 Seon 2017 Brasz ( In contrast, when the whole set of available data is compared with our model, the dispersion around Oh * explodes as anticipated in figure 8: this is illustrated in figure 9(b). The statistics of both fittings are also given, noting that the distributions are nearly invariant along the Oh range except at the vicinity of Oh * . This finding highlights once more the strong sensitivity of the focusing effect on the geometry and the initial conditions: (i) the artificial shape of the initial rim after bursting, which critically determines the initial surface energy content as the pioneering work of Duchemin et al. (2002) revealed; (ii) initial bubble shape; and (iii) the absence or presence of gravity along the evolution. In addition, the numerical accuracy as well as the models for free surface and surface tension reconstruction, etc., also affect the results. Note that the precision of our simulations does not avoid their maximum deviation: the fact that we used an initial edge shape thickness (initial minimum film thickness) equal to 0.1R o and initial hole radius 0.15R o , leads to the largest initial meridional radius of curvature among the published initial geometries (see figure 7), which shifts Oh * from 0.03 to 0.034.
In contrast, the actual initial surface energy content after film bursting in real experiments is expected to be fixed by a self-similar process (Eggers & Fontelos 2015, and references therein) with a universal result only in terms of Oh and Bo alone. In effect, despite the relative paucity of data around Oh * and the fact that all the experiments available in this work have been performed with air at atmospheric pressure, naturally including interaction with the gaseous environment, gas compressibility, surface effects (presence of contaminants, surfactants, particles, etc.) and nonlinear breakdown and local detachment of the first droplet, the scatter compared with that of the numerical data is minimal. This supports the general assumption that this problem is fundamentally determined in nature by two non-dimensional parameters only: Oh and Bo.
On the other hand, the considerably less numerous data available for the ejection velocity V could be similarly represented in the form ξ 1/2 V/V o , which tends to V/V o for very small Oh and Bo. However, since our model predicts the ejection velocity at the droplet formation point, while the available data correspond to the ejection velocity measured at a fixed point in space (usually, the initial level of the free surface), a direct comparison with our model predictions is not possible. To compensate for the effect of the different measurement point, which necessarily implies the additional action of viscous dissipation and the work of gravity, we propose a simple velocity correction of the form V/V o = k v Ψ −1/2 (1 + k 1 Bo + k 2 Oh) −1 , (2.28) which, under the same approach as the droplet size approach, provides an optimal fit to the experimental data of  for k 1 = 2.27 and k 2 = 16, with a minimum normalized standard deviation below 10 % and a prefactor k v = 0.39 (see figure 10a). Including the entire data set available for ejection velocity (numerical data from Berny et al. (2020) and our own), we observe the same drastic increase of the normalized standard deviation for droplet size. The statistical distribution of experimental data around our model, with an approximately normal distribution and a standard deviation below 10 %, provides solid statistical grounds (i.e. the null hypothesis is virtually excluded) for its consistency and  Figure 10. The velocity of the first emitted droplet, expressed as Ψ 1/2 (1 + k 1 Bo + k 2 Oh)V/V o , with k 1 = 2.27, with k 2 = 16, for (a) the experimental data of  and Duchemin et al. (2002), after optimal data collapse according to (2.28) and (b) the whole set of available experimental and numerical data, using the same fitting as in panel (a). correctness: our model agrees with the available experimental data for the whole (Oh, Bo) range. We emphasize that (3.1) is not a fitting polynomial expression: it reflexes the physics of the phenomenon as described in this work using just two non-dimensional parameters. The main advantage of (3.1) compared with other proposed models which provide good results but are not valid beyond Oh * (Gañán-Calvo 2017, 2018; Gordillo & Rodriguez-Rodriguez 2019) is its uniform validity and continuity along the whole (Oh, Bo) range. Among an ample variety of applications, this provides a fundamental tool to describe the statistical distribution of the sea spray as a function of the observed bubble size distribution near the sea surface. A similarly accurate expression for the ejection velocity is not currently possible from available data, due to its strong dependency on the point where it is measured. Based on data at the level of the unperturbed free surface, we propose the following expression from an optimum fitting: V = 0.39V o Oh − Oh * 2 + 0.004Oh + 0.00017Bo −1/2
(3.2) Of course, the complexity of the phenomenon can be augmented including surface viscosity (Ponce-Torres et al. 2017), Marangoni and non-Newtonian effects (Sanjay et al. 2021), the presence of immiscible liquids (Yang, Tian & Thoroddsen 2020) or the statistics of the emitted droplets following the first one. In this regard, Berny et al. (2020) presented a detailed analysis of all issued droplets from the scaling analysis of the first. Despite these are topics of subsequent studies, our results entirely capture the physics needed to describe this phenomenon.