# The atomising pulsed jet

Yash Kulkarni, Cesar Pairetti, Raphaël Villiers, Stéphane Popinet, Stéphane Zaleski

► To cite this version:

Yash Kulkarni, Cesar Pairetti, Raphaël Villiers, Stéphane Popinet, Stéphane Zaleski. The atomising pulsed jet. Journal of Fluid Mechanics, 2025, 1009, pp.A35. 10.1017/jfm.2025.218. hal-05104157

HAL Id: hal-05104157

https://hal.science/hal-05104157v1

Submitted on 9 Jun 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_2.png" /></figure>

# The atomising pulsed jet

Yash Kulkarni<fnref n="1" /> ID, Cesar Pairetti<fnref n="1" />,<fnref n="2" /> ID, Raphaël Villiers<fnref n="1" />, Stéphane Popinet<fnref n="1" /> and Stéphane Zaleski<fnref n="1" />,<fnref n="2" /> ID

**Corresponding authors:** Yash Kulkarni, kulkarniyash2398@gmail.com; Cesar Pairetti, paire.cesar@gmail.com; Stéphane Zaleski, stephane.zaleski@sorbonne-universite.fr

(Received 19 May 2024; revised 20 December 2024; accepted 31 January 2025)

Direct numerical simulations of the injection of a pulsed round liquid jet in a stagnant gas are performed in a series of runs of geometrically progressing resolution. The Reynolds and Weber numbers and the density ratio are sufficiently large for reaching a complex high-speed atomisation regime but not so large so that the small length scales of the flow are impossible to resolve, except for a very small liquid-sheet thickness. The Weber number based on grid size is then small, an indication that the simulations are very well resolved. Computations are performed using octree adaptive mesh refinement with a finite volume method and height-function computation of curvature, down to a specified minimum grid size $\Delta$. Qualitative analysis of the flow and its topology reveals a complex structure of ligaments, sheets, droplets and bubbles that evolve and interact through impacts, ligament breakup, sheet rupture and engulfment of air bubbles in the liquid. A rich gallery of images of entangled structures is produced. Most processes occurring in this type of atomisation are reproduced in detail, except at the instant of thin sheet perforation or breakup. We analyse droplet statistics, showing that as the grid resolution is increased, the small-scale part of the distribution does not converge, and contains a large number of droplets close in order of magnitude to the minimum grid size with a significant peak at $d = 3\Delta$. This non-convergence arises from the numerical sheet breakup effect, in which the interface becomes rough just before it breaks. The rough appearance of the interface is associated with a high-wavenumber oscillation of the curvature. To recover convergence, we apply the controlled 'manifold death' numerical procedure, in which thin sheets are detected, and then pierced by fiat before they reach a set critical thickness $h_c$ that is always larger than $6\Delta$. This allows convergence of the droplet frequency above a certain critical diameter $d_c$, above and close to $h_c$. A unimodal distribution is observed in the converged range. The number of holes pierced in the sheet is a free parameter in the manifold death procedure;

<footnote n="1">¹Sorbonne Université and CNRS, UMR 7190, Institut Jean Le Rond d'Alembert, 75005 Paris, France</footnote>

<footnote n="2">²Facultad de Ciencias Exactas, Ingeniería y Agrimensura, Universidad Nacional de Rosario, 2000 Rosario, Argentina</footnote>
however, we use the Kibble-Zurek theory to predict the number of holes expected on heuristic physical grounds.

**Key words:** aerosols/atomisation, breakup/coalescence, bubble dynamics

# 1. Introduction

Atomisation simulations have progressed at an amazing rate, however, the topic is still far from mature. As we show in this paper, the prediction of the breakup of liquid masses in a typical large-speed flow is marred by vexing numerical effects and profound physical uncertainties. New developments, such as better codes, rapidly increasing processing power and new numerical methods are however poised to mitigate the difficulties. In this paper we investigate the pulsed jet, a paradigmatic case of atomising flow inspired by diesel engine jets, although we do not aim at solving a particular applied problem or even any experimental configuration, but rather to investigate the potential and limitations of direct numerical simulation (DNS) of an atomising flow, and in particular, the ability of such DNS to reveal physically relevant properties of the flow through statistically converged numerical approaches. The emphasis here is on statistical convergence, rather than ‘trajectory’ convergence, since only the former is conceivable in very complex, irregular flows.

The study of such convergence already has a long history, despite rapid progress. We review specifically the history of round-jet atomising simulations and we refer the reader to reviews on the general topic of atomisation and its simulation (Gorokhovski & Herrmann 2008; Villermaux 2020). The first attempts at testing the convergence of the probability distribution function (PDF) of droplet sizes were those of Herrmann (2011). His simulations of a round jet used a variant of the level-set method and showed that the number of small droplets in the PDF was underestimated by coarse-grid simulations. This can be understood by the fact that the level-set methods tend to eliminate small droplets by ‘evaporating’ them. On the other hand, volume-of-fluid (VOF) methods keep too many droplets in an intermediate range around the grid size $\Delta$, and in that range the PDF of droplet sizes is overestimated. This can be seen, for example, in a study of the convergence of the droplet-size distribution in the round jet by Pairetti et al. (2020) where coarse grids overestimate the number of large droplets.

A graphical illustration of the contrasting effects of level-set and VOF methods on the droplet sizes in atomisation is shown on figure 1. Despite the inaccuracy of the distribution for small droplet diameters, both methods may converge as the grid is refined progressively. There is thus some hope that very large simulations will eventually produce converged distributions, a hope we want to explore in this paper.

Aside from the two investigations cited above, very few studies address the issue of PDF convergence, despite detailed analyses of aspects of the flow. Perhaps the first round, single-jet atomisation simulation in the conditions of a diesel jet was that of Bianchi et al. (2007) using the VOF method, followed by more detailed simulations by Ménard et al. (2007) using the combined level-set VOF method (CLSVOF). This was followed by other simulations using CLSVOF by Lebas et al. (2009), Chesnel et al. (2011) and Anez et al. (2019) and by simulations using the VOF method (Fuster et al. 2009). The latter studied both the coaxial jet cases in the so-called ‘assisted atomisation’ set-up and the conical jet. For the conical jet, Fuster et al. (2009) showed that the PDF changed drastically with
Journal of Fluid Mechanics

<figure><img src="image_3.png" /><figcaption>Figure 1. An illustration of the outcomes for the numerical simulation of a thinning liquid sheet; panel (a) shows the initial configuration, before breakup. The other three images (b–d) show a schematic view of the outcome either in reality or in various types of numerical simulation. It is arbitrarily assumed that there are two breakup locations. Both the VOF and the level-set methods yield topology changes when the sheet thickness reaches the grid size. In (b) fragments larger than the grid size are obtained because of mass conservation in the VOF method. (c) In reality, the sheet thinning continues until much later than in the numerics, unless extremely fine grids are used. The final size of some of the droplets is then much smaller than in the VOF simulation. (d) The level-set or diffuse-interface methods on the other hand evaporate the thin parts of the sheet and loses much more mass.</figcaption></figure>

grid resolution. The PDF never peaked as the droplet size was decreased towards the grid
size. Until 2022, the most detailed (or highest-fidelity) published simulations of a round jet
were performed by Shinjo & Umemura (2010) using a combination of level-set and VOF
techniques with a ratio of jet diameter to grid size of 286. These simulations, together
with those of Herrmann (2011) are the two outliers above the trend line in figure 2. In
addition to revealing a wealth of details about ligament and sheet formation, and the
perforation of sheets, Shinjo & Umemura have also shown distributions of droplet sizes,
but without investigating explicitly their convergence. As in earlier simulations of the
round jet, the droplet-size distribution has a single maximum close to the small droplet
end of the spectrum, almost at the grid size. Studies by Jarrahbashi & Sirignano (2014)
and Jarrahbashi et al. (2016) focused on a round, spatially periodic jet and analysed the
effect of vortex dynamics on the development of instabilities along with the jet core.
Studies by Zhang et al. (2020) and Pairetti et al. (2020) applied the VOF method with
octree adaptive mesh refinement. As already mentioned, the latter paper investigated the
distribution of droplet sizes showing some difficulties in reaching statistical convergence.
Several authors including Torregrosa et al. (2020) and Salvador et al. (2018) have focused
on the effect of injection conditions and turbulence. Khanwale et al. (2022) and Saurabh
et al. (2023), using a diffusive interface method with octree adaptive mesh refinement,
obtained the most refined simulations published at the time of writing with a special
Y. Kulkarni, C. Pairetti, R. Villiers, S. Popinet and S. Zaleski

<figure><img src="image_4.png" /><figcaption>Figure 2. The increase in two-phase round-jet grid resolution in time. The graph includes two simulations published only on the Gerris and Basilisk websites (and in other channels outside of academic journals) before 2017. The 2024 simulations are those reported in this paper.</figcaption></figure>

treatment of the refinement in thin regions. It should be noted that all of the round-jet
studies cited above involve a variety of physical parameters and grid resolutions, despite
the fact that they share similarities, such as moderate density ratios and liquid Reynolds
numbers (defined below) in a range $5000 \le Re_l \le 13\,400$ that are attainable by DNS for
single-phase flow. For example, the two octree studies of Zhang *et al.* (2020) and Pairetti
*et al.* (2020) have gas-based Weber numbers (defined below) of respectively 177 and 417.
As a general rule, it is better to select relatively small density ratios, Reynolds and Weber
numbers to increase the likelihood of reaching convergence in computation. In this work,
we thus chose the parameters in the lower range of *Re*, *We*, as it correspondens to Weber
number and density ratio for this type of flow. Similar attempts at obtaining convergence
were realised in the simulations of Ling *et al.* (2017), investigating assisted quasi-planar
atomisation after the experiment of Grenoble (Ben Rayana *et al.* 2006; Fuster *et al.* 2013).
Ling *et al.* (2017) performed computations on four grids of increasing resolution. Despite
the huge computational effort involved in that latter simulation, convergence was clearly
not reached, both quantitatively as droplet-size distributions had a very narrow region of
overlap, and qualitatively as clearly under-resolved structures in thin sheet-like regions
were observed. One of our objectives in this work is to see what happens if ever finer grids
are used, until very thin sheets are clearly resolved. In experiments as well as in some
simulations it is clear that thin sheets are the site of weak spots (Lohse & Villermaux
2020). The presence of these weak spots as one of the mechanisms leading to atomisation
forces the following somewhat sobering conclusion. Numerical simulations of diffuse-
interface, level-set or VOF type can never be converged if thin sheets break or perforate
only when they reach the grid size. Instead, a physical mechanism for weak spots or
perforations must be present to nucleate holes in a manner independent of the numerics.
Although such a mechanism may not be known yet, we suggest as a backstop procedure
to define a critical sheet thickness $h_c$ beyond which perforations occur with relatively high
frequency, that is, at the rate of one or more perforations per connected thin sheet region.
Chirco *et al.* (2022) have called manifold death (MD) this procedure (this choice of words
partially echoes that of Debrégeas et al. 1998). It should be noted that bags or membranes are observed at relatively low velocities in many atomisation configurations (droplets, jets etc.) and that at higher velocities, bags are not visible and perhaps non-existent. We note that the literature (Lasheras & Hopfinger 2000; Tolfts et al. 2023) distinguishes two kinds of atomisation, a bag or membrane-type atomisation at low Weber numbers and a fibre-type atomisation at high Weber numbers. The transition between the two behaviours was found by Tolfts et al. (2023) to lie between $We_g = 70$ and $We_g = 400$. From the connection between thin sheets and the lack of statistical convergence one may conjecture that in the fibre-type regime statistical convergence could be achieved. This is however contradicted by the lack of convergence observed by Herrmann (2011) at $We_g = 500$. A definite answer to this argument may have to wait for a study similar to the one in this paper in the fibre-type regime, a study which would probably be even more expensive.

In what follows we first describe the characteristics of our test case. We then describe our numerical method, including the MD procedure. Then we continue with our results that are both qualitative and quantitative. Qualitatively our results reveal new processes of numerical rupture. Quantitatively they allow an analysis of the droplet-size distributions or PDF. A discussion of the phenomenology of hole expansion in a thin sheet is given in Appendix A. We end with a conclusion including perspectives and discussion.

## 2. Mathematical model and numerical method

Our mathematical model is based on the mass and momentum conservation equations for incompressible and isothermal flow, i.e.

$$ \nabla \cdot \mathbf{u} = 0, \qquad (2.1) $$

$$ \frac{\partial \rho \mathbf{u}}{\partial t} + \nabla \cdot (\rho \mathbf{u} \mathbf{u}) = -\nabla p + \nabla \cdot (2\mu \mathbf{D}) + \mathbf{f}_{\sigma}, \qquad (2.2) $$

where $\mathbf{u}(\mathbf{x}, t)$ is the velocity field and $p(\mathbf{x}, t)$ is the pressure field. The tensor $\mathbf{D}$ is defined as $(1/2)[\nabla \mathbf{u} + (\nabla \mathbf{u})^T]$. The density and viscosity of the flow are noted as $\rho$ and $\mu$, respectively. The last term on the right-hand side of the Navier–Stokes equation (2.2) represents the surface tension force

$$ \mathbf{f}_{\sigma} = \sigma \kappa \mathbf{n} \delta_{s}, \qquad (2.3) $$

which depends on the surface tension $\sigma$ and the interface shape, particularly on its curvature $\kappa$ and unit normal vector $\mathbf{n}$. The Dirac distribution $\delta_s$ indicates that the force only acts at the free surface. We consider the phase distribution function $c(\mathbf{x}, t)$ that takes value unity in the reference phase and null outside of it. The transport equation for this function is

$$ \frac{\partial c}{\partial t} + \nabla \cdot (c \mathbf{u}) = c \nabla \cdot \mathbf{u}. \qquad (2.4) $$

In the context of incompressible flow, the right-hand side above is equal to zero. The $c$ function implicitly defines the interface at its discontinuity surface, defining also the $\delta_s$, $\mathbf{n}$ and $\kappa$ fields, e.g. $\nabla c = \mathbf{n} \delta_s$.

The VOF method represents the evolution of $c$ using the piecewise linear interface capturing method of Hirt & Nichols (1981) and DeBar (1974). In this context, the mean value of the colour function on a cell is

$$ f = \frac{1}{|\Omega|} \int_{\Omega} c(\mathbf{x}, t) dV, \qquad (2.5) $$
where $|\Omega|$ is the volume of the cell $\Omega$. Then $f$ is the volume fraction of the reference phase in the cell. The mixture properties of the cell may then be computed by arithmetic averages

$$ \rho_{\Omega} = f\rho_{l} + (1 - f)\rho_{g}, \quad \mu_{\Omega} = f\mu_{l} + (1 - f)\mu_{g}. \qquad (2.6) $$

Spatial discretisation of the above equation is realised on a network of cubic cells obtained by a tree-like subdivision of an initial cubic cell of size $L_0$. The subdivision is realised using a wavelet-based error estimate and is adapted dynamically as the simulation progresses. At all times the maximum subdivision level is a fixed number $\ell$ and the smallest grid size is

$$ \Delta_{\ell} = 2^{-\ell} L_0. \qquad (2.7) $$

From now on, we consider the algebraic equations derived from applying the finite volume method on each cell. In this context, the approximate projection method by Chorin (1968) can be used to solve the coupling between (2.1) and (2.2), considering that the velocity $\mathbf{u}$ is staggered in time with respect to the volume fraction $f$ and the pressure $p$. The discrete set of equations can then be expressed as

$$ \frac{\rho \mathbf{u}^* - \rho \mathbf{u}^n}{\Delta t} + \nabla \cdot \left( \rho^{n+\frac{1}{2}} \mathbf{u}^{n+1/2} \mathbf{u}^{n+1/2} \right) = \nabla \cdot \left[ \mu^{n+\frac{1}{2}} (\mathbf{D}^n + \mathbf{D}^*) \right] + \mathbf{f}_{\sigma}^{n+1/2}, \quad (2.8) $$

$$ \nabla \cdot \left( \frac{\Delta t}{\rho^{n+\frac{1}{2}}} \nabla p^{n+\frac{1}{2}} \right) = \nabla \cdot \mathbf{u}^*, \qquad (2.9) $$

$$ \mathbf{u}^{n+1} = \mathbf{u}^* - \frac{\Delta t}{\rho^{n+\frac{1}{2}}} \nabla p^{n+\frac{1}{2}}. \qquad (2.10) $$

In the above the expression, $\nabla \cdot (\rho^{n+\frac{1}{2}} \mathbf{u}^{n+1/2} \mathbf{u}^{n+1/2})$ must be interpreted as the use of a predictor-corrector scheme for advection of the velocity field. Both the predictor and the corrector use the Bell-Collela-Glaz scheme (Popinet 2003, 2009), and involve two projections. A full description of the predictor corrector scheme with projection is best found in the 'literate' source code at http://basilisk.fr/src/navier-stokes/centered.h. In addition to these equations, $f$ is advanced in time using the split-volume-fraction advection scheme of Weymouth & Yue (2010). We use the semi-implicit Crank-Nicholson time stepping to compute the diffusive flux due to the viscous term in (2.8). An important aspect of the method is the numerical approximation of the surface tension force $\mathbf{f}_\sigma$. We use the continuous surface force well-balanced formulation $\mathbf{f}_\sigma = \sigma \kappa \nabla c$ for (2.3), where we discretise $\nabla c$ at cell faces with the same scheme employed to compute the pressure gradient $\nabla p$. This is useful to reduce spurious currents provided curvature is accurately computed, as explained by Popinet (2009, 2018). The interface curvature is computed using a complex method described in Popinet (2009). When the grid resolution is sufficient, second-order stencils based on height functions are used to compute curvature. In some large-curvature cases the height functions are not defined and fits by quadratic functions are used instead. Equations (2.9) and (2.10) are the projection steps that will ensure continuity for the velocity field at time step $n+1$. The mesh is adapted dynamically, by splitting the grid cells to eight smaller cells whenever the estimated local discretisation error exceeds a set threshold (see Appendix B). Details of the procedure can be found on the web site basilisk.fr and, in particular, in the documented free code available at http://basilisk.fr/src/examples/atomisation.c.

In all the data reported in this paper, the velocities are expressed in multiples of the jet velocity $U$. All the numerical values for time are made dimensionless by
$\tau_0 = L_0/(3U)$ and all the reported lengths are made dimensionless by $L_0/3$. However in the theoretical discussions the dimensional values are used unless otherwise indicated. To avoid confusion, we note with a * exponent the dimensionless values, so that

$$t^* = t/\tau_0 \qquad x^* = 3x/L_0 \qquad \Delta^* = 3(2^{-\ell}). \qquad (2.11)$$

In addition to the aforementioned model and procedure, we implement the MD procedure as described by Chirco et al. (2022). This procedure involves the following two primary steps. (i) The identification of thin liquid sheets with a thickness of $h_c$ or less, achieved through a local integration of the characteristic function $c$. This integration results in the calculation of the signature of a quadratic form, based on the local main inertia moments. (ii) The creation of holes within the so identified thin sheets, using a probabilistic approach. This procedure is applied at a given frequency, defined by the user within this method. From a physical standpoint, this technique aligns with the weak spot model, wherein the likelihood of hole formation is exceedingly low for $h > h_c$ and substantially increases, rapidly approaching certainty, for $h < h_c$. Note that this procedure distinguishes thin sheets from other shapes such as ligaments, only the thin sheets are detected and perforated. The specific implementation of the MD method we used is called the signature method. The complete MD code with the signature method can be accessed at http://basilisk.fr/sandbox/lchirco/signature.h. In the present case, hole punching is attempted every time interval of $\tau_m = 0.01\tau_0$, that is, the raw dimensionless breakup frequency is $\tilde{f}_b^* = \tau_0/\tau_m = 100$. At these instants $N_b = 200$ holes are attempted. We loosely term the pair of parameters $(\tilde{f}_b^*, N_b)$ the breakup frequency. The breakup frequency chosen in this study is probably too high, as will appear in the discussion below. However, it is chosen to maintain the same parameters as in other studies, and possible improvements are discussed in the conclusion section. As we show in Appendix A, the characteristic time for sheet thinning is $\tau_c$ given by

$$\tau_c = \frac{D}{U} \left( \frac{\rho_l}{\rho_g} \right)^{1/2}, \qquad (2.12)$$

where $D$ is the jet diameter. From the definition of these quantities, we can obtain the breakup attempt frequency in units of the sheet thinning characteristic time $f_b = \tau_c/\tau_m$, i.e.

$$f_b = \frac{3\tilde{f}_b^* D}{L_0} \left( \frac{\rho_l}{\rho_g} \right)^{1/2}, \qquad (2.13)$$

that is, with the parameters of our simulation

$$f_b = \frac{50}{3} \left( \frac{\rho_l}{\rho_g} \right)^{1/2}. \qquad (2.14)$$

The numerical value is $f_b \simeq 87.9$, ensuring that punching attempts are sufficiently frequent as not to 'miss' transient thinning sheets. It is large enough to eliminate numerical breakup events but is likely to be too high as discussed below.

Using this numerical method, we analyse the atomisation of a circular jet injected at average velocity $U$ in a gas-filled cubical chamber of edge length $L_0$. We consider incompressible, isothermal flow. The dimensionless groups most relevant for this problem are the gas-based Weber number, the liquid-based Reynolds number and the ratios

$$We_g = \frac{\rho_g U^2 D}{\sigma}, \qquad Re_l = \frac{\rho_l U D}{\mu_l}, \qquad \rho^* = \frac{\rho_l}{\rho_g}, \qquad \mu^* = \frac{\mu_l}{\mu_g}, \qquad \frac{L_0}{D}, \qquad (2.15)$$
<table><thead><tr><th><i>Re</i><sub>l</sub></th><th><i>We</i><sub>g</sub></th><th>&rho;<sup>*</sup></th><th>&mu;<sup>*</sup></th><th><i>A</i><sub>p</sub></th><th>&omega;<sub>p</sub></th><th><i>L</i><sub>0</sub>/<i>D</i></th></tr></thead><tbody><tr><td>5800</td><td>200</td><td>27.84</td><td>27.84</td><td>0.05</td><td>&pi;/5</td><td>18</td></tr></tbody></table>
<p>Table 1. The dimensionless numbers characterising our simulation. The Reynolds number <i>Re</i><sub>l</sub> based on the liquid is rather moderate. The density and viscosity ratios are identical, which implies that <i>Re</i><sub>g</sub> = <i>Re</i><sub>l</sub>.</p>

where $\rho_l$ and $\rho_g$ are the densities of liquid and gas, respectively. The viscosity coefficients are noted by $\mu_l$ and $\mu_g$, and $\sigma$ is the surface tension coefficient. The boundary condition on the injection plane $x=0$ imposes the no-slip condition everywhere except on the liquid section ($y^2 + z^2 < D^2/4$), where the injection velocity is

$$u_x(t) = U \left[ 1 + A_p \sin \left( \omega_p U t / D \right) \right]. \qquad (2.16)$$

On the remaining cube sides, we allow free outflow: $\partial_n \mathbf{u}_\Gamma = 0$ and $p_\Gamma = 0$.

The dimensionless parameters of our simulation are given in **table 1**. These numbers are in the low range of the dimensionless numbers used in the above cited work. They are much smaller than the numbers of the spray G case of the engine combustion network (see Duke *et al.* 2017) used by Zhang *et al.* (2020). The Reynolds number and density ratio are identical to those of Ménard *et al.* (2007) and of the same order as those of Herrmann (2011). The equality of the gas and liquid Reynolds numbers stems from the fact that the viscosity and density ratios are equal, so the ratio of kinematic viscosities is one. However, while the moderate density ratio is realistic in the context of diesel jets, the viscosity ratio is not and would be much larger in a realistic set-up, which would lead to much larger values of $Re_g$ that are expensive to resolve by DNS. This motivates the choice of viscosity ratio. We note that the same choice of equal viscosity and density ratios and moderate $We_g$ was made by Ling *et al.* (2017) who also studied atomisation and PDF convergence in a simplified set-up. The Reynolds number is somewhat larger than that of Shinjo & Umemura (2010) while the gas Weber number $We_g$ is smaller than that of all the other papers in the literature. The rationale for selecting such a small Weber number is twofold. First the simplest prediction for the droplet size, based on the ratio of Bernoulli pressure $\rho_g U^2$ to surface tension, is $d = DWe_g^{-1}$. Indeed if a purely inviscid flow with a velocity jump is considered, the cutoff wavelength for the Kelvin-Helmholtz instability scales as $DWe_g^{-1}$. The wavelength thus obtained is a natural scale for the subsequent formation of sheets and ligaments. The number of grid points in the diameter of a droplet of this diameter $d$ is then $We_\Delta^{-1}$, where

$$We_l = \rho_g U^2 l / \sigma \qquad (2.17)$$

is the Weber number based on length scale $l$ and $\Delta$ is the grid size. Similarly,

$$Oh_l = \mu_l (\sigma \rho_l l)^{-1/2} \qquad (2.18)$$

is the Ohnesorge number based on scale $l$. The values of $We_\Delta$ for two of the grids we use are given in **table 2**. We also give the value of the Ohnesorge numbers $Oh_\Delta$ and $Oh_{hc}$, which characterise the interplay of viscous and surface tension effects. In particular, they typify the regime of the Taylor-Culick rims (Song & Tryggvason 1999; Savva & Bush 2009) occurring at the edge of a sheet of thickness $\Delta$ or $h_c$. The Oh number of the thin sheets at the moment of hole formation is relatively large, implying a viscous Taylor-Culick flow and the absence of unsteady fragmentation as in Wang & Bourouiba (2018) and Kant *et al.* (2023).
<table><caption>Table 2. Characteristic numbers related to the grid size. The level $ℓ$ is defined in (2.7).</caption><thead><tr><th>Level <i>ℓ</i></th><th><i>D</i>/<i>Δ</i></th><th><i>We</i><sub>Δ</sub></th><th><i>Oh</i><sub>Δ</sub></th></tr></thead><tbody><tr><td>11</td><td>114</td><td>1.75</td><td>0.137</td></tr><tr><td>15</td><td>1820</td><td>0.11</td><td>0.55</td></tr></tbody></table>

<table><caption>Table 3. Characteristic numbers based on the critical sheet thickness in our simulations. The MD level $m$ is defined in the text.</caption><thead><tr><th>MD level <i>m</i></th><th><i>D</i>/<i>h</i><sub>c</sub></th><th><i>We</i><sub><i>h</i><sub>c</sub></sub></th><th><i>Oh</i><sub><i>h</i><sub>c</sub></sub></th></tr></thead><tbody><tr><td>12</td><td>227</td><td>0.88</td><td>0.194</td></tr><tr><td>13</td><td>455</td><td>0.44</td><td>0.27</td></tr></tbody></table>

The second important consequence of selecting the relatively small value $We_g = 200$ is that it makes the simulation sit on the boundary between membrane type and fibre type defined in Lasheras & Hopfinger (2000). Membrane type in that later paper refers to the formation of bags and sheets. Thus, the prevalence of thin sheets discussed in this paper may be less marked at higher Weber numbers, the sheets being replaced by fibres, and also less important at smaller Weber numbers where the flow may be less unstable.

## 3. Results

### 3.1. *Uncontrolled perforation: numerical sheet breakup*

Figure 3 shows a global view of the jet evolution (see also supplementary movies 1–6 available at https://doi.org/10.1017/jfm.2025.218). The dynamics is initialised with a slightly penetrated jet at $t = 0$. At $t = 0.1$, a mushroom-like structure starts rolling up. The annular structure stretches and extends in the form of a ‘corona flap’ that develops behind the mushroom head. A rim forms at the edge of this flap then detaches. It is visible at $t = 0.2$ as a corrugated and partially fragmented annular ligament. As the jet evolves, pulsation results in a periodic series of corona flaps seen at $t = 0.8$. These flaps interact with the ligaments coming from the mushroom head. At $t = 2$, the jet has evolved further and we can identify an interval in the core structure where the effect of pulsation is lost. Several mechanisms of droplet formation like sheet rupture, ligament breakup and their interactions with the jet core and corona flaps are seen. At $t = 3$, the jet is fully developed and we see long axial ligaments with a marked velocity gradient made visible by the colour change along the axis, circular ligaments, corona fingers merging and sheet rupture. A large number of droplets and ligaments are produced in the mushroom tip region. Most of the phenomena illustrated above are also observed with controlled perforation, except the initiation of sheet rupture and annular ligament detachment.

Figure 4 shows a contrasting view of the jet with the camera direction aligned with the jet axis and the view positioned at the inlet. A full evolution of the jet, as seen from the inlet, is available in supplementary movie 6. In figure 4(a) the objects are coloured by curvature, so small droplets are red and small bubbles, with negative curvature, are blue. The bubbles have been trapped in the bulk or core of the jet or inside other liquid masses. The likely mechanism for the formation of a bubble is liquid mass collision, for example, droplets and wavelets of fluid mass impacting on the jet core. Many such
Y. Kulkarni, C. Pairetti, R. Villiers, S. Popinet and S. Zaleski

<figure><img src="image_5.png" /><figcaption>Figure 3. The advancing pulsed jet at various time instants *t* and level *ℓ* = 14. The fluid interface is coloured by the axial velocity and the background is coloured by the vorticity. The background also shows the mesh refinement. (a) The pulsed jet develops a mushroom head and a rim. In (b) the rim detaches. (c) Development of flaps coming from the sinusoidal pulsation. (d) A jet entering in a regime where the effect of pulsation is lost at the mushroom head. (e) A fully developed jet and a rich spectrum of droplets and ligaments. Since *Re*<sub>g</sub> is rather low at 5800, there is relatively little vorticity away from the interface unlike in the case of Kant *et al.* (2023).</figcaption></figure>

droplet or ligament impacts are seen in supplementary movie 1. A surprising aspect of the display is the relatively large number of droplets and bubbles seen. In figure 4(b) the objects are coloured by axial velocity, as with supplementary movies 1 and 3. A further interesting aspect of this view is to show the large dispersion in the axial velocities of the droplets.

At this point in the exposition it may be useful to define the number distribution and the PDF and to recall the connection between the two. We consider a fixed instant of time $t_0$. We partition the interval of droplet sizes ($0, d_M$) into a discrete set of $M$ intervals [$d_i, d_{i+1} = d_i + \delta_i$]. The index $i$ varies from 0 to $m$. The intervals are geometrically distributed, with $d_{i+1} = 1.047d_i$, so $\log(d_i)$ is evenly distributed in log space. The number of droplets with diameter in the interval $[d_i, d_i + \delta_i]$ is $n_{d,i}\delta_i$. The PDF $p(d)$
<figure><img src="image_6.png" /><figcaption>Figure 4. View from the inlet at time $t = 3.04$ showing the inner region of the central core liquid jet. Plot (a) is coloured by the curvature showing the encapsulation of gas bubbles identified by the negative curvature (blue) in the liquid core encircled in the black circle. The droplets have positive curvature (red). The entrained bubbles travel with the core jet velocity and could also result in the formation of a few compound droplets during atomisation or provide a physical breakup mechanism for thin sheets. Plot (b) is the same as (a) but coloured by the axial velocity. The simulation corresponds to level $\ell = 14$ with the MD method applied at level $m = 13$.</figcaption></figure>

is proportional to the expected value $\langle n_{d,i} \rangle$ as

$$p(d_i) = \frac{1}{\mathcal{N}} \langle n_{d,i} \rangle, \qquad (3.1)$$

where the normalising factor is

$$\mathcal{N} = \sum_{i=0}^{M} \langle n_{d,i} \rangle \delta_i. \qquad (3.2)$$

We show two kinds of plots, one with the number $N(d_i)$ of droplets in each bin $i$, that is,

$$N(d_i) = n(d_i) \delta_i, \qquad (3.3)$$

and the other with the PDF $p(d)$. Because of the difference between the definitions of $N$ (which is bin-size dependent) and $n$ and $p$ (which are both bin-size independent), there is a scaling relation

$$N(d) \sim d p(d) \sim dn(d) \qquad (3.4)$$

between the two distributions. (We note that in the literature both $n$ and $N$ are reported with the name ‘number distribution’.) The pulsed jet process is both deterministic (except for the random character of sheet perforations) and transient, so the value of $n_{d,i}$ obtained at $t = t_0$ in a simulation is well defined. The expected value $\langle n_{d,i} \rangle$ differs very little from the numerically obtained value $n_{d,i}$. Moreover, in such a deterministic process no difference is expected between two realisations, that is, two identical simulations should yield identical droplet counts. On the other hand, in a jet with a random upstream injection condition instead of a pulsed injection condition, the simulation is not deterministic and the value of $n_{d,i}$ obtained in a simulation at some time $t$ is fluctuating. The transient nature of the jet
Y. Kulkarni, C. Pairetti, R. Villiers, S. Popinet and S. Zaleski

<figure><img src="image_7.png" /><figcaption>Figure 5. (a) The bubble size distribution for the image shown in figure 4. The number N is defined in (3.3). (b) A 2-D histogram for the same image showing bubble size distribution on the transverse coordinates. The jet is advancing in the x direction, which is the axial direction. The solid white circle is the inlet region. The curly dotted white lines are drawn to highlight pale patches to indicate that some bubbles also exist outside the jet core, indicating possible compound drops.</figcaption></figure>

is also important. In a steady-state simulation (which should be obtained some time after the jet starts flowing out of the simulation domain), initial perturbations are amplified chaotically and the system is again probabilistic. Thus, our pulsating initial condition simplifies somewhat the analysis by yielding a deterministic instead of a probabilistic number density.

Figure 5 shows the number distribution of bubbles. The ordinate $N$ is the number $n_i$ of droplets defined in the text. We see that above the critical thickness, a distribution $N \sim r^{-\gamma}$ is seen. This distribution may have consequences in the sheet rupture and, hence, the droplet-size distribution. We have discussed the physical implications of bubble size distribution in Appendix A. In figure 5(b) we plot a two-dimensional (2-D) distribution of the bubble size on a y-z plane. The white circle represents the inlet. The plot shows that most of the bubbles are near the outer boundary of the jet core where we see wavelet interactions and drop impacts the most (figure 3(c) and supplementary movie 1). The distribution also shows eight characteristic high bubble-density zones. The occurrence of these zones and what decides their number is beyond the scope of the current paper. We also see some bubbles present outside the jet core encircled in dotted white lines. This indicates that we do have a relatively small number of compound droplets.

We now discuss the statistical properties of the jet, in particular, the droplet-size distribution in the case of uncontrolled sheet perforation or ‘no-MD’ case. As the jet advances, the droplet-size distribution evolves. This is shown in figure 6. One can see that the droplet-size distribution has started to converge at $t = 1.8$ and is approximately converged at $t = 2.8$. The histogram thus only shifts in the vertical direction (the number of droplets), while maintaining the same qualitative nature. This histogram is bimodal with two peaks or local maxima. There is a large amplitude peak corresponding to a smaller droplet size (referred to as peak 1) and a smaller amplitude one corresponding to a larger droplet size (referred to as peak 2).

We investigate the statistical convergence of the droplet-size distribution by plotting in figure 7 the histogram at $t = 3.5$ when the jet is fully developed. We plot various grid resolutions or levels on that same figure. The highest resolution in the no-MD case corresponds to the $\Delta_{14}$ line. This implies $D/\Delta = 920$ grid points per initial diameter. Figure 7 has two important characteristics: (i) the histograms at various resolutions are
<figure><img src="image_8.png" /><figcaption>Figure 6. (a) The droplet-size number distribution ($N$ is defined in (3.3)) and (b) the probability density function of the droplet diameter (defined in (3.1)) at various time instants. The plot shows that the PDF has converged in time $t \sim 2.5$. This time convergence determines the choice of the end time of the simulation at $t = 3.5$. Both (a) and (b) have 200 bins. The simulation corresponds to level 14 and the vertical dashed line represents the grid size.</figcaption></figure>

<figure><img src="image_9.png" /><figcaption>Figure 7. (a) The droplet-size distribution and (b) the probability density function of the droplet diameter $d$ at the final time $t = 3.5$ for various grid resolutions. The grid sizes are shown with vertical dashed lines and $\Delta\ell$ is defined as in (2.7). The inset in (a) shows the converged tail region. All the plots are done with 200 bins. A Pareto distribution $N(d) \sim d^{-2}$ is added to compare with the converged region of the distributions.</figcaption></figure>

qualitatively similar, that is, bimodal with a major peak 1 at small scales and a minor peak 2 at large scales; and (ii) despite such high resolution, it is clear that there is no convergence for the droplet-size distribution. However, the tail region of the distribution, to the right of peak 2, shows some degree of convergence. In that tail region a Pareto distribution $N(d) \sim d^{-2}$ is seen reminiscent of such distributions found in other contexts (Balachandar et al. 2020; Pairetti et al. 2021). The dependence of the two peaks with the grid size is shown explicitly in figure 8(b) where we plot the ratios $d_i/\Delta\ell$ corresponding to the peaks as a function of the maximum level $\ell$. We observe that the ratio for peak 1, $d_1/\Delta\ell$, remains constant around a value of 3, while for peak 2, $d_2/\Delta\ell$ weakly increases with $\ell$.
Y. Kulkarni, C. Pairetti, R. Villiers, S. Popinet and S. Zaleski

<figure><img src="image_10.png" /><figcaption>Figure 8. (a) Schematic of the bimodal droplet-size distribution. The arrows indicate the logarithmic distance from the grid size $log d_i/\Delta = log d_i - log \Delta$. (b) The values of $d_i/\Delta$ as a function of the maximum grid refinement level $\ell$ for the distribution of figure 7. The dashed line is at $3\Delta$, implying the droplets at peak 1 have a grid-dependent diameter of $d_1 \sim 3\Delta$.</figcaption></figure>

The grid dependence of these two peaks is synonymous with the absence of statistical convergence. We discuss possible explanations for the behaviour of these two peaks below and, in particular, the relation between peak 1 and the curvature ripples appearing just before sheet rupture. We also discuss a mechanism and a possible scaling for peak 2 in Appendix A. Since the grid dependence behaviour is maintained at all resolutions, the constant scaling of the diameter $d_1$ of peak 1 with the grid size implies that no matter how much we refine the grid, the distribution of droplet sizes is grid dependent.

The mechanism leading to peak 1 is related to the perforation of thin sheets. This perforation is observed to occur in two contrasted ways. A simple mechanism is the impact of small droplets on the liquid sheet. The droplets are formed by breakup previously occurring elsewhere in the simulation. The impacts yield in some cases the formation of an expanding hole, as shown on figure 9. We note that such impacts have already been observed in previous studies such as those of Ménard et al. (2007) and Shinjo & Umemura (2010). Not all droplet impacts create perforations. In figure 9, at time $t = 2.94$, we track the droplets indicated by arrows. We see that droplet impacts often create holes with a characteristic finger-like ligament inside the hole. Impacts are shown by arrows at time $t = 3.10$ in figure 9. In the lower panel at $t = 3.12$ and $t = 3.14$ some droplet impacts do not result in holes but instead droplets merge with the sheet. This hypothesis was confirmed by repeated observations of the sheets at various angles and by visioning supplementary movie 5. The red arrow from $t = 3.14$ and $t = 3.16$ captures an interesting moment, where the droplet impact and sheet rupture happen simultaneously and it is unclear if the hole was created by the numerical rupture process discussed below or by droplet impact.

Another mode of thin sheet perforation, much more complex and purely numerical, is illustrated on figure 10 as well as supplementary movie 5. The curvature colouring of the interface helps us identify the perforation spots. We see high-frequency oscillations of the curvature field in the form of alternating red and blue colours. These indicate that the ripples originate at the fluid interface before the sheet rupture happens. Finally, at $t = 2.82$, we see a dense network of small diameter ligaments, scaling with grid size. Another mechanism that may contribute to peaks 1 and 2 is the rupture of ligaments, also shown on figure 10. Small diameter ligaments have large curvature and are thus easily
<figure><img src="image_11.png" /><figcaption>Figure 9. Droplet impacts on the frontal liquid sheet. The points of interest are indicated by the arrows. Droplet impacts can result in holes with characteristic ligaments as seen at $t = 3.10$. In some cases droplets coalesce into the sheet seen at $t = 3.12$ and $t = 3.14$. At $t = 3.16$, the sheet is ruptured but the droplet is still identifiable. All simulations are for the maximum level $l = 13$. The MD method is not applied. The interface is coloured by axial velocity. The artefacts or corrugated surfaces particularly visible at $t = 3.02$ and $t = 3.10$ at the top right are not aliasing artefacts (due to an approximate isosurface interpolation) but are representative of the curvature oscillations seen in the next figures.</figcaption></figure>

identified by their colour. They break physically and not numerically by the Rayleigh–Plateau instability. The origin of these ligaments themselves is often an expanding hole whose rim collides into other rims.

We note that the high-frequency pressure oscillations are due in part to the use of the height-function method to compute curvature. The height-function method ceases to work when the sheet is too thin. However, other methods will also fail as represented in figure 1. To investigate the small-scale ligament networks, we zoom around such a breakup event in figure 11. We see that a weak spot develops and curvature ripples appear at $t = 2.78$ as soon as the sheet reaches a thickness of a few times the cell size, as shown in the inset. These curvature ripples then give rise to ligament networks inside the expanded hole
Y. Kulkarni, C. Pairetti, R. Villiers, S. Popinet and S. Zaleski

<figure><img src="image_12.png" /><figcaption>Figure 10. Appearance and evolution of the curvature ripples on the interface. The interface is coloured by the interface curvature. These ligaments are coloured a darker red and the interfaces on the thin sheets are closer to white. Ligament rupture is encircled in blue. The rupture or ligament pinch off appears in a darker red colour. Sheet rupture, encircled in green, displays curvature ripples in the form of red-blue oscillations. At time $t = 2.82$, the green-circled rupture region displays a grid-dependent ligament network that has evolved from these ripples. These images correspond to a level $ℓ = 14$ simulation. (See also supplementary movie 5.)</figcaption></figure>

eventually leading to the grid-dependent droplets. To summarise, as the sheet reaches a thickness of order $3Δ$, the rupture mechanism is as follows: thin sheet regions > curvature ripples > ligament networks and tiny droplets ($d \sim 3Δ$). Since this sequence of events is inherent to the numerical breakup of thin sheets, it is impossible to escape its occurrence by increasing the grid refinement, although refinement can delay it. The numerical sheet breakup identified here is the direct cause of the lack of convergence of peak 1. Indirectly, the formation of droplets and ligaments of a given size in a non-converged manner leads to all the droplet size counts being unconverged as these other size droplets are formed through coalescence and breakup events from the peak 1 droplets and similarly scaled ligaments.

We also note that the point of view and region of figure 9 were selected to show a large number of drop impact situations, which are however scarcer than the numerical breakup of sheets.

To give an overall idea of the evolution of the jet, we show the fully evolved jet at $ℓ = 13$ in figure 12 and the same fully evolved jet but at $ℓ = 14$ in figure 13. The mist of tiny droplets near the mushroom head and near the inlet flaps is more pronounced in the $ℓ = 14$ case than in the $ℓ = 13$ case. The $ℓ = 14$ case retains more large drops and deformed ligaments at mid-length compared with the $ℓ = 13$ case. The inset zoom shows the numerical sheet rupture. The ligament network structure is more clearly seen at $ℓ = 14$ than at $ℓ = 13$.

## 3.2. Controlled perforation by the MD method

We now present the results of atomisation when we apply the MD method. Here we adopt the following conventions.

(i) The octree maximum refinement level $ℓ$ defines the finest grid size as given in (2.7).
<figure><img src="image_13.png" /><figcaption>Figure 11. The fully developed jet at $t = 2.8$. The inset shows a zoom in at the sheet rupture spot showing the curvature ripples in the weak spot about to be punctured and the resulting ligament network. This ligament network eventually produces grid-dependent droplets. The simulation shown here is at level $ℓ = 14$.</figcaption></figure>

(ii) The MD level $m$ implies that the critical thickness for punching holes is given by
$$h_c = 3\Delta_m.$$

As an example, the notation $(ℓ, m) = (13, 12)$ would mean that the finest cell size in the simulation is $\Delta_{13}$ and the critical thickness at which holes are being punched is $h_c = 3\Delta_{12} = 3 \times (2\Delta_{13}) = 6\Delta_{13}$. This convention for the definition of the ‘MD level $m$’ maintains consistency with the comments in the code available at http://basilisk.fr/sandbox/lchirco/signature.h.

The extent of the integrated region in the quadratic form computation leads to the correspondence $h_c = 3\Delta_m$. The overall evolution of the jet shows important variations with the no-MD case. Comparing figure 14 MD (14, 13) to figures 12 ($ℓ = 13$) and 13 ($ℓ = 14$), we note that (i) the mist of tiny droplets near the mushroom head has disappeared in figure 14, (ii) the mushroom-head corona-flap sheet is more affected by holes and breakups and is closer to figure 13 $ℓ = 14$ than to figure 12 $ℓ = 13$, (iii) the cloud of
Y. Kulkarni, C. Pairetti, R. Villiers, S. Popinet and S. Zaleski

<figure><img src="image_14.png" /><figcaption>Figure 12. The fully developed jet at $t = 3.03$. The inset shows a zoom in at the sheet rupture spot showing the curvature ripples in the weak spot about to be punctured. The simulation shown here is at level $ℓ = 13$.</figcaption></figure>

droplets at mid-length is similar to both no-MD cases but is more uniform and closer to the figure 12 $ℓ = 13$ case.

Figure 15 shows the evolution of the PDF of the droplet-size distribution with time. Similar to the no-MD plots, we see that here also the statistics are converged in time as soon as $t = 2.5$. When compared with figure 7, we see that the peak corresponding to the smaller droplet size has a smaller amplitude.

We show the number density $n(d)$ in figure 16. Figure 16(a) shows the histograms at various grid resolutions for a fixed MD threshold $h_c = 3\Delta_{12}$. The distribution is again bimodal, but the peak corresponding to the smaller droplets is weaker. We also observe a clear converged region in the interval $[d_c, d_{im+1})$, where $d_c$ is the minimum converged equivalent droplet diameter or ‘converged diameter’ for simplicity. The extent of the converged region increases, that is, $d_c$ decreases as the grid size $\Delta_l$ decreases. In figure 16(a) the vertical line marks $d_c$ for two consecutive grid sizes. We have shaded the onset of the non-converged region between two consecutive grid refinements for clarity.

In practice, the ‘converged diameter’ $d_c$ is determined from the comparison of a pair of distributions from simulations of two different grid sizes $\Delta_l$. We note this pair $\Delta_{l1} - \Delta_{l2}$.
<figure><img src="image_15.png" /><figcaption>Figure 13. The fully developed jet at $t = 3.01$. The inset shows a zoom in at the sheet rupture spot showing the curvature ripples in the weak spot about to be punctured. The simulation shown here is at level $ℓ = 14$.</figcaption></figure>

In figure 16(a) the orange vertical line marks $d_c$ for the pair $\{\Delta_{14} - \Delta_{13}\}$. Similarly, the purple vertical line marks $d_c$ for the pair $\{\Delta_{15} - \Delta_{14}\}$. Some additional information may be obtained by investigating a no-MD case. In such a case, there is a typical sheet thickness at breakup $h_e$ that is comparable to the controlled breakup for $h_c \simeq h_e$. Indeed, the critical breakup film thickness $h_e$ for no-MD is $h_e \simeq 1.5\Delta$, Then the above case $h_c = 3\Delta_{12}$ with MD has a critical sheet thickness equal to, per expression (2.7), $h_e = 1.5\Delta_{11}$ in the no-MD case at level 11. Thus, the no-MD simulation at $ℓ = 11$ is a rough equivalent to the three simulations with MD level $m = 12$. The vertical black line is the converged thickness $d_c$ for the pair $\{\Delta_{13}, \Delta_{11}\}$ where the first simulation is MD and the second is no-MD. One sees on figure 16(a) that the converged diameters, indicated by the solid vertical lines, decrease and tend to approach the critical sheet thickness $h_c$ for MD breakup, represented by the dashed vertical line. We now show the droplet-size distribution and converged region for a smaller critical sheet thickness $h_c = 3\Delta_{13}$. In figure 16(b) we see similar results but with only two instead of three levels, $l = 14$ and 15. Since a simulation at $\Delta_{15}$ is already quite expensive (it implies an equivalent of 35 trillion cells for a uniform grid with $\Delta_{15}$), we did
Y. Kulkarni, C. Pairetti, R. Villiers, S. Popinet and S. Zaleski

<figure><img src="image_16.png" /><figcaption>Figure 14. The (14, 13) jet at $t = 3.1$ with MD applied. The basilisk level $\ell = 14$ and MD level is 13. The insets and the boxes indicate droplets that are near the maximum size of the distribution in the converged region.</figcaption></figure>

not perform a third simulation at $l = 16$ in that case. A zoomed-in view for the converged region $d_c$ for each consecutive pair is shown in Appendix E.

The closely packed peaks around $log(d) = -3$ in figure 16(a), mostly visible for $(\ell, m) = (15, 12)$, are related to a process that produces a narrow spectrum of small droplets of size $\bar{d}$ with $\Delta_{15} < \bar{d} < h_c$. These small droplets coalesce at least two or three times to produce additional small droplets of size $2\bar{d}$, $3\bar{d}$, etc. We are unsure of the mechanism producing these initial droplets of well-defined diameter $\bar{d}$ and why it appears most clearly for the pair (15, 12).

The variation of the converged diameter $d_c$ as a function of $h_c$ and $\Delta_\ell$ is shown explicitly in figure 17. Figure 17(a) illustrates schematically how the converged diameter approaches $h_c$ as $\Delta_\ell$ is reduced. Figure 17(b) shows the data in the $d_c/h_c$ plane. Each data point in that plane corresponds to a simulation pair discussed above such as $\{\Delta_{14} - \Delta_{13}\}$, and to a vertical line in figure 16. Following a vertical line in figure 17(b), one can see $d_c$ tending to $h_c$ as the grid is refined for a fixed $h_c$.
<figure><img src="image_17.png" /><figcaption>Figure 15. The probability density function at various times t for the droplet-size distribution when MD is applied. The PDF has converged in time at $t = 2.5$. (a) The PDF at maximum level $ℓ = 13$ and MD level $m = 12$ and (b) the PDF at maximum level $ℓ = 14$ and MD level $m = 12$. The critical hole punching thickness $h_c = 3Δ_{m=12}$ is shown as the black dashed line and is the same for both plots.</figcaption></figure>

Figure 18 shows a zoom in at a hole expansion that was initiated by the MD method. The phase modification region is a cube but in less than 5–10 time steps the hole takes a rounded shape due to surface tension and starts expanding. When these holes expand, we do not see the noisy curvature ripples that resulted in grid-dependent droplets in the no-MD simulation of figure 11. This qualitative improvement in sheet collapse helps to obtain grid convergence.

The MD procedure we use is tuned to a relatively large frequency of hole formation. This results in features such as those shown in figure 19. We see that once we detect a thin sheet and the MD process generates holes, these expand. However, it often happens that before the holes have expanded over the entire sheet more holes are punched just outside the rim. These kind of holes are identified as black circles in figure 19. This effect is also seen in supplementary movie 4 that shows the evolution of thin sheets with MD. Experimental photographs and the reasoning in Appendix A point to a smaller number of holes. This smaller number could be achieved by an improved MD algorithm that would involve a more realistic model of the physical process of hole formation.

The final result of our manifold-death-enhanced study of droplet sizes is summarised on figure 20. In this figure we select two distributions of droplet sizes at two MD levels, $m = 12$ and $m = 13$, and we keep the ratio of the grid size to the MD threshold $h_c/Δℓ$ constant so the level difference is always $ℓ - m = 2$. The two obtained distributions are similar to each other, with a sudden increase in droplet number at the edge of the converged region. Moreover, these plots show that a sufficiently thin grid produces a bump-shaped distribution with a clear maximum (this is the case at $m = 13$ but not at $m = 12$). This maximum is located at $d_{max}/D = 4.3 10^{-2}$ in our case. It is interesting to try to understand the origin of this bump. By dimensional analysis the maximum should be expressed as

$$d_{max} = Df(h_c/D, N_1, \dots, N_n), \qquad (3.5)$$

where the $N_i$ are the dimensionless numbers of the problem other than $h_c/D$. It is quite tempting to assume that in some range of $h_c/D$ and of the dimensionless numbers a proportionality relation holds. From the existing data it would be

$$d_{max} = 4.3 h_c. \qquad (3.6)$$
Y. Kulkarni, C. Pairetti, R. Villiers, S. Popinet and S. Zaleski

<figure><img src="image_18.png" /><figcaption>Figure 16. The droplet-size distribution from simulations using the MD method. Each curve is labelled $\Delta\ell - \Delta_{MD=m}$, where $\ell$ is the grid level and $m$ is the MD level. The dashed black vertical line represents $h_c = 3\Delta_{m}$, that is, the critical thickness of punching holes. The solid vertical lines indicate the values of the converged diameter $d_c$, corresponding to the smallest diameter above which the difference between subsequent distributions is not significant. All curves in (a) have fixed $h_c = 3\Delta_{12}$ and all curves in (b) have fixed $h_c = 3\Delta_{13}$. Note that the $\Delta_{11}$ curve in (a) and the $\Delta_{12}$ curve in (b) indicate the no-MD method and roughly correspond to respectively $h_c = 3\Delta_{MD=12}$ and $h_c = 3\Delta_{MD=13}$, given the breakup due to a lack of sheet resolution in VOF simulations. The shading highlights the departure from the converged region. We see that $d_c$ approaches $h_c$ as the grid size is reduced. The inset shows the important change in the typical sizes of droplets in the first peak above $h_c$ as MD is applied. Peak A corresponds to peak 2 in the no-MD distributions of figure 7 while peak B, with much smaller droplets, corresponds to droplets seen in the MD figure 14.</figcaption></figure>
<figure><img src="image_19.png" /></figure>

<figure><img src="image_20.png" /><figcaption>Figure 17. (a) Schematic of the comparison of two simulations and how it allows us to define the ‘converged diameter’ $d_c$. It is seen that $d_c$ approaches $h_c$ as the simulations are refined. (b) Plot of $d_c/h_c$ versus $h_c$ corresponding to the $d_c$ values of figure 16. Along a fixed $h_c$, the $d_c/h_c$ points move closer to unity as the maximum level is increased. The black squares corresponds to the no-MD case, while all circles are for the cases with MD.</figcaption></figure>

<figure><img src="image_21.png" /><figcaption>Figure 18. Sheet perforation when the MD method is applied. Unlike figure 11, the holes are not preceded by large regions containing curvature ripples. Holes are punched at $t = 3.01$ and expand with time as shown. The interface in the left part is coloured by curvature and on the right is coloured by the axial velocity. The simulation shown here corresponds to $(\ell, m) = (14, 13)$.</figcaption></figure>

We note that this is to our knowledge the only example of a *converged* bump-shaped (unimodal) distribution in simulations of jet atomisation at a large Weber number. We attempt to fit this bump shape to standard distributions in the literature, namely the gamma and the log-normal. Both distributions are acceptable in the large *d* range, above 0.005, as
Y. Kulkarni, C. Pairetti, R. Villiers, S. Popinet and S. Zaleski

<figure><img src="image_22.png" /><figcaption>Figure 19. Further evolution of the holes of figure 18. The image is coloured by axial velocity. We see that near the rims of the main expanding holes, additional holes are punched, encircled in black. (Manifold death is applied.)</figcaption></figure>

<figure><img src="image_23.png" /><figcaption>Figure 20. Two distributions of droplet sizes at two MD levels, $m = 12$ and $m = 13$, with the same ratio of the grid size to the MD threshold $h_c/\Delta_l$. We only plot the converged region, that is, $d > d_c(h_c)$. The most refined simulation ($l, m = (15, 13)$) yields a distribution somewhat but not exactly similar to the less refined. Both distributions show a sudden increase in droplet number (a 'ridge') as the droplet size is decreased at the edge of the converged region. The vertical line shows the limit of convergence $d_c$. Plot (a) shows the un-rescaled number distribution, while plot (b) shows the PDFs. Plot (a) is clipped to show the PDFs $p(d)$ for $d > d_c(m = 12)$, where $d_c$ is the edge of the converged region for $m = 12$. The Pareto distribution $n(d) \sim d^{-2}$ is also plotted.</figcaption></figure>

seen in figure 21. The fitting distributions are defined as follows:

$$ \text{Log-normal: } P(x; \mu, \sigma) = \frac{1}{x\sigma\sqrt{2\pi}} \exp\left[-\frac{1}{2}\left(\frac{\ln(x) - \mu}{\sigma}\right)^2\right], \quad (3.7) $$

$$ \text{Gamma: } P(x; \alpha, \beta) = \frac{\beta^{\alpha}}{\Gamma(\alpha)} x^{\alpha-1} \exp(-\beta x). $$

## 4. Conclusions

We have studied an atomising, pulsed round jet at a succession of grid sizes with unprecedented resolution. Visual inspection of the resulting interface topology reveals a new numerical phenomenon preceding hole formation in the VOF method: high-frequency curvature oscillations followed by the appearance of a small-scale network of ligaments. While the experimental observation and theoretical analysis of droplet formation by expansion of holes is not new, this is the first time that these high-frequency oscillations
<figure><img src="image_24.png" /><figcaption>Figure 21. Same as figure 20 but with a linear abcissa and a fit to gamma and log-normal distributions.</figcaption></figure>

of curvature are observed in VOF simulations. The inspection of droplet-size frequency distributions reveals the presence of two large peaks tied to the grid size. It is clear that the small-$d$ peak and probably the second one are related to the unphysical breakup.

To mitigate this effect, we apply the MD procedure of Chirco et al. (2022). The procedure forces hole formation when sheets reach a set thickness $h_c$. Number frequency plots show convergence in a range of droplet sizes starting at a critical diameter $d_c$. The extent of this converged range increases when either the grid size $\Delta_l$ or the critical thickness $h_c$ are decreased. For a fixed $h_c$, when we reduce $\Delta_l$, the converged diameter $d_c$ approaches the critical thickness $h_c$. We thus recover the convergence properties observed by Chirco et al. (2022) and by Tang et al. (2023). We characterise them in the $h_c-d_c$ plot of figure 17(b). We note that the statistical accuracy is much stronger in the present case since our pulsed jet produces three orders of magnitude more droplets than the phase inversion case of Chirco et al. (2022). We do not have a definite explanation for the lack of convergence below the diameter $h_c$, that is, we are not sure why $d_c$ does not decrease below $h_c$ as $\Delta_l$ is decreased. We can however offer several avenues for improvement of this feature and for future exploration. One is the smoothness of the MD procedure itself and the choice of its parameters such as the number of holes punched or the punching frequency $1/\tau_m$. For the first time in VOF simulations of atomisation, a bump-shaped, converged droplet-size distribution is obtained. Even without convergence, bump-shaped distributions are harder to obtain in VOF simulations than in level-set or diffuse-interface simulations where they are common (Herrmann 2010; Khanwale et al. 2022; Saurabh et al. 2023), because these two other methods tend to ‘evaporate’ small droplets, as shown on figure 1. A scaling relation for the position of the bump is suggested. The bump only appears at the highest level. Whether the proposed scaling (3.6) for the bump would persist if yet more refined simulations are performed is an interesting open question.

Beyond the physics of hole formation and expansion, the simulations reveal a rich spectrum of elementary phenomena that are the building blocks of atomisation. These include networks of well-resolved ligaments in and around holes, ligaments detached from the main liquid core and either aligned or at an angle with it, bubbles trapped in liquid
Y. Kulkarni, C. Pairetti, R. Villiers, S. Popinet and S. Zaleski

regions, droplets impacting on sheets, etc. This rich physics is seen in increasing detail as grid resolution is increased. For example, the sharpness and apparent realism of the various expanding holes, with attached ligaments or networks of ligaments inside, was not seen in previous simulations, for example, in Ling et al. (2017). Moreover, a large number of each type of object or ‘building block’ are visible at each instant of time, leading to more significant and converged statistics.

A tantalising question is that of the physical, as opposed to numerical, origin of the perforations. We attempt some discussion of this issue in Appendix A. Two types of rate processes for heterogeneous nucleation are postulated, one based on in-sheet germs or nuclei, such as small bubbles, the other on external perturbations such as floating droplets or particles. Our simulations do not allow us to decide at present between these two processes for hole formation. However, they may in the future help us to understand how to discriminate between them.

Unlike most of the other round-jet simulations in the literature, ours is bimodal and even trimodal after the application of MD. It seems however that a single mode emerges in the converged range of the best resolved MD simulations, indicating that the appearance of a maximum of droplet size well in the converged region is possible only if the physics is correctly captured down to scales much smaller than the integral scales such as the jet diameter *D*.

Another important characteristic of simulations and experiments are the state of the fluids at the inlet. The primary mixing-layer instability in our simulations is initiated not by upstream turbulence, as in most round-jet simulations in the literature, but by a pulsating flow. In future investigations it would be interesting to compare the visual aspects of the fragmenting interface as well as number distributions and PDFs, between simulations with or without noisy/turbulent inflow conditions.

The other perspectives opened by this work are numerous. Clearly, simulations of comparable quality at a higher Weber number may allow investigation of the ‘fibre-type’ regime observed in many experiments. Currently, our simulations sit on the boundary between membrane type and fibre type. Simulations of comparable resolution should also be performed in the planar mixing-layer case investigated experimentally by the Grenoble group (Ben Rayana et al. 2006; Fuster et al. 2013) and numerically by Jiang & Ling (2021). The latter authors performed very high resolution simulations with hole formation and interestingly obtained a bimodal distribution, but they did not apply a controlled perforation scheme as in this paper. Simulations of secondary atomisation, that is, the atomisation of a droplet suddenly plunged in a high-speed flow, has already been performed together with controlled perforation by Tang et al. (2023) achieving convergence of the PDF. The secondary atomisation case is probably the most promising for an eventual comparison of the distribution of droplet sizes obtained by simulation with the rich experimental data reviewed and modelled, for example, by Jackiw & Ashgriz (2022).

More distant perspectives or more difficult problems remaining to solve are in our opinion twofold. One is clearly the improvement of the MD procedure, to have a better control of the number of holes punched in each thin sheet region, and a smoother manner to ‘cut out’ the initial hole. The breakup frequency and number of attempted breakups were chosen in a rather arbitrary way. These values are too high as shown above. An adequate breakup frequency would correspond to smaller values of *f*<sub>b</sub> and *N*<sub>b</sub>. A caveat is however that if the breakup frequency (*f*<sub>b</sub>, *N*<sub>b</sub>) is too low, some bags may escape MD and breakup numerically. An expensive trial-and-error procedure could be attempted to tune the breakup frequency. Another avenue of improvement would be to modify the MD method to perform an analysis of the bag region to identify the connected thin
sheet regions. Then a controlled number of perforations could be performed in any of these regions. This controlled number of perforations per thin region could be inspired by experimental observations or by a theory such as the one in Appendix A. A final, intriguing possibility is that if a level-set or diffuse-interface method is used in conjunction with MD, and the frequency is on the low side, numerical breakup will be avoided in the regions that would not be sufficiently perforated because of the low frequency. There would still be mass loss because of the evaporation of thin sheets, but no curvature oscillations nor an excessive number of small droplets scaling with grid size.

Another challenging future task is to use this type of simulation to infer, by careful matching of simulation results and experiments, the physical mechanisms for weak spots and sheet perforations. Given the wide variety of such mechanisms that have been advanced, this is a formidable enterprise for which we hope this paper could provide partial guidance.

**Supplementary movies.** Supplementary movies are available at https://doi.org/10.1017/jfm.2025.218.

**Acknowledgements.** We thank Marco Crialesi and Daniel Fuster for fruitful discussions. S.Z. would like to thank Ousmane Kodio for mentioning the Kibble-Zurek theory in discussions on pattern formation.

**Funding.** This project has received funding from the European Research Council (ERC) under the European Union's Horizon 2020 research and innovation programme (grant agreement nr 883849). We thank the European PRACE group, the Swiss supercomputing agency CSCS, the French national GENCI supercomputing agency and the relevant supercomputer centers for their grants of CPU time on massively parallel machines, and their teams for assistance and the use of Irene-Rome at TGCC and ADASTRA at CINES. This project was provided in particular with HPC computing and storage resources by GENCI at TGCC thanks to the grants 2023-A0152B14629 and 2022 SS012B15405 on the supercomputer Joliot Curie's SKL partition.

**Declaration of interests.** The authors report no conflict of interest.

## Appendix A. A phenomenological theory for hole nucleation and expansion

In what follows we describe hole nucleation and expansion as a random process. Its consequences on ligament size are given based on the hole nucleation rate in a manner independent on the exact physical mechanism of nucleation. An analogy between this description and the Kibble-Zurek theory is made.

In atomisation processes, sheet rupture may have several causes. Without deciding between these causes, we distinguish two mathematical forms of the hole forming process (Villermaux 2020). In one it is in-sheet germs or nuclei, such as small bubbles or small droplets of an immiscible liquid (Vernay et al. 2017) that create the holes. In the other one, external nuclei such as floating droplets, particles or hot puffs are the cause of perforation. The probability of punching a hole per unit area of thin sheet can be written to model both processes as

$$d\mathcal{P} = \mathcal{A}(h)dt + \mathcal{B}(h)dh. \qquad (A.1)$$

Here $\mathcal{A}(h)dt$ is the probability that an external object that would be able to punch a hole in a sheet of thickness $h$ actually hits the sheet during the time interval $dt$, while $\mathcal{B}(h)dh$ is the probability that an in-sheet heterogeneity is present that will nucleate a hole while the sheet thins from $h$ to $h-dh$. In our present study we do not have external objects such as solid particles, but we still have holes appearing due to the drop impact. Thus, in our simulations, the atomisation driven by drop impact is mostly contained in $\mathcal{A}(h)$, while those appearing due to numerical perforation of sheets are contained in $\mathcal{B}(h)$. Our simulation captures an overall picture and it is difficult to isolate the effect of any one of the two. In the MD method, we focus on tuning the hole appearance due to the in-sheet
Y. Kulkarni, C. Pairetti, R. Villiers, S. Popinet and S. Zaleski

heterogeneity part having a simplistic model for $\mathcal{B}(h)$ as a step function with critical size being $h_c$. However, in fluids, heterogeneities such as bubbles do not have a single size but a distribution of sizes, and correlations obtained experimentally may be obtained that yield $n_b(r)$ (Brujan 2010) such that $n_b(r) dr$ is the number of bubbles with size between $r$ and $r+dr$ per unit volume. The number of bubbles between $r$ and $2r$ then scales like $n_b(r)r$ and the average distance between such bubbles is $l_1 = [rn_b(r)]^{-1/3}$. A natural requirement is that this distance between bubbles is much larger than their radius. Then one must have $[rn_b(r)]^{-1/3} \gg r$. Since we probe the distribution at length scales $r \sim h$ that are relatively small, it makes sense to require that $n_b(r) \ll r^{-4}$ in the limit of small $r$. A simple model for the distribution of nuclei or number of nuclei per unit volume is

$$n_b(r) \sim c_n r^{-\gamma}, \qquad (A.2)$$

where $\gamma$ is an exponent to be determined by measurement and $c_n$ is a constant. By the requirement above the exponent $\gamma < 4$ and indeed Gavrilov (1969) found $\gamma \simeq 3.5$ while measurements by Shima & Sakai (1987) give $\gamma$ between 2 and 4. As a result, the volume fraction $r^4 n_b(r)$ of bubbles between $r$ and $2r$ vanishes at small $r$. Now suppose that bubbles of radius $r$ are the cause of weak spots in a sheet of thickness $r$. Then from (A.1) the probability of ‘in-sheet’ nucleation in a region of extent $\lambda^2$ while the sheet thickness decreases from $h$ to $h/2$ scales like $\mathcal{B}(h)h\lambda^2$. This nucleation is caused by small bubbles in the volume $\lambda^2 h$ then $\mathcal{B}(h)h\lambda^2 = n_b(h)\lambda^2 h^2$, and the ‘in-sheet’ nucleation probability per unit area is just

$$\mathcal{B}(h) = n_b(h)h. \qquad (A.3)$$

The typical distance $\lambda$ between weak spots in a sheet of size $h$ thinning to $h/2$ is given by the condition that the expected number of weak spots in area $\lambda^2$ is one, which reads

$$n_b(h)\lambda^2 h^2 \sim 1. \qquad (A.4)$$

Then

$$\frac{\lambda}{h} \sim (n_b(h)h^4)^{-1/2} \sim h^{-2+\gamma/2}, \qquad (A.5)$$

so for $2 < \gamma < 4$, one has $\lambda \gg h$, while $\lambda$ still tends to vanish at small $h$. The distance between holes $\lambda$ is much larger than the thickness of the sheet but is still small compared with the total size of the sheet. We shall call $L$ the total size of the sheet so the above reads $L \gg \lambda$. In other words, the holes are rather numerous in each sheet.

This is not what is observed in experiments, for example, by Kant et al. (2023), and one instead seems to see one or two holes per sheet. To resolve this contradiction, one must take into account the speed of expansion of the holes that form during sheet thinning. Indeed, expansion of the holes brings the process of hole formation to a halt after a short time $\tau_1$ that we determine below. We can now give a physical meaning to the critical thickness $h_c$ in the MD algorithm. The expected value of the number of nucleations in a thin sheet of area $L^2$ while the sheet retracts from $h_c$ to $h_c/2$ is

$$N_0(h_c) = L^2 \int_{h_c/2}^{h_c} \mathcal{B}(h)dh \sim h_c^{2-\gamma}. \qquad (A.6)$$

The critical thickness $h_c$ is then that for which

$$N_0(h_c) \sim 1 \qquad (A.7)$$

is of order unity. Independently of the above prediction for $\mathcal{B}(h)$, and only on the basis of a model where the nucleation probability increases as $h(t)$ decreases, it is possible to predict
the effective number of holes that will be observed. First consider the typical size $L(t)$ of an expanding bag. Experimental observations (Opfer et al. 2014; Kant et al. 2023) show that bags inflate exponentially. We consider only large-scale bags governed by the large-scale properties of the flow that thus inflate at a rate $\omega = 1/(2\tau_c)$. Then $L(t) \sim L_0 \exp(\omega t)$ and the bag thickness decreases as

$$h(t) \sim h_0 e^{-t/\tau_c}. \qquad (A.8)$$

Theory and experimental observations (Villermaux & Bossa 2009; Opfer et al. 2014; Marcotte & Zaleski 2019; Kant et al. 2023) indicate the scaling

$$\tau_c \sim \frac{D}{U} \left( \frac{\rho_l}{\rho_g} \right)^{1/2}. \qquad (A.9)$$

Once weak spots are activated in such a bag they expand at the Taylor–Culick velocity (Taylor 1959; Culick 1960)

$$V_T = \left( \frac{2\sigma}{\rho_l h_c} \right)^{1/2}. \qquad (A.10)$$

If $N$ holes are activated in a typical bag of size $L$, the average distance between hole nuclei positions scales like $\lambda = LN^{-1/2}$. This distance is covered by the expanding hole edges in a time $\tau_1 = \lambda/V_T$ and

$$\tau_1(N, h) = LN^{-1/2} \left( \frac{2\sigma}{\rho_l h} \right)^{-1/2}. \qquad (A.11)$$

We now distinguish between two extreme cases. In the first case $\tau_1(1, h_c) \ll \tau_c$. This means that after the first hole has formed, it will expand to the whole bag size before the sheet has significantly thinned. During the expansion of that first hole, $h$ will vary from the initial thickness noted $h_0$ to $h_0e^{-\tau_1/\tau_c} \simeq h_0 - h_0\tau_1/\tau_c$ and the expected number of other holes that could form during that time is

$$N_1(h_c) = L^2 \int_{h_0-h_0\tau_1/\tau_c}^{h_0} \mathcal{B}(h)dh \simeq N_0(h_c)\tau_1/\tau_c, \qquad (A.12)$$

and from (A.7) we obtain $N_1(h_c) \ll 1$, that is, the expected number of additional holes is much less than one, which means that in practice typically a single hole will be observed. Another way of deriving this result is to say that the hole expansion time is much shorter than the hole nucleation time. In the opposite extreme, the hole expansion time for a single hole is much longer than the hole nucleation time. One then needs to consider the expansion time $\tau_1(N)$ for a large number of holes $N$. At this point it is useful to borrow an argument from the theory of bifurcations or instabilities in time-varying environments. The Kibble–Zurek theory (Kibble 1976; Zurek 1985) argues that the time scale of the environment variation is equal to the time scale of the growth of the instability. We thus argue that the multi-hole expansion time $\tau_1(N)$ scales like the characteristic bag expansion time $\tau_c$, that is,

$$\tau_1(N) = \tau_c. \qquad (A.13)$$

To obtain this rule, consider how nucleation occurs as time progresses from the first perforation at time $t_0$ to a time where most of the surface is covered by expanded holes. At first the nucleation rate is near the threshold given by (A.7). As time progresses, previously nucleated holes expand and the thickness decreases. The decrease in thickness implies a higher nucleation rate since $\mathcal{B}(h)$ increases as $h$ decreases. Hole density increases and the
Y. Kulkarni, C. Pairetti, R. Villiers, S. Popinet and S. Zaleski

typical distance $\lambda(t)$ between holes decreases. When $\lambda(t)$ becomes of the same magnitude as the distance covered by the expanding rims, one has $\lambda(t) \sim V_T(t - t_0)$ and the process stops with $t - t_0 \sim \tau_1(N)$. Since we are in the second extreme case $N \gg 1$, the hole density and nucleation rate have increased considerably from that at $t_0$, which means that the sheet thickness has decreased significantly. A significant decrease of the sheet thickness requires a time of order $t - t_0 \sim \tau_c$, hence, (A.13). Then from the balance of time scales (A.13) the number of holes in a bag of size $L$ can be obtained as

$$N_1(h_c) \sim \frac{h_c \rho_g L^2 U^2}{D^2 \sigma}. \qquad (A.14)$$

There is a degree of uncertainty for the ratio $L/D$, however, $L$ and $D$ are both integral scales of the turbulent flow and can be considered of the same order of magnitude. (We disregard the possible case where eddies of smaller scale than the integral scale of turbulence are sufficiently energetic to generate smaller-scale bags.) The number of holes per bag in this mechanism is then the ratio of two typically large numbers, $D/h_c$ and $We_g$. Equation (A.14) is valid in the second extreme case when the right-hand side is much larger than unity. A very rough extrapolation between the two extreme cases is then

$$N(h_c) \sim \text{Max} [1, N_1(h_c)]. \qquad (A.15)$$

In most cases one expects $L$ to be of the same order of magnitude as $D$, which leads to the simplified prediction

$$N_1(h_c) = We_{h_c}, \qquad (A.16)$$

where $We_l = \rho_g l U^2 / \sigma$. Table 2 gives some typical values for $We_{h_c}$ and $We_\Delta$, allowing a prediction of the expected number of holes in the corresponding simulations. We note that the estimate above is only based on the assumption that $B(h)$ increases as $h$ decreases and does not depend on a specific model such as the distribution of bubbles (A.2).

A more sophisticated model of the number of holes can be obtained as follows. Considering the time $\tau_2(N)$ taken by the rate of nucleation $B(h)$ to increase from the initial rate, producing a single hole $N_0 \sim 1$, to the final nucleation rate, producing a large density of holes, the final number of holes $N(h_f)$ is given by

$$N(h_f) = L^2 \int_{h_f}^{h_c} B(h) dh. \qquad (A.17)$$

Then

$$N(h_f) \sim L^2 [(h/h_c)^{2-\gamma}]_{h_f}^{h_c} \qquad (A.18)$$

and

$$N(h_f) \sim L^2 (h_f / h_c)^{2-\gamma}. \qquad (A.19)$$

On the other hand, using (A.8),

$$h_f = h_c e^{-\tau_2(N)/\tau_c}, \qquad (A.20)$$

where, for simplicity, we write $N$ for $N(h_f)$. Hence,

$$N \sim L^2 e^{-\tau_2(N)(2-\gamma)/\tau_c}, \qquad (A.21)$$

and thus,

$$\tau_2(N) = \frac{\tau_c \ln N}{\gamma - 2}. \qquad (A.22)$$
<figure><img src="image_25.png" /><figcaption>Figure 22. A schematic view of how expanding holes in a thin sheet merge into ligaments.</figcaption></figure>

Here we use again an equality of time scales assumption so that $\tau_1 \sim \tau_2$ and (A.13) becomes

$$ \tau_1(N) = \frac{\tau_c \ln N}{\gamma - 2}. \qquad (A.23) $$

This brings a $\gamma$-dependent and logarithmic correction to the prediction (A.13) and, hence, to (A.14).

Beyond the number of holes, another interesting prediction is that of the final size of ligaments obtained from hole expansion. There are at least two possible regimes for hole expansion, either unsteady fragmentation or smooth expansion. In the unsteady fragmentation regime the advancing rims bordering the expanding hole continuously shed small droplets (Wang & Bourouiba 2018). In the smooth expansion regime, the rim advances without shedding droplets and volume is conserved. At the end of expansion a network of ligaments of diameter $\delta$ is formed, as schematised in figure 22 and shown numerically by Agbaglah (2021). The number of ligaments is of the same order as the number of holes $N$ and their length is of order of the diameter of the hole at the end of expansion $\lambda$. The volume of ligaments is equal to the initial sheet volume so that equating the initial and final volumes

$$ N\pi\lambda^2h_c/4 \sim N\pi^2\lambda\delta^2/4, \qquad (A.24) $$

from which we can deduce the ligament diameter

$$ \delta \sim \sqrt{h_c\lambda}. \qquad (A.25) $$

On the other hand, the number of holes is related to $\lambda$ by

$$ N_1 = \frac{L^2}{\lambda^2}. \qquad (A.26) $$

Then using (A.14) and (A.26) we obtain

$$ \delta \sim h_c^{3/4}L^{1/4}. \qquad (A.27) $$

It is possible, although this is a purely heuristic guess, that the second peak observed in the simulations with uncontrolled sheet perforation is caused by the breakup of this type of ligament into droplets. To obtain a simple scaling, we consider $L \simeq D$ and $h_c \simeq \Delta$ when sheet perforation is uncontrolled. The Rayleigh–Plateau instability of the ligament will yield $d_2 \sim \delta$ for the second peak. One then obtains

$$ \frac{d_2}{\Delta} \sim \Delta^{-1/4}. \qquad (A.28) $$

When plotted alongside our numerical data (figure 8b), a fair agreement is obtained.
Y. Kulkarni, C. Pairetti, R. Villiers, S. Popinet and S. Zaleski

<figure><img src="image_26.png" /><figcaption>Figure 23. Thin structure detection and MD perforation in an example 2-D case. The cells tagged as a thin structure at the MD level are filled with a red colour. The upper image corresponds to the fluid structure right before perforation when a thin sheet is detected and the lower part shows the fluid structure after perforation. Note that the images are at the same time stamp to explicitly show the perforation process. The image is done for case $(\ell, m) = (9, 8)$. Adaptive mesh refinement can be seen in the background.</figcaption></figure>

## Appendix B. Grid adaptation

As discussed in the main text, the mesh is adapted dynamically, by splitting the parent cells whenever the local discretisation error estimated by a wavelet approximation exceeds a threshold. This threshold is set at $\epsilon_c = 0.01$ for the error on the volume fraction variable $c$ and to $\epsilon_u = 0.1$ for the error on the velocity field. To investigate the influence of the relatively large error level on the velocity, we plot the histograms of droplet sizes at various error thresholds in figure 24. In that figure it is seen that although the total number of droplets changes significantly when $\epsilon_u$ decreases from 0.1 to 0.01, the position of the relative maxima and, hence, the typical droplet sizes vary relatively little. It is interesting to illustrate how grid adaptation interacts with ligament detection and perforation in the signature and MD methods. The signature method (Chirco et al. 2022) typifies the phase distribution with an index $i_s$. This index $i_s$ is determined by the signature of the quadratic form of Chirco et al. (2022). In summary, $i_s$ can take four values. The index $i_s = -1$
<figure><img src="image_27.png" /><figcaption>Figure 24. The droplet-size frequency for maximum level $l = 13$ and no MD at various error thresholds for the velocity.</figcaption></figure>

designates a gas phase, the index $i_s = 2$ designates a liquid phase while the index $i_s = 0$ designates an interface. The fourth value that the index takes is $i_s = 1$, which designates a thin film region. Note that this detection is performed at the MD level and is later prolonged to the maximum level. Moreover, the thin structure detection procedure acts without affecting the existing mesh refinement. Hence, the mesh refinement is controlled only by the thresholds on the volume fraction and velocity field and in no way depends on the thinsheet detection procedure. We illustrate the process of controlled sheet perforation using the MD method in figure 23. We see that we have detected a thin region (red) in the upper image and a part of this region is eventually punctured in the lower image.

## Appendix C. Estimates of the mass loss due to the MD

In the MD method we artificially punch holes in thin structures. This results in some fluid disappearance and, hence, mass loss. Here we give an estimate of the mass loss for the case $(l, m) = (14, 13)$, that is, the simulation corresponding to the red line in figure 16(b). Figure 25(a) shows how the total mass behaves as a function of time. It is a linear rise in time with a tiny sinusoidal pulsation caused by the injection condition. Since we punch holes in these structures every $t+=0.01$, a tiny decrease from this linear rise is expected every $+=0.01$ interval. The inset zoomed at $t = 1.35$ shows this decrease in mass. The decrease is so small that it is not visible on the scale of the outer plot. To quantify this decrease, we calculate two quantities illustrated in the schematic of figure 25(b). The $\Delta M$ is the difference in the total mass at the time stamp just before punching a hole and the time stamp just after punching the hole. The plot as a function of time is shown in figure 25(c), where the $\Delta M$ values are indicated every $t+=0.01$. Note that as mass is constantly being injected, this implies that the number of thin structures will increase with time and, hence, the $\Delta M$ is seen increasing with an oscillating behaviour. To get a non-dimensional estimate, we calculate the DIFF defined as

$$ \mathrm{DIFF} = 2M(i) - M(i-1) - M(i+1) \qquad (\mathrm{C.1}) $$

(see also the sketch in figure 25b) and calculate the percentage mass loss as

$$ \text{Mass loss (\%)} = \frac{\text{DIFF}}{M[i]} \times 100. \qquad (\text{C.2}) $$
Y. Kulkarni, C. Pairetti, R. Villiers, S. Popinet and S. Zaleski

<figure><img src="image_28.png" /></figure>

<figure><img src="image_29.png" /></figure>

<figure><img src="image_30.png" /></figure>

<figure><img src="image_31.png" /><figcaption>Figure 25. Estimates of the mass loss with MD performed at levels $(\ell, m) = (14, 13)$. (a) Total mass as a function of time. Inset shows the tiny droplet at time $t = 1.35$ due to artificial punching of thin structures every $t+=0.01$. (b) Schematic of the data points and the procedure to estimate the mass loss. The value at $i$ represents the mass just before punching the holes and the value at $i + 1$ represents the value just after punching the holes. Here $\Delta M$ and DIFF are used to do mass loss estimates. Panel (c) shows $\Delta M$ at every $t+=0.01$. (d) Mass loss in percentage calculated as $(\text{DIFF}/M[i]) \times 100$ and the red line is the best fit line done for $t \in (0.8, 1.4)$.</figcaption></figure>

The rationale of using DIFF over $\Delta M$ is that since we are injecting some liquid constantly, we want to estimate how the mass differs not only from the value just before punching holes but from the expected value that it was supposed to have. The percentage loss is shown in figure 22(d), where we can see a decrease in time. We fit an exponential power law of mass loss percentage as $a \exp(-bt)$ and in the steady-state time interval, we see that $b \sim \mathcal{O}(1)$. The mass loss relative to the total mass is going to zero as the jet advances. The question of the effect of the mass loss on the PDF also arises. Since the mass loss occurs in very thin sheets, it affects the very small-scale droplets produced in these sheets and has a negligible effect on the part of the PDF above the converged radius $d_c$.

## Appendix D. Evolving initial mushroom: MD in action

Figure 26 shows the evolution of the initial mushroom jet and the first sheet rupture. Figure 26(a) shows the no-MD case at $\ell = 14$. We can see in the inset at $t = 0.15$, numerical sheet rupture happens in the flap connecting the rolled-up mushroom rim with the mushroom head. The inset shows how a spectrum of tiny droplets forms along with a circular ligament rim at $t = 0.17$. Figure 26(b) shows the first sheet rupture for the case with MD method applied, corresponding to $(\ell, m) = (15, 13)$, and the one corresponding
<figure><img src="image_32.png" /><figcaption>Figure 26. Evolution of the initial mushroom and the first breakup in (a) no MD and (b) with MD. The interface is coloured by the curvature field in both cases. Image (a) shows a detailed view of the first numerical sheet rupture for an $ℓ = 14$ simulation. Image (b) corresponds to the simulation with the MD method applied corresponding to $(ℓ, m) = (15, 13)$. The images shown in (b) are at time $t = 0.1, 0.11, 0.12, \dots, 0.17$. The detachment of a circular ligament rim is seen in both cases at $t = 0.17$.</figcaption></figure>

to the blue line in figure 16(b). One can see the sheet rupture at $t = 0.11$ is much more clear now and the curvature oscillations are not present. In this case too we see detachment of the rolled-up tip in forming a circular ligament at $t = 0.17$. Simple analysis suggests that sheet thickness near the tip at the moment of rim formation should scale as $D/We_g$ (see, for example, Marcotte & Zaleski 2019). The condition $D/We_g \gg \Delta$ is then necessary for proper resolution of the sheet. This condition is equivalent to $We_\Delta \gg 1$, which is verified only for $ℓ = 15$ and perhaps 14 according to table 2. At lower levels, it is clear that numerical breakup may appear soon behind the rim in no-MD simulations.

## Appendix E. Determination of $d_c$ for MD statistics

The determination of $d_c$, that is, the diameter marking the lower bound of the converged region is done by visual inspection of the histogram showing the droplet-size distribution. Figure 27 shows a zoom in explaining how we determine $d_c$ for each consecutive level. Note that the black curve in the no-MD case is not overlapping in the converged region of the various MD simulations. This is an important fact since it means that an MD simulation is not ‘equivalent’ to a no-MD simulation with a similar ‘effective’ critical thickness $h_c$. Hence, it makes more sense to compare the MD plots to each other and treat the black no-MD curves only for the purpose of analysing the properties of classical non-converging simulations.

## Appendix F. Forensics of the Chirco et al. (2022) study

The MD method appeared first in the work of Chirco et al. (2022). Here we revisit this previous study with the droplet-size distribution comparison method described above. In Chirco et al. (2022) a region of convergence was identified at $8\Delta_9$, where $\Delta_9$ was the grid size of the coarsest simulation, also used as the MD level. In our study, we have shown that the region of statistical convergence can be reached up to the critical sheet thickness
<figure><img src="image_33.png" /><figcaption>Figure 27. A zoom in for the converged region corresponding to figure 16 showing how $d_c$ (indicated by vertical lines) is obtained. Panels (a) and (b) correspond to a zoom in for figure 16(a) and 16(b), respectively.</figcaption></figure>

<figure><img src="image_34.png" /><figcaption>Figure 28. Replotting figure 12 of Chirco *et al.* (2022) with additional vertical lines showing the critical thickness for MD sheet perforation threshold $h_c$, along with two vertical lines for $d_c$. The shaded regions illustrate the departure of the measured distribution from each other below $d_c$.</figcaption></figure>
of the MD threshold $h_c$, given that the mesh is refined sufficiently. In figure 12 of Chirco *et al.* (2022), the authors argue that above $d_c = 8\Delta\eta$, the simulations seem to converge, where the coarsest mesh had the equivalent of $512^3$ cells ($\ell = 9$) and the finest mesh was the equivalent of $2048^3$ cells ($\ell = 11$). The case that the authors simulated was a phase inversion case that produces a much smaller number of droplets than our pulsed jet case does. This implies that the number frequencies are affected by a large statistical error. Hence, the conclusions presented in Chirco *et al.* (2022) are qualitative: they demonstrate the rupture of the sheet by MD in a clear manner. To attempt a more quantitative analysis in line with our analysis of the pulsed jet, we replot figure 12 of Chirco *et al.* (2022) in our figure 28. There we mark the MD threshold $h_c = 3\Delta\eta$ and two values of $d_c$. We see that $d_c$ again moves towards $h_c$. This movement is in line with our results. Due to computational cost, Chirco *et al.* (2022) stopped at the $2048^3$ mesh and concluded that the boundary of the converged region lies at $d_c \sim 8\Delta\eta$. It is likely that the converged region could be pushed more towards $h_c$ with increased resolution.

REFERENCES

AGBAGLAH, G.G. 2021 Breakup of thin liquid sheets through hole-hole and hole-rim merging. *J. Fluid Mech.* **911**, A23.

ANEZ, J., AHMED, A., HECHT, N., DURET, B., REVEILLON, J. & DEMOULIN, F.X. 2019 Eulerian-Lagrangian spray atomization model coupled with interface capturing method for diesel injectors. *Intl J. Multiphase Flow* **113**, 325-342.

BALACHANDAR, S., ZALESKI, S., SOLDATI, A., AHMADI, G. & BOURGUIBA, L. 2020 Host-to-host airborne transmission as a multiphase flow problem for science-based social distance guidelines. *Intl J. Multiphase Flow* **132**, 103439.

BEN RAYANA, F., CARTELLIER, A. & HOPFINGER, E. 2006 Assisted atomization of a liquid layer: investigation of the parameters affecting the mean drop size prediction. In *Proceedings of ICLASS 2006, August 27- September 1, Kyoto, Japan*. Academic Publishing and Printings,

BIANCHI, G.M., MINELLI, F., SCARDOVELLI, R. & ZALESKI, S. 2007 3D large scale simulation of the high-speed liquid jet atomization. *SAE Tech. Pap.* 2007-01-0244.

BRUJAN, E. 2010 *Cavitation in Non-Newtonian Fluids: with Biomedical and Bioengineering Applications*. Springer Science & Business Media.

CHESNEL, J., MÉNARD, T., RÉVEILLON, J. & DEMOULIN, F.-X. 2011 Subgrid analysis of liquid jet atomization. *Atomiz. Sprays* **21** (1), 41-67.

CHIRCO, L., MAAREK, J., POPINET, S. & ZALESKI, S. 2022 Manifold death: a volume of fluid implementation of controlled topological changes in thin sheets by the signature method. *J. Comput. Phys.* **467**, 111468.

CHORIN, A.J. 1968 Numerical solution of the Navier-Stokes equations. *Maths Comput.* **22** (104), 745-762.

CULICK, F.E.C. 1960 Comments on a ruptured soap film. *J. Appl. Phys.* **31** (6), 1128-1129.

DEBAR, R. 1974 Fundamentals of the KRAKEN code. Tech. Rep. UCIR-760.LLNL.

DEBRÉGEAS, G.D., DE GENNES, P.-G. & BROCHARD-WYART, F. 1998 The life and death of ‘bare’ viscous bubbles. *Science* **279** (5357), 1704-1707.

DUKE, D.J. *et al.* 2017 Internal and near nozzle measurements of Engine Combustion Network ‘Spray G’ gasoline direct injectors. *Expl Therm. Fluid Sci.* **88**, 608-621.

FUSTER, D., BAGUÉ, A., BOECK, T., LE MOYNE, L., LEBOISSETIER, A., POPINET, S., RAY, P., SCARDOVELLI, R. & ZALESKI, S. 2009 Simulation of primary atomization with an octree adaptive mesh refinement and VOF method. *Intl J. Multiphase Flow* **35** (6), 550-565.

FUSTER, D., MATAS, J.P., MARTY, S., POPINET, S., HOEPFFNER, J., CARTELLIER, A. & ZALESKI, S. 2013 Instability regimes in the primary breakup region of planar coflowing sheets. *J. Fluid Mech.* **736**, 150-176.

GAVRILOV, L.R. 1969 On the size distribution of gas bubbles in water. *Sov. Phys. Acoust.* **15** (1), 22-24.

GOROKHOVSKI, M. & HERRMANN, M. 2008 Modeling primary atomization. *Annu. Rev. Fluid Mech.* **40** (1), 343-366.

HERRMANN, M. 2010 A parallel Eulerian interface tracking/Lagrangian point particle multi-scale coupling procedure. *J. Comput. Phys.* **229** (3), 745-759.

HERRMANN, M. 2011 On simulating primary atomization using the refined level set grid method. *Atomiz. Sprays* **21** (4), 283-301.
Y. Kulkarni, C. Pairetti, R. Villiers, S. Popinet and S. Zaleski

HIRT, C.W. & NICHOLS, B.D. 1981 Volume of fluid (VOF) method for the dynamics of free boundaries. *J. Comput. Phys.* **39** (1), 201–225.

JACKIW, I.M. & ASHGRIZ, N. 2022 Prediction of the droplet size distribution in aerodynamic droplet breakup. *J. Fluid Mech.* **940**, A17.

JARRAHBASHI, D. & SIRIGNANO, W.A. 2014 Vorticity dynamics for transient high-pressure liquid injection. *Phys. Fluids* **26** (10).

JARRAHBASHI, D., SIRIGNANO, W.A., POPOV, P.P. & HUSSAIN, F. 2016 Early spray development at high gas density: hole, ligament and bridge formations. *J. Fluid Mech.* **792**, 186–231.

JIANG, D. & LING, Y. 2021 Impact of inlet gas turbulence on the formation, development and breakup of interfacial waves in a two-phase mixing layer. *J. Fluid Mech.* **921**, A15.

KANT, P., PAIRETTI, C., SAADE, Y., POPINET, S., ZALESKI, S. & LOHSE, D. 2023 Bag-mediated film atomization in a cough machine. *Phys. Rev. Fluids* **8** (7), 074802.

KHANWALE, M.A., SAURABH, K., ISHII, M., SUNDAR, H. & GANAPATHYSUBRAMANIAN, B. 2022 Breakup dynamics in primary jet atomization using mesh-and interface-refined Cahn–Cilliard Navier–Stokes. arXiv e-prints, arXiv:2209.

KIBBLE, T.W.B. 1976 Topology of cosmic domains and strings. *J. Phys. A: Math. Gen.* **9** (8), 1387–1398.

LASHERAS, J.C. & HOPFINGER, E.J. 2000 Liquid jet instability and atomization in a coaxial gas stream. *Annu. Rev. Fluid Mech.* **32** (1), 275–308.

LEBAS, R., MENARD, T., BEAU, P.A., BERLEMONT, A. & DEMOULIN, F.X. 2009 Numerical simulation of primary break-up and atomization: DNS and modelling study. *Intl J. Multiphase Flow* **35** (3), 247–260.

LING, Y., FUSTER, D., ZALESKI, S. & TRYGGVASON, G. 2017 Spray formation in a quasiplanar gas-liquid mixing layer at moderate density ratios: A numerical closeup. *Phys. Rev. Fluids* **2** (1), 014005.

LOHSE, D. & VILLERMAUX, E. 2020 Double threshold behavior for breakup of liquid sheets. *Proc. Natl Acad. Sci. USA* **117** (32), 18912–18914.

MARCOTTE, F. & ZALESKI, S. 2019 Density contrast matters for drop fragmentation thresholds at low Ohnesorge number. *Phys. Rev. Fluids* **4** (10), 103604.

MÉNARD, T., TANGUY, S. & BERLEMONT, A. 2007 Coupling level set/VOF/ghost fluid methods: validation and application to 3D simulation of the primary break-up of a liquid jet. *Intl J. Multiphase Flow* **33** (5), 510–524.

OPFER, L., ROISMAN, I.V., VENZMER, J., KLOSTERMANN, M. & TROPEA, C. 2014 Droplet-air collision dynamics: evolution of the film thickness. *Phys. Rev. E* **89** (1), 013023.

PAIRETTI, C., VILLIERS, R. & ZALESKI, S. 2021 On shear layer atomization within closed channels: Numerical simulations of a cough-replicating experiment. *Comput. Fluids* **231**, 105125.

PAIRETTI, C.I., DAMIAN, S.M., NIGRO, N.M., POPINET, S. & ZALESKI, S. 2020 Mesh resolution effects on primary atomization simulations. *Atomiz. Sprays* **30** (12), 913–935.

POPINET, S. 2003 Gerris: a tree-based adaptive solver for the incompressible Euler equations in complex geometries. *J. Comput. Phys.* **190** (2), 572–600.

POPINET, S. 2009 An accurate adaptive solver for surface-tension-driven interfacial flows. *J. Comput. Phys.* **228** (16), 5838–5866.

POPINET, S. 2018 Numerical models of surface tension. *Annu. Rev. Fluid Mech.* **50** (1), 49–75.

SALVADOR, F.J., RUIZ, S., CRIALESI-ESPOSITO, M. & BLANQUER, I. 2018 Analysis on the effects of turbulent inflow conditions on spray primary atomization in the near-field by direct numerical simulation. *Intl J. Multiphase Flow* **102**, 49–63.

SAURABH, K., ISHII, M., KHANWALE, M.A., SUNDAR, H. & GANAPATHYSUBRAMANIAN, B. 2023 Scalable adaptive algorithms for next-generation multiphase flow simulations. In *2023 IEEE International Parallel and Distributed Processing Symposium (IPDPS)*, pp. 590–601. IEEE.

SAVVA, N. & BUSH, J.W.M. 2009 Viscous sheet retraction. *J. Fluid Mech.* **626**, 211–240.

SHIMA, A. & SAKAI, I. 1987 On the equation for the size distribution of bubble nuclei in liquids (Report 2). Rep. Inst. High Speed Mech. Tohoku Univ, 54, 51–59.

SHINJO, J. & UMEMURA, A. 2010 Simulation of liquid jet primary breakup: dynamics of ligament and droplet formation. *Intl J. Multiphase Flow* **36** (7), 513–532.

SONG, M. & TRYGGVASON, G. 1999 The formation of thick borders on an initially stationary fluid sheet. *Phys. Fluids* **11** (9), 2487–2493.

TANG, K., ADCOCK, T.A.A. & MOSTERT, W. 2023 Bag film breakup of droplets in uniform airflows. *J. Fluid Mech.* **970**, A9.

TAYLOR, G.I. 1959 The dynamics of thin sheets of fluid. III. Disintegration of fluid sheets. *Proc. R. Soc. Lond. A: Math. Phys. Sci.* **253** (1274), 313–321.

TOLFTS, O., DEPLUS, G. & MACHICOANE, N. 2023 Statistics and dynamics of a liquid jet under fragmentation by a gas jet. *Phys. Rev. Fluids* **8** (4), 044304.
TORREGROSA, A.J., PAYRI, R., SALVADOR, F.J. & CRIALESI-ESPOSITO, M. 2020 Study of turbulence in atomizing liquid jets. *Intl J. Multiphase Flow* **129**, 103328.

VERNAY, C., RAMOS, L., WÜRGER, A. & LIGOURE, C. 2017 Playing with emulsion formulation to control the perforation of a freely expanding liquid sheet. *Langmuir* **33** (14), 3458–3467.

VILLERMAUX, E. 2020 Fragmentation versus cohesion. *J. Fluid Mech.* **898**, P1.

VILLERMAUX, E. & BOSSA, B. 2009 Single-drop fragmentation determines size distribution of raindrops. *Nat. Phys.* **5** (9), 697–702.

WANG, Y. & BOUROUIBA, L. 2018 Unsteady sheet fragmentation: droplet sizes and speeds. *J. Fluid Mech.* **848**, 946–967.

WEYMOUTH, G.D. & YUE, D.K.P. 2010 Conservative volume-of-fluid method for free-surface simulations on Cartesian-grids. *J. Comput. Phys.* **229** (8), 2853–2865.

ZHANG, B., POPINET, S. & LING, Y. 2020 Modeling and detailed numerical simulation of the primary breakup of a gasoline surrogate jet under non-evaporative operating conditions. *Intl. J. Multiphase Flow* **130**, 103362.

ZUREK, W.H. 1985 Cosmological experiments in superfluid helium? *Nature* **317** (6037), 505–508.