2.1. Governing equations

We investigate the collapse of an open bubble cavity at the interface in a viscoelastic medium (of figure 1) using an axisymmetric domain with incompressible fluids. Length scales are normalised using the initial bubble radius giving $\mathcal {L} = \tilde {\mathcal {L}}R_0$ as characteristic length, and the time is normalised using the inertiocapillary time scale $\tau _\gamma = \sqrt {\rho _s {R_0}^3/\gamma }$ giving $t = \tilde {t}\tau _{\gamma }$ . These normalisations yield an inertiocapillary velocity scale $u_{\gamma } = \sqrt {\gamma / \rho _{s} R_0}$ for the velocity field $\boldsymbol {u} = \tilde {\boldsymbol {u}}u_\gamma$ . Lastly, all stresses are normalised using the Laplace pressure scale, $\boldsymbol {\sigma } = \tilde {\boldsymbol {\sigma }}\sigma _\gamma$ , where $\sigma _\gamma = \gamma /R_0$ . Here, as usual, non-dimensionalised quantities are denoted with a tilde, although from here onwards, we drop the tilde, and all equations are thus dimensionless in the current section. Throughout the manuscript, we use the subscripts $s$ , $p$ and $g$ to denote liquid solvent, polymer and gas, respectively. The governing mass and momentum conservation equations for the liquid phase read as

(2.1) \begin{equation} \boldsymbol {\nabla \cdot u} =0,\,\text {and} \end{equation}

(2.2) \begin{equation} \frac {\partial \boldsymbol {u}}{\partial t} + \boldsymbol {\nabla \cdot } \left (\boldsymbol {u}\boldsymbol {u}\right ) = -\boldsymbol {\nabla }p + \boldsymbol {\nabla \cdot }\left (\boldsymbol {\sigma _s}+\boldsymbol {\sigma _p}\right ), \end{equation}

where the Newtonian contribution (coming from the solvent) $\boldsymbol {\sigma _s}$ is

(2.3) \begin{align} \boldsymbol {\sigma _{s}} = 2 Oh_s \boldsymbol {\mathcal {D}}, \end{align}

with $\boldsymbol {\mathcal {D}} = (\boldsymbol {\nabla u} + ( \boldsymbol { \nabla u} ) ^T )/2$ representing the symmetric part of the velocity gradient tensor – equal to half of the rate-of-strain tensor. The non-Newtonian contribution $\boldsymbol {\sigma _{p}}$ arises from the presence of polymers in the fluid. We emphasise that, although we refer to $\boldsymbol {\sigma _{p}}$ as ‘polymeric stresses’ in the context of dilute polymer liquids, this concept extends to any deformable microstructure within the fluid that responds to flow (Saramito Reference Saramito2007; Snoeijer et al. Reference Snoeijer, Pandey, Herrada and Eggers2020; Balasubramanian et al. Reference Balasubramanian, Sanjay, Jalaal, Vinuesa and Tammisola2024; França et al. Reference França, Jalaal and Oishi2024). To characterise the deformation of these microstructures, we introduce the conformation tensor $\boldsymbol {\mathcal {A}}$ , an order parameter that evolves from an initial identity state $\boldsymbol {\mathcal {A}} = \boldsymbol {\mathcal {I}}$ (figure 1 a-ii). Here, we employ the Oldroyd-B model, which represents the simplest conformation tensor-based constitutive equation for viscoelastic fluids (Oldroyd Reference Oldroyd1950; Bird, Armstrong & Hassager Reference Bird, Armstrong and Hassager1977; Snoeijer et al. Reference Snoeijer, Pandey, Herrada and Eggers2020; Stone, Shelley & Boyko Reference Stone, Shelley and Boyko2023; Boyko & Stone Reference Boyko and Stone2024). This model assumes a linear relationship between elastic stresses and polymeric deformation

(2.4) \begin{align} \boldsymbol {\sigma _{p}} = Ec \left (\boldsymbol {\mathcal {A}} - \boldsymbol {\mathcal {I}}\right ), \end{align}

where $Ec$ is the elastocapillary number (1.3), representing the strength of the polymers analogous to a dimensionless elastic modulus. Note that, even though the polymeric stresses $\boldsymbol {\sigma _{p}}$ grow linearly with $\boldsymbol {\mathcal {A}}$ , the polymeric deformations $\boldsymbol {\mathcal {A}}$ can be highly nonlinear. Naturally, in the limit of $Ec = 0$ , the polymeric stress would vanish, and the system will give a viscous Newtonian dictated by the solvent Ohnesorge number $Oh_s$ (see (2.3)).

Additionally, the conformation tensor $ \boldsymbol {\mathcal {A}}$ relaxes to its base state $\boldsymbol {\mathcal {I}}$ over time due to thermal effects. Once more, using the Oldroyd-B model, $ \boldsymbol {\mathcal {A}}$ follows a linear relaxation law (i.e. the rate of change of $\boldsymbol {\mathcal {A}}$ in the Lagrangian frame is linear in $\boldsymbol {\mathcal {A}}$ )

(2.5) \begin{align} \stackrel {\kern5pt{\boldsymbol {\nabla }}}{\boldsymbol {\mathcal {A}}} \,= - \frac {1}{De} \left ( \boldsymbol {\mathcal {A}} - \boldsymbol {\mathcal {I}} \right ), \end{align}

where

(2.6) \begin{align} \stackrel {\kern5pt{\boldsymbol {\nabla }}}{\boldsymbol {\mathcal {A}}} \, \equiv \frac {\partial \boldsymbol {\mathcal {A}}}{\partial t} + \left (\boldsymbol {u\cdot \nabla }\right )\boldsymbol {\mathcal {A}} - \boldsymbol {\mathcal {A}\cdot }\left (\boldsymbol {\nabla u}\right ) - \left (\boldsymbol {\nabla u}\right )^T\boldsymbol {\cdot \mathcal {A}} ,\end{align}

is the frame-invariant upper convected Oldroyd derivative of second-rank tensor $\boldsymbol {\mathcal {A}}$ , and $De = \lambda /\tau _\gamma$ (defined in (1.4)) is the Deborah number, representing the ratio of the polymer relaxation time $\lambda$ to the process time scale $\tau _\gamma$ . We note that, while the Oldroyd-B model is nonlinear in terms of the velocity field and its gradient, both the stress term and its relaxation law remain linear in $\boldsymbol {\mathcal {A}}$ . This characteristic contrasts with models such as the Giesekus model, which involves a quadratic term $\boldsymbol {\mathcal {A}\cdot \mathcal {A}}$ (Giesekus Reference Giesekus1982), or the finite extensible nonlinear elastic (FENE) models, which include a nonlinear term involving a finite-extensibility parameter $L$ (Bird, Dotson & Johnson Reference Bird, Dotson and Johnson1980). Therefore, the Oldroyd-B model is often referred to as ‘quasi-linear’ (Davoodi et al. Reference Davoodi, Lerouge, Norouzi and Poole2018; Alves et al. Reference Alves, Oliveira and Pinho2021).

The Deborah number characterises the polymeric liquid’s memory. It is instructive to note that, in the limit of $De \to \infty$ , polymeric liquids have permanent memory and the dissolved polymers undergo affine motion (see (2.5) and Snoeijer et al. Reference Snoeijer, Pandey, Herrada and Eggers2020; Stone et al. Reference Stone, Shelley and Boyko2023; Boyko & Stone Reference Boyko and Stone2024)

(2.7) \begin{align} \stackrel {\kern5pt{\boldsymbol {\nabla }}}{\boldsymbol {\mathcal {A}}} \,= 0, \end{align}

indicating that they follow the flow and deform according to the velocity field. In this limit, for finite $Ec$ , the Oldroyd-B model is equivalent to the damped neo-Hookean model (also known as the Kelvin–Voigt model) for solids (Snoeijer et al. Reference Snoeijer, Pandey, Herrada and Eggers2020). Conversely, at $De = 0$ , polymeric liquids have no memory of their initial condition and subsequent deformations, relaxing immediately to the base state. For non-infinite $Ec$ values, polymeric stresses vanish, resulting in a Newtonian response (2.4) governed by the solvent Ohnesorge number $Oh_s$ (see (2.3)). It is, therefore, surprising that both $Ec = 0$ and $De = 0$ (figure 1 b) represent Newtonian responses, irrespectively of the corresponding other parameter.

Equations (2.4) and (2.5) can be combined to get

(2.8) \begin{align} De \stackrel {\kern-6pt{\boldsymbol {\nabla }}}{\boldsymbol {\sigma _p}} +\,\,\boldsymbol {\sigma _p} = 2Oh_p {\boldsymbol {\mathcal {D}}}, \end{align}

where $Oh_p = Ec \times De$ is the polymeric Ohnesorge number (1.5). Consequently, in the limit $De \to 0$ at fixed $Oh_p$ (e.g. moving along constant $Oh_p$ lines in figure 1 b), the system exhibits a viscous Newtonian response with a total dimensionless viscosity of $Oh_s+Oh_p$ .

The Oldroyd-B model, despite its widespread use due to its simplicity, fails to capture several important physical phenomena (Snoeijer et al. Reference Snoeijer, Pandey, Herrada and Eggers2020). It is inadequate to describe shear-thinning behaviour in polymeric liquids (Yamani & McKinley Reference Yamani and McKinley2023) and erroneously predicts unbounded stress growth in strong extensional flows (McKinley & Sridhar Reference McKinley and Sridhar2002; Eggers, Herrada & Snoeijer Reference Eggers, Herrada and Snoeijer2020). The numerical discretisation of Oldroyd-B (§ 2.2) also features an implicit stress regularisation due to the finite grid size (Renardy & Thomases Reference Renardy and Thomases2021) – similar in spirit to the implicit slip regularisation of the contact line singularity (Afkhami et al. Reference Afkhami, Buongiorno, Guion, Popinet, Saade, Scardovelli and Zaleski2018; Fullana et al. Reference Fullana, Kulkarni, Fricke, Popinet, Afkhami, Bothe and Zaleski2024). These limitations can be addressed by incorporating finite polymer extension, for example, by increasing the effective $Ec$ as the polymer approaches full extension (Hinch & Harlen Reference Hinch and Harlen2021; Zinelis et al. Reference Zinelis, Abadie, McKinley and Matar2024). Various extensions of the Oldroyd-B equations have been developed to account for such nonlinearity, either in (2.4) and (2.5) or in the solvent contribution in (2.3) (de Gennes Reference de Gennes1974; Tanner Reference Tanner2000; McKinley & Sridhar Reference McKinley and Sridhar2002; Alves et al. Reference Alves, Oliveira and Pinho2021). In this study, we employ the Oldroyd-B model to include the two primary effects of the polymer addition: the additional stress ( $Ec$ ) and polymeric liquid memory ( $De$ ) (Snoeijer et al. Reference Snoeijer, Pandey, Herrada and Eggers2020). Our aim is to provide a comprehensive understanding of the entire $Ec$ – $De$ parameter space (figure 1 b). However, it is crucial to note that the Oldroyd-B model, while serving as a useful baseline, cannot accurately reproduce the finite-time breakup of viscoelastic filaments (Eggers et al. Reference Eggers, Herrada and Snoeijer2020) or the full complexity of interface rupture (Lohse & Villermaux Reference Lohse and Villermaux2020). These limitations warrant caution when interpreting the final stages of jet thinning and droplet formation, particularly in scenarios involving strong polymer stretching.

2.2. Methods

We employ the open-source software Basilisk C (Popinet & collaborators Reference Popinet2013–Reference Popinet2024; Popinet Reference Popinet2015) to solve the governing equations outlined in § 2.1. To solve the Oldroyd-B viscoelastic constitutive relation (2.8), Basilisk C uses the log-conformation method (Fattal & Kupferman Reference Fattal and Kupferman2004) implemented by López-Herrera et al. (Reference López-Herrera, Popinet and Castrejón-Pita2019) which has been used extensively at finite $De$ (Turkoz et al. Reference Turkoz, Lopez-Herrera, Eggers, Arnold and Deike2018, Reference Turkoz, Stone, Arnold and Deike2021). To explore the entire $Ec$ – $De$ parameter space (figure 1 c), we have extended the log-conformation formulation to solve (2.4) and (2.5). In the spirit of Basilisk C, this code is detailed open source in Sanjay (Reference Sanjay2024). The rest of the governing equations are solved using the one-fluid approximation (Tryggvason, Scardovelli & Zaleski Reference Tryggvason, Scardovelli and Zaleski2011), with surface tension incorporated as singular body force at the liquid–gas interface (Brackbill, Kothe & Zemach Reference Brackbill, Kothe and Zemach1992). To account for the gas phase, in addition to the dimensionless parameters described in §§ 1 and 2.1, we maintain constant density and viscosity ratios of $\rho _{r} = \rho _{g}/\rho _{s} = 10^{-3}$ and $\eta _{r} = \eta _{g}/\eta _{s} = 2 \times 10^{-2}$ , respectively. The liquid–gas interface is tracked using the volume of fluid (VoF) method, governed by the advection equation

(2.9) \begin{align} \frac {\partial \Psi }{\partial t} + \boldsymbol {\nabla \cdot }\left (\Psi \boldsymbol {u}\right ) = 0, \end{align}

where $\Psi$ represents the VoF colour function. We implement a geometric VoF approach, reconstructing the interface at each time step and applying surface tension forces as singular forces (Brackbill et al. Reference Brackbill, Kothe and Zemach1992; Popinet Reference Popinet2009)

(2.10) \begin{align} \boldsymbol {f}{\gamma } \approx \kappa \boldsymbol {\nabla } \Psi , \end{align}

with curvature $\kappa$ calculated using the height-function method (Popinet Reference Popinet2018). The explicit treatment of surface tension imposes a time step constraint based on the smallest capillary wave oscillation period (Popinet Reference Popinet2009). Yet another time step restriction, usually more relaxed than the surface tension one, comes from the explicit treatment of the polymeric stress term $\boldsymbol {\sigma }_p$ . We impose no-penetration and free-slip conditions at wall boundaries to avoid wall-shear effects, with outflow conditions at the top boundary to prevent droplet rebound. Pressure gradients are set to zero at domain boundaries for both liquid and gas phases.

The initial bubble shape is determined by solving the Young–Laplace equations for quasi-static equilibrium (Toba Reference Toba1959; Princen Reference Princen1963; Sanjay Reference Sanjay2022; Villermaux et al. Reference Villermaux, Wang and Deike2022). While the shape’s asymmetry increases with the Bond number $Bo$ , we focus on the limit $Bo \to 0$ , setting $Bo = 0.001$ to regularise the singularity at the sphere–plane intersection. This results in a near-spherical initial cavity shape (figure 1 a-i). We stress that, here, we assume that the bubble has resided at the liquid–gas interface for a duration far exceeding the polymeric medium’s relaxation time, ensuring that elasticity does not influence the initial configuration (Balasubramanian et al. Reference Balasubramanian, Sanjay, Jalaal, Vinuesa and Tammisola2024). During the bubble cap bursting, the film cap retracts almost instantaneously (once again, we neglect the influence of elasticity), after which the capillary waves are generated. As we are interested only in the bubble cavity collapse, the simulations begin with an open cavity without the thin cap (figure 1 a-ii), as also done similarly in recent studies (Deike et al. Reference Deike, Ghabache, Liger-Belair, Das, Zaleski, Popinet and Séon2018; Gordillo & Rodríguez-Rodríguez Reference Gordillo and Rodríguez-Rodríguez2019; Sanjay et al. Reference Sanjay, Lohse and Jalaal2021). The computational domain spans $8R_0 \times 8R_0$ , discretised using quadtree grids with adaptive mesh refinement (Popinet Reference Popinet2009). Error tolerances for the VoF colour function, curvature, velocity and order parameter $\boldsymbol {\mathcal {A}}$ are set to $10^{-3}$ , $10^{-6}$ , $10^{-3}$ and $10^{-3}$ , respectively.

In this work, following our earlier study (Sanjay et al. Reference Sanjay, Lohse and Jalaal2021), most simulations maintain a minimum grid size of $\unicode{x1D6E5} = R_0/512$ , which dictates that, to get consistent results, 512 cells are required across the bubble radius while using uniform grids. We have also used an increased resolution ( $\unicode{x1D6E5} = R_0/1024$ for high $De$ cases and $\unicode{x1D6E5} = R_0/2048$ near transitions) as needed. These resolutions are consistent with previous studies by Berny et al. (Reference Berny, Deike, Séon and Popinet2020, Reference Berny, Popinet, Séon and Deike2021) on bubble bursting and Turkoz et al. (Reference Turkoz, Lopez-Herrera, Eggers, Arnold and Deike2018, Reference Turkoz, Stone, Arnold and Deike2021) on visco-elastic thinning with a maximum level of resolution of 14 (for $\unicode{x1D6E5} = R_0/2048$ and domain size $L_0 = 8R_0$ ). We have carried out extensive grid independence studies to ensure that changing the grid size does not influence the results (see Appendix C). We refer the readers to Popinet (Reference Popinet2015), Sanjay (Reference Sanjay2022) and Sanjay & Dixit (Reference Sanjay and Dixit2024) for further details of the numerical method used in this work.

3.1. Polymeric liquids with permanent memory

We begin our analysis by considering the limit of $De \to \infty$ , where the polymeric solutions feature affine motion (2.7) and maintain a permanent memory of their initial condition and subsequent deformations without relaxation during the process time scale. Figure 2 illustrates representative cases in viscoelastic media for $Oh_s = 0.025$ and varying elastocapillary numbers ( $Ec$ ). The figure presents a temporal evolution of the interface profile (green line) alongside the velocity magnitude on the left and the trace of elastic stress $\boldsymbol {\sigma }_p$ on the right. Remarkably, despite all cases exhibiting a total Ohnesorge number of infinity ( $Oh_s + Oh_p \to \infty$ ), which typically implies highly viscous behaviour (see figure 13), low $Ec$ scenarios demonstrate a dynamics qualitatively resembling Newtonian fluids. In these cases, capillary waves drive the collapse of a bubble cavity, converging at its bottom to form a Worthington jet that subsequently fragments into droplets (see figure 2 a). Intuitively, the elastic stresses are concentrated near the axis of symmetry where the strain is maximum (Turkoz et al. Reference Turkoz, Lopez-Herrera, Eggers, Arnold and Deike2018; Eggers et al. Reference Eggers, Herrada and Snoeijer2020). The process concludes within a finite time scale ( $\sim \tau _\gamma$ ), resulting in a regular limit as $Ec \to 0$ . As a result, the system’s behaviour deviates gradually from the Newtonian case at $Ec = 0$ , exhibiting a continuous transition as the elasticity increases. This absence of singularity contrasts with elastic Taylor–Culick-type retractions, where an infinite process time scale allows the elastic stresses to develop, leading to distinct behaviours for $Ec = 0$ and $Ec \to 0$ (Bertin et al. Reference Bertin, Sanjay, Oratis and Snoeijer2024), i.e. a singular limit.

Figure 2.

Temporal evolution of the bubble cavity collapse at $De \to \infty$ and $Oh_s = 0.025$ for $Ec =$ (a) $0.0001$ , (b) $0.01$ and (c) $0.1$ . The colour scheme in the left panel of each snapshot represents the magnitude of the velocity field normalised by the inertiocapillary velocity, while the right panel of each snapshot shows the trace of the elastic stress $\boldsymbol {\sigma }_p$ that represents twice the elastic energy stored in polymeric deformations on a $\log _{10}$ scale. See also the supplementary movies SM1.

We stress that in this limit, the jet breakup occurs due to finite grid resolution in our numerical code (Lohse & Villermaux Reference Lohse and Villermaux2020; Chirco et al. Reference Chirco, Maarek, Popinet and Zaleski2022; Kant et al. Reference Kant, Pairetti, Saade, Popinet, Zaleski and Lohse2023). We cannot differentiate between a case of drop detachment from the jet or the case when they are still connected through a thin filament – also known as the beads-on-a-string structure (Clasen et al. Reference Clasen, Eggers, Fontelos, Li and McKinley2006; Pandey et al. Reference Pandey, Kansal, Herrada, Eggers and Snoeijer2021; Hosokawa et al. Reference Hosokawa, Kamamoto, Watanabe, Kusuno, Kobayashi and Tagawa2023; Zinelis et al. Reference Zinelis, Abadie, McKinley and Matar2024). Although current simulations fully resolve other aspects, they cannot resolve these finest threads, which may have subgrid cell sizes depending on $Ec$ . At higher grid resolutions, we expect to recover the beads-on-a-string configuration, as the Oldroyd-B model does not yield a finite time breakup singularity in the infinite $De$ regime, instead converging to a finite filament (Turkoz et al. Reference Turkoz, Lopez-Herrera, Eggers, Arnold and Deike2018, Reference Turkoz, Stone, Arnold and Deike2021; Eggers et al. Reference Eggers, Herrada and Snoeijer2020). To prevent infinite thread thinning, a nonlinear elastic model could also be employed (see § 2.1 for further discussions).

As $Ec$ increases, we observe jet formation without droplet ejection (figure 2 b). At higher $Ec$ values, even jet formation is suppressed due to elevated elastic resistance (figure 2 c). Notably, while polymeric effects significantly influence the dynamics after the convergence of capillary waves (figure 2, $t/\tau _\gamma = 0.8, 1.2$ ), the propagation of capillary waves (figure 2, $t/\tau _\gamma = 0.1, 0.4$ ) remains largely unaffected. Figure 3(a) quantifies the trajectories of these capillary waves across three orders of magnitude variation in $Ec$ at two different $Oh_s$ . The capillary wave speed is independent of both liquid and polymeric control parameters, mirroring the behaviour observed in Newtonian media (Gordillo & Rodríguez-Rodríguez Reference Gordillo and Rodríguez-Rodríguez2019) and contrasting those for viscoplastic media (Sanjay et al. Reference Sanjay, Sen, Kant and Lohse2022). The independence of capillary wave speed on the polymeric control parameters has also been reported in experiments (Cabalgante-Corrales et al. Reference Cabalgante-Corrales, Muñoz-Sánchez, López-Herrera, Cabezas, Vega and Montanero2025). Following capillary wave collapse, the Worthington jet initially elongates to a maximum length ( $L_{{max}}$ ) before retracting. As shown in figure 3(b) for $Oh_s = 0.04$ , $L_{{max}}$ decreases with increasing $Ec$ due to stronger resistive stresses.

Figure 3.

(a) Trajectory of the maximum curvature capillary wave parameterised using the angle $\theta _c(t)$ as depicted in the inset at $De \to \infty$ for different $Oh_s$ and $Ec$ . (b) Evolution of the jet length $L(t)$ at $Oh_s = 0.04$ and $De \to \infty$ for different $Ec$ .

Figure 4(a) presents a phase map of $L_{{max}}$ , compiled from approximately 100 simulations. For Newtonian liquids, $L_{{max}}$ peaks near $Oh_s \approx 0.03$ , corresponding to the value of observed hydrodynamic singularities (Zeff et al. Reference Zeff, Kleber, Fineberg and Lathrop2000; Lohse Reference Lohse2003; Eggers & Fontelos Reference Eggers and Fontelos2015; Yang, Tian & Thoroddsen Reference Yang, Tian and Thoroddsen2020), before decreasing at higher $Oh_s$ (Duchemin et al. Reference Duchemin, Popinet, Josserand and Zaleski2002; Deike et al. Reference Deike, Ghabache, Liger-Belair, Das, Zaleski, Popinet and Séon2018; Gordillo & Rodríguez-Rodríguez Reference Gordillo and Rodríguez-Rodríguez2019). Jet formation ceases altogether beyond a critical value of $Oh_c = 0.11$ (Sanjay et al. Reference Sanjay, Lohse and Jalaal2021) (defined here when $L_{{max}} \lt 0.3R_0$ ). As $Ec$ increases, viscoelastic effects become significant. The value of $L_{{max}}$ decreases monotonically with $Ec$ due to increased elastic resistance, with jet formation suppressed beyond $Ec=0.086$ . Unlike the non-monotonic relationship between $L_{{max}}$ and $Oh_s$ , where increasing $Oh_s$ initially produces thinner and faster jets, the $L_{{max}}(Ec)$ relationship remains consistently monotonic. Even the $Oh_s$ -sensitive singular Worthington jets disappear with increasing $Ec$ . Notably, the critical $Ec$ values for these transitions appear to be largely independent of $Oh_s$ , in contrast to the $Oh_s$ -dependent behaviour observed in the Newtonian limit.

Figure 4.

(a) The maximum jet length $L_{{max}}$ at $De \to \infty$ in the $Ec$ - $Oh_s$ phase space, depicted by the colour map, where the lighter region corresponds to higher values. For the Newtonian liquid $(Ec \to 0)$ , the jetting transition occurs at $Oh_s = 0.11$ , denoted by the horizontal dotted line. Due to the elastic effects, this transition occurs at $Ec = 0.086$ , as depicted by the vertical dotted line. (b) The size of the first droplet at $De \to \infty$ in the $Ec$ - $Oh_s$ phase space. For the Newtonian liquid, the dropping transition is observed at $Oh_s = 0.0375$ , denoted by the horizontal dotted line. Further, the transition due to elastic effects is very sensitive to $Oh_s$ and is shown by the inclined dotted line.

The emerging Worthington jet may eject multiple droplets. For Newtonian liquids, predictions, for the first droplets’ size $r_d$ are well understood (see Appendix A and Gañán-Calvo Reference Gañán-Calvo2017; Blanco-Rodríguez & Gordillo Reference Blanco-Rodríguez and Gordillo2020). Here, $r_d$ decreases with $Oh_s$ until $Oh_s \approx 0.0375$ , beyond which the droplet breaks from the jet due to the Rayleigh–Plateau instability and falls downwards. Our analysis focuses on droplets propagating away from the source, excluding those with downward velocity upon breakup (observed in Newtonian media for $0.0375 \lt Oh_s \lt 0.045$ ). For elastic cases, despite unresolved filaments connecting droplets and jets, we have rigorously verified the convergence of the first droplet’s size to at least 10 % accuracy. Figure 4(b) illustrates a phase map of the first droplet’s size $r_d$ , revealing intriguing differences from the jet behaviour. While $r_d$ follows the same trend with $Oh_s$ observed at Newtonian limits and remains invariant of $Ec$ below critical values, the critical $Ec$ for droplet suppression differs from that of jet suppression. As the jet width is determined solely by $Oh_s$ , independently of $Ec$ , the first emerging droplet’s size also remains independent initially. However, as $Ec$ increases further, rising elastic stresses suppress droplet formation more abruptly than jet formation. The critical values $Ec_d$ for the transition between jet formation with and without droplet breakup (dropping transition) are sensitive to $Oh_s$ , with the critical $Ec_d$ decreasing as $Oh_s$ increases. This trend is in stark contrast with the transition from jet formation to jet suppression (jetting transition), where critical $Ec$ values remain largely $Oh_s$ -independent.

3.2. Polymeric liquids with poor memory

This section examines the dynamics in media with a poor memory of the initial conditions and subsequent deformations ( $De \to 0$ ). For sufficiently small Deborah numbers $De$ , the polymers relax rapidly, resulting in elastic stresses of the polymeric liquid that are considerably lower than those observed in cases where $De \to \infty$ . The stress relaxation also results in the dissipation of elastic energy stored in stretched polymers. Figure 5 illustrates representative cases for $De = 0.01$ , showcasing three distinct regimes as a function of the elastocapillary number ( $Ec$ ). The figure presents a temporal evolution of the interface profile (green line) alongside velocity magnitude on the left and the trace of elastic stress $\boldsymbol {\sigma }_p$ on the right for $Ec = 1$ , $5$ and $10$ . For $Ec = 1$ (figure 5 a), we observe a slender Worthington jet that forms a droplet. As $Ec$ increases to $5$ (figure 5 b), the jet persists but fails to produce a droplet. At $Ec = 10$ (figure 5 c), jet formation is completely suppressed, with the interface showing only slight deformations during cavity relaxation. The qualitative trends with respect to the elastocapillary number ( $Ec$ ) remain consistent as compared with those in § 3.1. However, the critical $Ec$ values for different regimes differ markedly from those observed at $De \to \infty$ . Notably, jet formation and droplet production persist at $Ec = 1$ (figure 5 a), despite this value being an order of magnitude higher than the critical $Ec$ for the jetting transition at infinite $De$ . This difference underscores the dependence of transition thresholds on $De$ .

Figure 5.

Temporal evolution of bubble cavity collapse at $De = 0.01$ and $Oh_s = 0.025$ for $Ec =$ (a) $1$ , (b) $5$ and (c) $10$ . The colour scheme in the left panel of each snapshot represents the magnitude of the velocity field normalised by the inertiocapillary velocity, while the right panel of each snapshot shows the trace of the elastic stress $\boldsymbol {\sigma }_p$ that represents twice the elastic energy stored in polymeric deformations on a $\log _{10}$ scale. See also the supplementary movies SM2.

To further interpret the jetting dynamics and drop formation, figure 6 presents phase maps illustrating the behaviour of maximum jet lengths ( $L_{{max}}$ ) and first droplet sizes ( $r_d$ ) for $De = 0.01$ . Figure 6(a) shows $L_{{max}}$ across a range of $Ec$ and $Oh_s$ values. For low $Ec$ , $L_{{max}}$ shows Newtonian-like $Oh_s$ dependence. As $Ec$ increases, $L_{{max}}$ decreases monotonically until jet formation ceases beyond an $Oh_s$ -independent critical $Ec_j$ , mirroring the infinite $De$ limit behaviour. Figure 6(b) maps $r_d$ , showing $Ec$ -independent droplet sizes that are equal to values at the Newtonian limit until near the transition point, where droplet formation is suppressed. For $De \ll 1$ , the critical $Ec_d$ for the dropping transition exhibits minimal $Oh_s$ -dependence, contrasting with the $Oh_s$ -sensitive behaviour at infinite $De$ . Comparing these results with the $De \to \infty$ limit reveals persistent fundamental regimes across different $De$ values, but the transition thresholds are highly sensitive to the polymeric liquid’s relaxation time. Critical $Ec$ values for both jet and droplet suppression are significantly higher at low $De$ compared with the infinite $De$ limit, indicating that rapid relaxation of polymeric stresses allows jet and droplet formation at higher $Ec$ values. This low $De$ behaviour suggests an interplay between elastic and viscous effects, explored further in § 4.

Figure 6.

(a) The maximum jet length $L_{{max}}$ at $De = 0.01$ in the $Ec$ - $Oh_s$ phase space, depicted by the colour map, where the lighter region corresponds to higher values. For the Newtonian liquid, the jetting transition occurs at $Oh_s = 0.11$ , denoted by the horizontal dotted line. Due to the elastic effects, this transition occurs at $Ec = 9.3$ , as depicted by the vertical dotted line. (b) The size of the first droplet at $De = 0.01$ in the $Ec$ - $Oh_s$ phase space. For the Newtonian liquid, the dropping transition is observed at $Oh_s = 0.0375$ , denoted by the horizontal dotted line. Further, the $Oh_s$ -independent transition due to elastic effects occurs at $Ec= 2.5$ , as shown by the vertical dotted line.

4.1. Summary of the different regimes

The transitions between these regimes depend on both $Ec$ and $De$ , exhibiting markedly different characteristics in two limiting cases: $De \to \infty$ and $De \to 0$ . Figure 7 maps these transitions in the elastocapillary–Deborah number ( $Ec$ – $De$ ) phase space, delineating the boundaries between droplet-forming jets and jets without droplets. Figure 8 complements this by illustrating the transition to complete jet suppression. Notably, the infinite $De$ asymptotic behaviour extends down to $De \approx 1$ , reflecting that polymers lack sufficient time to relax when relaxation times exceed the process time scale.

Figure 7.

The elastocapillary–Deborah number ( $Ec$ – $De$ ) phase map delineating the transition between the regimes: (i) jets forming droplets and (ii) jets without droplet formation. The data points represent the critical elastocapillary number $Ec_d(De, Oh_s)$ at which this transition occurs. The transition behaviour exhibits distinct characteristics in different limits: as $De \to \infty$ , the transition occurs at a constant $Ec$ which is highly sensitive to $Oh_s$ (see the grey dashed line showing $Ec_d \sim De^0$ ), while for $De \to 0$ , the transition is $Oh_s$ -independent and occurs at constant $Oh_p$ (see the grey solid line showing $Ec_d \sim De^{-1}$ , i.e. $Oh_{p,d} \sim De^0$ ).

Figure 8.

(a) The elastocapillary–Deborah number ( $Ec$ – $De$ ) and (b) the polymeric Ohnesorge–Deborah number ( $Oh_p$ - $De$ ) phase map delineating the transition between the regimes: (ii) jets without droplet formation and (iii) absence of jet formation. The data points represent the $Oh_s$ -independent critical elastocapillary number $Ec_j(De)$ at which this transition occurs. The transition behaviour exhibits distinct characteristics in different limits: as $De \to \infty$ , the transition occurs at a constant $Ec$ (see grey dashed line showing $Ec_d \sim De^0$ ), while for $De \to 0$ , the transition occurs at constant $Oh_p$ (see grey solid line showing $Ec_d \sim De^{-1}$ , i.e. $Oh_{p,d} \sim De^0$ ).

For polymeric liquids with long relaxation times ( $De \gg 1$ ), we observe that:

1. (i) the dropping transition occurs at $Ec_d(Oh_s)$ , with strong $Oh_s$ dependence (figures 4b and 7); and

2. (ii) the jetting transition occurs at $Ec_j \approx 0.086$ , independent of $Oh_s$ (figure 8a).

Conversely, for polymeric liquids with short relaxation times ( $De \ll 1$ ), we find that both transitions are $Oh_s$ -independent and occur at constant polymeric Ohnesorge number $Oh_p = Ec \times De$ :

1. (i) the dropping transition occurs at $Oh_{p,d} \approx 0.048$ (figure 7); and

2. (ii) the jetting transition occurs at $Oh_{p,j} \approx 0.129$ (figure 8 b).

These behaviours reflect fundamentally different physical mechanisms: at high $De$ , depending on $Oh_s$ , the medium behaves like an elastic solid ( $Oh_s \to 0$ ) or Kelvin–Voigt solid (finite $Oh_s$ ). However, at low $De$ , polymer addition manifests as an enhanced viscous effect characterised by $Oh_p$ . The trend of dropping transition in the small $De$ regime is qualitatively similar to recently reported experimental observation (Cabalgante-Corrales et al. Reference Cabalgante-Corrales, Muñoz-Sánchez, López-Herrera, Cabezas, Vega and Montanero2025). Although, a quantitative comparison cannot be made due to significant differences in $Bo$ . We further investigate the jetting transition using the slender jet equations in § 4.2 following similar approaches by Driessen et al. (Reference Driessen, Jeurissen, Wijshoff, Toschi and Lohse2013), Gordillo, Onuki & Tagawa (Reference Gordillo, Onuki and Tagawa2020), Zinelis et al. (Reference Zinelis, Abadie, McKinley and Matar2024) and Sen et al. (Reference Sen, Lohse and Jalaal2024).

4.2. What sets the different transitions, and what do we learn from these transitions?

To understand the mechanisms governing bubble cavity collapse, we analyse the jet dynamics using a control volume approach (figure 9). Employing the slender jet approximation (Shi, Brenner & Nagel Reference Shi, Brenner and Nagel1994; Driessen et al. Reference Driessen, Jeurissen, Wijshoff, Toschi and Lohse2013; Eggers & Fontelos Reference Eggers and Fontelos2015), given the small radial-to-axial length scale ratio, the vertical momentum equation for the jet reads

(4.1) \begin{align} \rho _s \left ( \frac {\partial v}{\partial t} + v\frac {\partial v}{\partial z} \right ) = - \gamma \frac {\partial \kappa } {\partial z} + \frac {1}{h^2} \frac {\partial }{\partial z}\left [ h^2 \left ( 3\eta _s \frac {\partial v}{\partial z} + G\left (\mathcal {A}_{zz}-1\right ) \right ) \right ]. \end{align}

Figure 9.

Temporal evolution of the Worthington jet for a representative case, where the jet emerges, reaches a maximum, and is pulled to merge with the liquid bath. The control volume contains the jet region, as shown by the region within the grey lines. Here, $h(z,t)$ is the width of the jet, which becomes $h_{{base}}$ at the base of the jet. The capillary force at the jet base is $F_\gamma = \gamma (2 \pi h_{{base}})$ that acts radially outwards. At the same time, the elastic and viscous stresses act at the base of the jet as $F_\eta + F_p = ( \sigma _{\eta , {base}} + \sigma _{p, {base}}) \pi h_{{base}}^2$ .

Here, $v(z,t)$ is the radially averaged jet velocity, and the shape of this jet is $h(z,t)$ . We define a control volume containing the emerging jet that is always bounded by the inflection points at the interfaces, see figure 9(b). Integrating over this control volume (with differential volume element $d\Omega = \pi h(z,t)^2\mathrm {d}z$ ) yields the force balance (Trouton Reference Trouton1906)

(4.2) \begin{align} \frac {\textrm {d} \mathcal {M}_{{jet}}}{\textrm {d} t} = 3 \eta _sh^2\frac {\partial v}{\partial z}\Bigg |_{{base}} + Gh^2(\mathcal {A}_{zz}-1)\Bigg |_{{base}} = \left (\sigma _{\eta ,{base}} + \sigma _{p,{base}}\right )\pi h_{{base}}^2 ,\end{align}

where $\mathcal {M}_{{jet}}(t) = \int _{\Omega (t)}\rho _sv(z,t)\pi h(z,t)^2\mathrm {d}z$ denotes the momentum of the jet. The capillary stress (first term on the right-hand side of (4.1)) integral vanishes due to orthogonal interface intersection with the control volume (see Marchand et al. Reference Marchand, Weijs, Snoeijer and Andreotti2011; Munro Reference Munro2019, pp. 16–21). We chose this control volume because of its vanishing integral feature. Furthermore, the integral of the second term on the right-hand side forms an exact integral which vanishes at the tip where it is zero owing to $h(z = L_{{max}}(t)) = 0$ . Consequently, jet evolution depends solely on stresses at the base: viscous ( $\sigma _{\eta ,{base}}(t)$ ) and elastic ( $\sigma _{p,{base}}(t)$ ). For relevant $Oh_s$ values, $\sigma _{\eta ,{base}}(t)$ is too weak to suppress the Worthington jet. Numerical simulations allow us to estimate $\sigma _{p,{base}}(t)$ . As the capillary waves collapse, the base elastic stress reaches a global maximum, before decreasing again at later times. Jet formation occurs if inertial flow focusing is sufficiently strong at the peak elastic stress. We will now evaluate this competition for the two limits of $De$ .

4.2.1. The limit of $De \to \infty$

Figure 10(a) shows that, for $De \gt 1$ , the maximum elastic stress $\text {max} (\sigma _{p,{base}}(t) )$ reaches a plateau, dependent only on $Ec$ . This $De$ -independence coincides with the extent of infinite $De$ asymptotes featured in the transitions discussed in § 4.1. The upper limit of elastic resistance competes with inertial flow focusing to inhibit jet formation. We quantify the inertial stresses at peak elastic stress using

(4.3) \begin{align} \sigma _{I,{base}} = \frac {2}{h_{{base}}^2}\int _{o}^{h_{{base}}} \rho _s v^2 h\mathrm {d}h, \end{align}

where $h_{{base}}$ is the jet width at its base (see figure 9). Figure 10(b) reveals that the ratio of elastic to inertial stresses is largely independent of $Oh_s$ . As $Ec$ increases, this ratio reaches a maximum beyond which jet suppression occurs. It is important to note that the apparent decrease in this stress ratio with increasing $Ec$ and $Oh_s$ near the jetting transition in figure 10(b) occurs due to a decrease in both the polymeric and inertial stresses in this region of the parameter space.

Figure 10.

(a) Evolution of the maximum elastic stress at the jet base ( $\text {max}(\sigma _{p,{base}}(t))$ ), normalised by the Laplace pressure scale $\sigma _\gamma = \gamma /R_0$ , as a function of $De$ for different $Ec$ at $Oh_s = 0.001$ . Note that $Oh_p = Ec\times De$ . (b) Comparison of the resistive elastic stress $\text {max}(\sigma _{p,{base}}(t))$ in the high $De$ regime $(\to \infty )$ against the inertial stresses $\sigma _{I,{base}}$ , plotted against $Ec$ for different $Oh_s$ .

4.2.2. The limit of $De \to 0$

In the zero $De$ limit, polymeric liquids exhibit additional viscous effects characterised by the polymeric Ohnesorge number $Oh_p$ (also see § 4.1). The maximum elastic stress $\text {max}(\sigma _{p,{base}}(t))$ , when normalised by the Newtonian-like viscous stress $\sigma _{N,{base}}$ , collapses for all $Oh_p$ as $De \to 0$ , where

(4.4) \begin{align} \sigma _{N,{base}} = \frac {2}{h_{{base}}^2}\int _{o}^{h_{{base}}} G\lambda \frac {\partial v}{\partial z} h\mathrm {d}h. \end{align}

As $De$ approaches unity, marking the onset of the infinite $De$ asymptotic regime, the elastic stress scales as $\text {max}(\sigma _{p,{base}}(t)) \sim De \times \sigma _{N,{base}}$ . This scaling remarkably resembles that predicted by Boyko et al. (Reference Boyko, Hinch and Stone2024) for flow in a slowly varying contraction at the infinite $De$ asymptote, despite significant geometric differences. While our study focuses on free surface flows and Boyko et al. (Reference Boyko, Hinch and Stone2024) examined contraction geometries, this unexpected similarity hints at a potentially universal behaviour near the infinite $De$ asymptote. To further examine this intriguing connection, a similar closed-form $De$ expansion for free surface flows is necessary. However, we caution that this scaling approach to the infinite $De$ asymptote could be system-dependent (Hinch et al. Reference Hinch, Boyko and Stone2024).

At zero $De$ , the elastic stress reduces to a Newtonian-like viscous stress with polymeric viscosity $\eta _p$ , yielding $\sigma _p \approx 2G\lambda \boldsymbol {\mathcal {D}}$ for the Oldroyd-B rheology. The force balance in (4.2) becomes

(4.5) \begin{align} \frac {d \mathcal {M}_{{jet}}}{d t} = \left (3\eta _s + 2G\lambda \right )h^2\frac {\partial v}{\partial z}\Bigg |_{{base}}, \end{align}

which depicts the balance of jet inertia with viscous forces. Using characteristic scales for jet momentum $\mathcal {M}_{{jet}} \sim \rho V_{{jet}} h_{{base}}^3$ , velocity gradient $\partial _zv \sim V_{{jet}}/\delta _\eta$ and time $\tau _i \sim h_{{base}}/V_{{jet}}$ , the force balance yields

(4.6) \begin{align} \rho V_{{jet}}^2 \sim \eta _{{effective}}\frac {V_{{jet}}}{\delta _\eta }. \end{align}

Here, $\delta _\eta$ represents the viscous length scale and the effective viscosity is

(4.7) \begin{align} \eta _{{effective}} = 3\eta _s + 2G\lambda . \end{align}

Since polymers do not affect the flow before jet formation (§ 3), the jet Weber number remains constant at inception (Blanco-Rodríguez & Gordillo Reference Blanco-Rodríguez and Gordillo2021)

(4.8) \begin{align} We_{{jet}} = \frac {\rho V_{{jet}}^2 \delta _\eta }{\gamma } =\,\text {constant}. \end{align}

Combining (4.6) and (4.8), we get

(4.9) \begin{align} V_{{jet}} \sim \frac {\gamma }{\eta _{{effective}}} ,\end{align}

analogous to Newtonian media but with modified viscosity (Gordillo & Rodríguez-Rodríguez Reference Gordillo and Rodríguez-Rodríguez2019; Blanco-Rodríguez & Gordillo Reference Blanco-Rodríguez and Gordillo2020).

Figure 11(b) illustrates the jet velocity as a function of the effective Ohnesorge number

(4.10) \begin{align} Oh_{{effective}} = 3Oh_s+2Oh_p ,\end{align}

(reflecting (4.7)) at different $De$ . We stress that the jet velocity varies in time (Deike et al. Reference Deike, Ghabache, Liger-Belair, Das, Zaleski, Popinet and Séon2018; Sanjay et al. Reference Sanjay, Sen, Kant and Lohse2022; Gordillo & Blanco-Rodríguez Reference Gordillo and Blanco-Rodríguez2023) and is maximum at its inception, which is the value that we report here. For sufficiently large $Oh_{{effective}}$ and small $De$ , we recover the scaling predicted in (4.9). However, as $De$ increases, the added elastic stresses cannot be directly substituted with Newtonian-like viscous stresses, and the underlying assumption fails, evident in the deviation of $V_{{jet}}$ from the prediction.

Figure 11.

(a) Evolution of the maximum elastic stress at the jet base ( $\text {max}(\sigma _{p,{base}}(t))$ ), normalised by the Newtonian-like viscous stress $\sigma _{N,{base}}$ with viscosity $\eta _{p} = G\lambda$ , as a function of $De$ for different $Oh_p$ at $Oh_s = 0.001$ . The grey dashed horizontal line represents $\text {max}(\sigma _{p,{base}}(t)) \approx \sigma _{N,{base}}$ while the black dashed line serves as a guide to the eye representing $\text {max}(\sigma _{p,{base}}(t))/\sigma _{N,{base}} \sim De$ . Note that $Ec = Oh_p/De$ . (b) The variation of jet’s tip velocity $V_{{jet}}$ , normalised by the inertiocapillary velocity $u_\gamma = \sqrt {\gamma /\rho _sR_0}$ , with $Oh_{{effective}} = 3Oh_s + 2 Oh_p$ at different $De$ and $Oh_s = 0.01$ . The grey dashed line represents $V_{{jet}} \sim \gamma /\eta _{{effective}}$ .

On the other hand, for small $Oh_{{effective}}$ , $V_{{jet}}$ for all $De$ closely matches the corresponding speed in Newtonian liquids, as observed in figure 11(b) for $Oh_p = 0$ . As $Oh_{p}$ increases, $V_{{jet}}$ also increases, reaching a maximum before decreasing and following (4.9). Although the capillary wave speed remains unaffected in the polymeric medium, increasing $Oh_p$ triggers elastic stresses in smaller-wavelength capillary waves, which are promptly dissipated due to small $De$ . Consequently, improved flow focusing occurs as the strongest undamped capillary wave survives, thus increasing $V_{{jet}}$ . This behaviour is analogous to the non-monotonicity observed and well understood for Newtonian liquids at small $Oh_s$ (Duchemin et al. Reference Duchemin, Popinet, Josserand and Zaleski2002; Deike et al. Reference Deike, Ghabache, Liger-Belair, Das, Zaleski, Popinet and Séon2018; Gordillo & Rodríguez-Rodríguez Reference Gordillo and Rodríguez-Rodríguez2019; Yang et al. Reference Yang, Tian and Thoroddsen2020; Sanjay et al. Reference Sanjay, Sen, Kant and Lohse2022; Gordillo & Blanco-Rodríguez Reference Gordillo and Blanco-Rodríguez2023), further supporting the observation that polymeric liquid exhibit a Newtonian-like viscous response in the zero $De$ limit.

To further quantify this behaviour, figure 12(a) illustrates jet features at inception for different $Oh_p$ at $De = 0.001$ , while figure 12(b) shows jet radius as a function of $Oh_p$ at different $De$ . At small $Oh_p$ , we observe that the jet radius maintains a value comparable to the Newtonian reference case (figure 12 a: $Oh_p = 0, 0.005$ ). This behaviour is consistent with the $De \to 0$ limit, where polymeric additives primarily contribute enhanced effective viscosity. Since the jet radius determines the resulting drop size (Gañán-Calvo Reference Gañán-Calvo2017; Blanco-Rodríguez & Gordillo Reference Blanco-Rodríguez and Gordillo2020), this independence of jet radii in the low $Oh_p$ regime suggests minimal variation in droplet size distribution compared with Newtonian cases. As $Oh_p$ increases, we observe a pronounced reduction in jet width until reaching $Oh_{p,c} \approx 0.017$ (figure 12 a: $Oh_p = 0.013, 0.018$ ). At this critical value, the system transitions to a bubble entrainment regime (figure 12 a: $Oh_p = 0.018, 0.025$ Gordillo & Rodríguez-Rodríguez Reference Gordillo and Rodríguez-Rodríguez2019; Blanco-Rodríguez & Gordillo Reference Blanco-Rodríguez and Gordillo2020; Rodríguez-Díaz et al. Reference Rodríguez-Díaz, Rubio, Montanero, Gañán-Calvo and Cabezas2023). Interestingly, the prediction for Newtonian liquids applies well to viscoelastic liquids by substituting $Oh_p$ for $Oh_s$ (figure 12 b), particularly in the $De \to 0$ limit. Beyond $Oh_{p,c}$ , the jet radius becomes ill defined as the jet gradually widens (figure 12 a: $Oh_p = 0.035$ ), first reaching the dropping transition at $Oh_{p,d} \approx 0.048$ (figures 7 and 12 b) and ultimately vanishing at $Oh_{p,j} \approx 0.129$ (figure 8 b).

Figure 12.

The capillary waves focus and collapse at the bottom of the cavity. (a) The inception of the jet after the collapse at different $Oh_p$ at $De = 0.001$ and $Oh_s = 0.01$ . The radius of the jet at the base $R_{{jet}}$ decreases with $Oh_p$ until $Oh_{p,c} = 0.017$ , beyond which bubbles are entrained and the jet radius increases. (b) Radius of jet $R_{{jet}}$ with $Oh_p$ at $Oh_s = 0.01$ and different $De$ . Here, $R_{{jet}}$ remains close to the value at the Newtonian limit $Oh_p = 0$ , and decreases sharply as it approaches $Oh_{p,c}$ . Beyond $Oh_{p,j}$ jets are no longer observed.