Experimental observations and modeling of sub-Hinze bubble production by turbulent bubble break-up

We present experiments on large air cavities spanning a wide range of sizes relative to the Hinze scale $d_\mathrm{H}$, the scale at which turbulent stresses are balanced by surface tension, disintegrating in turbulence. For cavities with initial sizes $d_0$ much larger than $d_\mathrm{H}$ (probing up to $d_0 / d_\mathrm{H} = 8.3$), the size distribution of bubbles smaller than $d_\mathrm{H}$ follows $N(d) \propto d^{-3/2}$, with $d$ the bubble diameter. The capillary instability of ligaments involved in the deformation of the large bubbles is shown visually to be responsible for the creation of the small ones. Turning to dynamical, three-dimensional measurements of individual break-up events, we describe the break-up child size distribution and the number of child bubbles formed as a function of $d_0 / d_\mathrm{H}$. Then, to model the evolution of a population of bubbles produced by turbulent bubble break-up, we propose a population balance framework in which break-up involves two physical processes: an inertial deformation to the parent bubble that sets the size of large child bubbles, and a capillary instability that sets the size of small child bubbles. A Monte Carlo approach is used to construct the child size distribution, with simulated stochastic break-ups constrained by our experimental measurements and the understanding of the role of capillarity in small bubble production. This approach reproduces the experimental time evolution of the bubble size distribution during the disintegration of large air cavities in turbulence.


Introduction
1.1. Broader context Gas bubbles dispersed in liquids provide surface area through which mass can be exchanged by diffusion. Ocean-atmosphere exchanges of CO 2 , for example, are enhanced by bubble-mediated transfer in regions of the globe where high winds lead to high rates of wave breaking, as entrained air cavities break apart into small bubbles in the turbulent field under the breaking wave (Deike & Melville 2018;Reichl & Deike 2020;Deike 2022). Further, many industrial processes involve facilitating gas transfer to a liquid through bubble interfaces (Schludieter et al. 2021). In both environmental and industrial scenarios, the breakage of bubbles by the turbulence of the bulk flow increases the total surface area through which transfers may occur and modulates the bubbles' dynamics.
Despite the ubiquity of bubble break-up across disciplines, the physics of bubble breaking in turbulence remains to be fully understood, as turbulent effects are often accompanied by buoyant effects and shear in the mean structure of the flow (Risso & Fabre 1998). Further, the fast dynamics of bubble pinching have, until recently, been difficult to measure experimentally, leaving open questions regarding the final portion of the break-up process (Ruth et al. 2019). These various challenges have led to a wide variability in the predictions of models for both the rate at which bubbles break and the sizes of bubbles they break into.

Bubble break-up in turbulence
We consider the break-up of a bubble with an effective diameter d 0 , taken to be the diameter of a sphere with the same volume. Before considering the turbulent nature of the liquid around it, the bubble in a liquid is described by the density of the liquid and gas phases, ρ and ρ g , their viscosities µ and µ g , the acceleration due to gravity g, and the surface tension of the liquid-gas interface σ. When the carrier flow in which the bubbles are dispersed (with velocity u) is turbulent, it is characterized by the dissipation rate of the turbulence , which is the rate at which kinetic energy in turbulent fluctuations is dissipated to heat. The turbulence is comprised of fluctuating motions existing over a range of length scales, extending from larger motions near the integral length scale L int (beyond which the velocity field becomes uncorrelated) down to the Kolmogorov scale η, at which turbulent motions are dissipated by the viscosity of the fluid (Pope 2000).
With nine independent physical parameters which span three physical dimensions, we require six dimensionless parameters to describe the problem of bubble break-up in turbulence, for which we choose where the subscript "0" indicates a quantity refers to an initial condition. The size of the parent bubble relative to the capillary length scale l cap = σ/(ρg) describes the relative importance of gravity and surface tension effects for the parent bubble. The large-scale turbulence Reynolds number Re t represents the separation of length scales in the turbulence. The bubble size relative to the integral length scale d 0 /L int , along with Re t , describes the spatial separation between the bubble and the turbulence scales. With Re t 1 and ρ/ρ gas and µ/µ gas both fixed constants 1 for common liquid-gas configurations, we will neglect their impact in the rest of the experimental study. The Weber number of the parent bubble We 0 , which parameterizes the balance between turbulent stresses and surface tension, will be the main parameter of focus.
For a bubble in the inertial subrange of the turbulence (η d 0 L int ), the ratio of the inertial stresses arising from velocity gradients in the turbulence and surface tension stresses defines the Weber number, We(d) = C 2 ρ 2/3 d 5/3 /σ, with C 2 = 2, and is central in the analysis of bubble break-up (Risso & Fabre 1998;Rivière et al. 2021;Perrard et al. 2021). The definition of a critical Weber number for break-up We c yields the Hinze scale (Hinze 1955), and we typically use the ratio d/d H = (We/We c ) 3/5 in place of We. Estimations of We c vary, and generally involve either considerations of how likely a bubble is to break apart over some physically-relevant time or within some spatial observation window (Hinze 1955;Martínez-Bazán et al. 1999b;Risso & Fabre 1998;Rivière et al. 2021), or considerations of the shape of the bubble size distribution resulting from break-ups (Deane & Stokes 2002). Since We c is influenced by factors like the buoyancy and specificity of the turbulent flow, and since the turbulent stresses on a bubble are stochastic in nature, the Hinze scale as defined in eq. (1.2) represents a soft limit for break-up. Different experimental and computational setups will lead to a range of reported or inferred critical Weber numbers, which typically vary from 1-5 Risso & Fabre 1998;Hinze 1955;Martínez-Bazán et al. 1999b;Vejrazka et al. 2018). In this paper, we will use We c = 1, consistent with our results and similar experiments in a turbulent flow forced by underwater pumps (Vejrazka et al. 2018). We note that the inertial stresses on a bubble that arise from the velocity slip between the bubble and the surrounding liquid can induce stresses comparable to those associated with the turbulence's inherent velocity gradients at the bubble scale (Masuk et al. 2021), that eddies smaller than the bubble can also contribute to deformation and break-up (Luo & Svendsen 1996;Qi et al. 2022), and that the turbulent flow can trigger bubble shape oscillations (Risso & Fabre 1998;Ravelet et al. 2011). These factors will contribute to bubble deformation and break-up in ways that are not directly parameterized in the definition of d H . The bubble size distribution N (d) gives the number density of bubbles with diameter d, and given the nature of experiments reported in this paper, we define it such that N (d)dd gives the total number of bubbles with diameters ∈ (d, d + dd). Garrett et al. (2000) proposed that, for bubbles larger than the Hinze scale, a power-law scaling N (d) ∝ d −10/3 describes the steady-state bubble size distribution, assuming that the break-up rate scales with the turbulent frequency at the bubble size. This regime has since been reported in several experiments (Deane & Stokes 2002;Blenkinsopp & Chaplin 2010;Rojas & Loewen 2007) and simulations (Deike et al. 2016;Wang et al. 2016;Gao et al. 2021;Chan et al. 2021;Soligo et al. 2019;Rivière et al. 2021;Mostert et al. 2022). For smaller bubbles, the size distribution typically exhibits a shallower slope (Deane & Stokes 2002;Blenkinsopp & Chaplin 2010), with fewer studies resolving this range of scales and some variation in the values that have been reported. The N (d) ∝ d −3/2 distribution for d < d H has been observed experimentally (Deane & Stokes 2002) and numerically (Wang et al. 2016;Mostert et al. 2022) for bubbles under breaking waves, though the identification of a sub-Hinze power-law slope is additionally complicated by the transient nature of bubble disintegration ) and breaking wave  events. Recent work has identified the capillary pinching of gas ligaments created by turbulent deformations as an origin of sub-Hinze bubbles, with theoretical arguments relating to the timescale over which such pinching occurs supporting the N (d) ∝ d −3/2 sub-Hinze scaling (Rivière et al. 2022). Relating measured size distributions to theoretical scalings derived from break-up physics is complicated by the fact that bubbles' motions, and hence their residence time in some experimental domain, are dependent on their size and the characteristics of the turbulence they encounter (Garrett et al. 2000). Smaller bubbles or bubbles in regions of more intense turbulence will rise slower than others, for example (Ruth et al. 2021); accounting for these effects requires detailed knowledge of the size dependencies of the bubbles' motions.

Child size distribution and break-up time scales
In this work, we will employ experimental observations to describe bubble breakup over a range of spatial scales: we consider parent bubbles ranging in size from the Hinze scale to d 0 = 8.3d H , and investigate how they break up to produce child bubbles that may be orders of magnitude smaller than the Hinze scale. As volume is conserved in any break-up, we will work with bubble volumes V = πd 3 /6 when discussing bubble break-up, denoting parent bubble volumes by V = ∆ and child bubble volumes by V = δ.
Expressions for a break-up kernel f (δ; ∆), for which f (δ; ∆)dδ gives the rate at which a parent bubble of volume ∆ will break into a child bubble with volume ∈ (δ, δ + dδ) in some turbulent flow, are informed by experiments and simulations on break-up. Most experimental studies have involved air bubbles in water under Earth's gravitational acceleration, with turbulence in the water generated by one or more jets (Martínez-Bazán et al. 1999b;Vejrazka et al. 2018;Qi et al. 2020), rotating blades (Ravelet et al. 2011), or by turbulent flow through a reactor or channel (Andersson & Andersson 2006). Risso & Fabre (1998) performed experiments on bubble break-up in microgravity to remove the effects of buoyancy, which also contributes to bubble deformation and break-up, and more recently, Rivière et al. (2021) performed direct numerical simulations (DNSs) of bubble break-up without gravity, solving the full two-phase Navier Stokes equations for a bubble subjected to homogeneous, isotropic turbulence.
These studies have confirmed that the time over which a break-up occurs is controlled by both the turbulent scales and the bubble's oscillatory scales. Rivière et al. (2021) showed that, as a bubble of size d 0 d H is introduced to turbulence, it first breaks up after a time comparable to eddy turn-over time at its scale, T turb (d 0 ) = −1/3 d 2/3 0 . Experimental studies have shown that the time over which deformation occurs prior to break-up scales similarly (Qi et al. 2020;Risso & Fabre 1998). As the deformation of moderately-sized bubbles is also impacted by the surface tension, capillary dynamics remain important, as a bubble's natural oscillation frequency remains apparent in its shape oscillations (Risso & Fabre 1998;Ravelet et al. 2011;Perrard et al. 2021). Further, the turbulent turnover time is typically comparable to the capillary oscillation time at the parent bubble scale for air bubbles in water at moderate d 0 /d H , which can lead to a resonance which aides break-up (Risso & Fabre 1998;Ravelet et al. 2011).
The break-up frequency ω is defined as the inverse of the typical time until Focus on Fluids articles must not exceed this page length a bubble undergoes a break-up, and is distinct from the (necessarily shorter) typical duration over which a break-up occurs. Ravelet et al. (2011) showed that the distribution of the times until a bubble breaks mirrors the distributions of the times between severe shape deformations and the times between large instantaneous Weber numbers. The most energetic scales capable of deforming a bubble are those at the scale of the bubble, and experiments from which ω was extracted suggested that the break-up frequency initially increases with bubble size as the turbulence becomes more capable of counteracting surface tension, and then decreases for even larger bubbles, as the time required for a turbulent eddy to act across the bubble scale becomes longer (Martínez-Bazán et al. 1999b), though this analysis may have missed break-ups in which one child bubble size is close to the parent size (Lehr et al. 2002). Recent experiments from Qi et al. (2022) showed that eddies smaller than d 0 can also cause break-up, and other theoretical analyses have considered the action of a range of turbulent scales which may cause break-up. In such models, the product of the rate at which eddies of a given size interact with a bubble and each interaction's likelihood of causing break-up are integrated over a range of eddy sizes (Prince & Blanch 1990;Tsouris & Tavlarides 1994;Luo & Svendsen 1996;Lehr et al. 2002;Aiyer et al. 2019;Yuan et al. 2021), causing the break-up frequency to increase with the bubble size as more turbulent scales contribute to break-up.
Various models for the child size distributions p(δ; ∆) have been proposed, most of which assume that each break-up produces two bubbles. The child size distribution has been described with a ∩-shaped dependence on δ-that is, the most likely outcome is to produce child bubbles that are comparable in size to the parent bubble (Martínez-Bazán et al. 1999a; or with a ∪-or W-shaped child size distributions, in which small bubbles are more likely to be produced than moderately-sized ones (Qi et al. 2020;Rivière et al. 2021;Vejrazka et al. 2018;Andersson & Andersson 2006;Tsouris & Tavlarides 1994;Luo & Svendsen 1996;Lehr et al. 2002;Yuan et al. 2021;Qi et al. 2020). Experimental and numerical evidence suggests that break-ups often produce just two child bubbles when d 0 /d H is close to 1 (Vejrazka et al. 2018;Rivière et al. 2021). However, break-ups at larger d 0 /d H are more severe and often result in more than two child bubbles being formed in a single coherent event (Vejrazka et al. 2018;Hinze 1955;Rivière et al. 2021). Hill & Ng (1996) developed generalized expressions for p(δ; ∆) as products of power-law relations (each ∝ δ α ) for α > −1 and integer numbers of child bubbles, which by design satisfy constraints relating to the sizes of the bubbles formed. Their analysis was extended to break-ups with a non-integer average number of child bubbles by Diemer & Olson (2002).
In the work discussed so far, the role of capillarity has been to counteract the turbulent stresses and prevent severe deformation, while also providing a resonance mechanism at moderate d 0 /d H . However, more recent work has shown that capillarity also plays an important role late in the break-up process, even after a turbulent stress has decidedly overcome it. Andersson & Andersson (2006) showed that asymmetries in a deformed bubble shape can become more pronounced as a bubble breaks apart due to the variation in capillary pressure associated with the deformation. More recently, Rivière et al. (2022) showed that very small bubbles originate not from turbulent motions at very small scales, but rather from the capillary instabilities of ligaments arising from much larger-scale deformations.
1.4. Outline of the paper In this work we address the problem of bubbles breaking up in forced turbulence, which is applicable to break-up under breaking waves and in industrial reactors. We probe a wide range of scales, with bubbles ranging in size from the Hinze scale to d = 8.3d H (corresponding to We 0 = 34). Further, we resolve the size distribution down to approximately an order of magnitude smaller than d H , enabling us to identify the way in which the sub-Hinze size distribution scales when there is a large separation between the Hinze scale and the bubbles which break.
The experiment set-up, including the turbulence generation, is detailed in Section 2. The results on the disintegration of large air cavities are given in Section 3, spanning a wide range of d 0 /d H . We demonstrate experimentally that a N (d) ∝ d −3/2 distribution below the Hinze scale is observed when the initial cavity size is much larger than the Hinze scale, supporting the notion that the capillary pinching dynamics proposed by Rivière et al. (2022) are effective at producing sub-Hinze bubbles. The dynamically-tracked individual bubble breakups with moderate d 0 /d H and resulting child size distributions are discussed in Section 4. In Section 5 we develop a model for turbulent bubble break-up that unifies the turbulent inertial dynamics with the faster, capillary pinching dynamics responsible for sub-Hinze bubble production, ascribing these physical mechanisms to various components of a modeled child size distribution. The model is informed by both experimental observations of the disintegrations of air cavities of various sizes and by experimental and numerical observations of individual break-up events. Concluding remarks are given in Section 6.

Experimental setup
This paper presents the results of two separate, complementary experiments, both involving air bubbles breaking apart in forced water turbulence. In the first, we generate large cavities of air with sizes much larger than the Hinze scale (with d 0 /d H betwen 2.1 and 8.3) and measure the transient evolution of the bubble size distribution as the cavity disintegrates in successive break-ups. In the second experiment, we introduce moderately-sized bubbles (with d 0 /d H between 1 and 3) into the turbulence, and track the outcomes of their individual break-ups. The turbulence generation is identical in both set-ups.

Turbulence generation and characterization
Turbulence in a 0.37 m 3 water tank is generated by the convergence of eight turbulent jets created by four submerged water pumps, as sketched in Figure 1 (a) and described in greater detail in Ruth et al. (2021). The flow from each pump is split into two parallel jets at a Y, with each outlet separated by 7.8 cm, with the centers of the Y forming the vertices of a 25 cm square in the horizontal plane. Figure 1 (b) presents properties of the flow as characterized in the central plane (y = 0) of the experiment with two-dimensional, two-component particle image velocimetry (PIV). The background gives the local fluctuation velocity u = (u x 2 + u z 2 )/2, where u i = (u i − u i ) 2 and overbars denote averaging in time. u tends to be largest in the plane of the jets (z ≈ 0.01 cm) and in the region below their convergence zone (x ≈ y ≈ 0). PIV is performed in nine parallel planes, enabling the three-dimensional interpolation of turbulence quantities at any location within the measurement domain. As described in Ruth et al. (2021), we compute the integral length scale L int locally at each point in the flow by integrating the spatial autocorrelation function. It changes throughout the experiment, being the shortest where the turbulence is the strongest. The cyan lines in Figure 1 (b) denote the value of L int at various locations in the central plane of the experiment: L int is shortest near the convergence of the jets, and grows at lower and higher depths. With u and L int calculated from the PIV data, we can then compute the local turbulence dissipation rate under the assumption of isotropy with = C u 3 /L int , with C = 0.7 (Sreenivasan 1998), and the Kolmogorov microscale with η = ((µ/ρ) 3 / ) 1/4 (Pope 2000). The Hinze scale d H , calculated using eq. (1.2), is denoted at various locations by the diameter of the black circles drawn in Figure 1 (b). The Hinze scale is smaller where the turbulence is more intense, meaning that more bubbles will be larger than the Hinze scale and susceptible to break-up at these locations. We refer to Ruth et al. (2021) for more details on the structure of the turbulence field and for maps of turbulent quantities outside of the central plane.

Large cavity disintegration experiment
For the experiment on large cavity break-ups, air cavities were produced following Landel et al. (2008) by placing a hollow hemispherical cup with R = 5 cm underwater, sketched in Figure 2 (a), and bubbling a known volume of air The characteristic length scales η, dH, Lcap, and Lint taken in analyzing the data, the pixel size ∆x and the minimum bubble size considered dmin, and the diameters of the air cavities studied (circles). Distributions of the turbulence quantities in the field of view in the center of the tank (within the green rectangle in Figure 1) are also given in gray.
V 0 = πd 3 0 /6 into it. Once bubbles in this cup have coalesced into a single air cavity, the cup is then inverted by rotating it rapidly half a revolution, such that the air inside is suddenly no longer constrained by the curved cup surface. The top surface of the initial volume of air roughly conforms to the curved inner surface of the cup. The large air cavity, having been suddenly exposed to stresses from the surrounding turbulence and its buoyant rise through the water, deforms and starts a complex sequence of break-ups, leading to its disintegration. The surface of the cup rotates with a speed around 0.4 m/s to 0.9 m/s, and we have checked that this speed does not systematically impact the early stages of the bubble size distribution. Further, similar experiments run without turbulence yield very little break-up, as evidenced in Appendix B.
The turbulent flow in the region of the tank imaged in this experiment is denoted by the green rectangle in Figure 1 (b). The turbulence varies spatially, so to simplify the analysis, we take u ≈ 0.2 m/s, L int ≈ 1.5 cm, and η ≈ 37 µm as characteristic values, which set d H = 3.2 mm and Re t = 3400. These length scales are denoted in Figure 2 (c), which also gives the distribution of the length scales present in the field of view in the middle of the tank. The mean flow is downwards with W ≈ −0.25 m/s, largely counteracting the buoyant rise speed of larger bubbles. This enables the bubble population to linger in the measurement region for a sufficient period of time to image it over multiple large-scale eddy turn-over times T int = L int /u ≈ 0.075 s. The cavities range in size between d 0 /d H = 2.12 and 8.30. Data for each condition, as well as the number of runs recorded at each, are given in Table 1.
One image is shown in Figure 2 (b). The cup is visible in the bottom of the image as it has not yet fully rotated out of the field of view. The measurement region, which spans 15.8 cm in the x direction and 8.9 cm in the z direction, is illuminated from the back, and the disintegration of the cavity is filmed with a high-speed camera at 500 Hz with a spatial resolution of 38 µm/pixel. The field of view is much larger than all the bubbles considered, so it does not introduce a significant bias related to bubbles whose images extend partially outside the field of view. Bubbles are detected with an image processing method described in Appendix A.1, and their effective diameters d are determined as the equivalent diameter of a circle with the same area as the projected bubble image.
In analyzing the data, we consider only bubbles for which d 400 µm, for which the detection is less sensitive to the chosen image intensity threshold. Given the typical severe deformation and overlapping images of larger bubbles (d 6 mm), we note that their sizes will tend to be over-estimated by this method. The air void fraction in the vicinity of the cavity is high enough that we are unable to track the dynamics of individual break-ups, so we restrict our analysis to the resulting bubble size distribution.
To account for the limited field of view in our experiments, we adjust the measured size distributions by keeping a record of bubbles which have left and entered the field of view. Those which leave are "locked" into the bubble record used in computing the size distributions, while those that enter the field of view are excluded from the calculation of the size distribution. This, along with a slight smoothing in d and t to account for the limited number of bubbles observed at early times or with small cavities, is described in greater detail in Appendix A.2, and has only a limited impact on the results reported in this paper, as we do not consider the size distribution at late times.

Individual break-up tracking experiment
In the second set of experiments of bubble break-up, we dynamically track the individual break-ups of bubbles in the turbulent region. As sketched in Figure 3 (a), bubbles are introduced to the bottom of the tank through a needle and rise to the turbulent region. Two cameras, which are synchronized with a function generator, film at 1000 fps. They are oriented 90°from each other and their fields of view overlap in a measurement volume of approximately 200 cm 3 . The cameras are calibrated by mapping their pixels to the paths of the light rays reaching the pixels, following the method presented by Machicoane et al. (2019). Then, following a method similar to that used in Ruth et al. (2021), we identify the 3-D location of the bubbles which are simultaneously captured by each camera. The spatial resolution of each camera varies with the position of the bubble, but the typical value of the two cameras are 28.9 µm/pixel and 57.1 µm/pixel. An approximate lower bound for the size of the smallest resolved bubble is then d min ≈ 200 µm. The trajectories taken by the bubbles are then determined using the Python package Trackpy (Allen et al. 2021), which implements the algorithm from Crocker & Grier (1996). Such trajectories are shown in Figure 3 (b). Using the three-dimensional map of the turbulence statistics obtained from PIV, we compute the bubble's size relative to the local Hinze scale d/d H (computed with the local value of ) at each bubble location, which is encoded in the color in the figure. The mean dissipation rate at the break-up locations is = 0.52 m 2 /s 3 , with a standard deviation of 0.21 m 2 /s 3 . The mean values and standard deviations of quantities describing the initial conditions for the break-ups studied in this experiment are given in Table 1.
From the bubble trajectories, we identify each time a bubble breaks apart, which occurs when a new trajectory (or trajectories) appears in the vicinity of a previously-existing bubble. As the tracking algorithm will initially link the parent bubble to only one of the child bubbles, the parent bubble trajectory is then split at this time, and both child bubbles are treated equally. These events are denoted by the gray lines connecting the "end" of one bubble to the "beginning" of another in Figure 3 (b). Given the complex deformations involved in some break-ups, the method does not always resolve the fast splitting dynamics accurately; the breakup child size distributions we report, however, are not sensitive to the order of events occurring within one break-up event. 3. Size distribution evolution during the disintegration of a large air cavity Here, we present experimental results on the disintegration of air cavities of various sizes from the experiment described in Section 2.2. First, we qualitatively discuss the break-up of cavities in two illustrative cases, one close to the critical size for break-up, and one with a large separation of scales between the cavity and the Hinze scale. Then, we analyze the transient evolution of the bubble size distributions.

Disintegration of cavities of increasing sizes
The break-ups of two air cavities, one with d 0 = 0.68 cm and one with d 0 = 2.3 cm, are shown in Figure 4 and Figure 5, respectively. These correspond to nondimensional sizes of d 0 /d H = 2.1 and 7.0, d 0 /L int = 0.46 and 1.5, and d 0 /l cap = 2.51 and 8.28. As a reference, the constant values taken for L int and d H and the initial size of the cavity d 0 are denoted in the top-left corner of the first image. In both cases, the hemispherical cup which had constrained the bubble is visible at early times as it is rotated away.
In the disintegration of the smaller cavity, with d 0 /d H = 2.1 (shown in Figure 4), the bubble emerges from the cup with a moderate deformation caused by buoyancy and the surrounding turbulence. Eventually, within approximately one integral-scale turn-over time, the bubble becomes more elongated and breaks into two bubbles. One is near the parent bubble in size, and other is slightly smaller than the Hinze scale. These two bubbles persist without breaking for at least ∼ 2 more integral-scale turn-over times, at which point the smaller of the two bubbles is advected out of the field of view by the downwards mean flow.
The deformation to the larger cavity shown in Figure 5 is more severe, leading to a more complex sequence of events during its disintegration. Upon emerging from the rotating cup, the cavity is flattened due to buoyancy (as d 0 /l cap = 8.28 for this case), and turbulent deformations to the cavity shape on the order of the cavity size itself quickly develop. By t/T int ≈ 0.4, the cavity consists of two lobes (each of which is significantly deformed), separated by a shrinking neck of air. By the time the neck has pinched apart (t/T int ≈ 0.7), the two larger child bubbles stemming from the two lobes are accompanied by much smaller child bubbles (some with d d H and d d 0 ) which were formed during the collapse of the air neck. The larger child bubbles themselves go on to further break apart in a chain of break-ups, some of which similarly involve small bubble production via the collapse of elongated air necks. Many small bubbles which are more than an order of magnitude smaller than the initial one are eventually visible. At much later times, the largest bubbles have risen out of the field of view, and the total air volume imaged is significantly decreased.

Transient evolution of the bubble size distributions
The experiment was carried out with six values of d 0 /d H between 2.1 and 8.3, with 10-20 runs taken at each condition, as given in Table 1. Note that the largest cavities exceed the integral length scale in size, so the typical turbulent stress at their spatial scale will be saturated relative to that predicted by the Kolmogorov scaling employed in the definition of the Hinze scale. Figure 6 shows the transient  Figure 6: Bubble size distributions during the disintegration of cavities with d0/dH between 2.1 and 8.3 and times up to 4Tint after the cavity is released into the turbulence. The size of the parent bubble is denoted by the dashed vertical line. Each distribution integrates to the average number of bubble observed at that condition at that time. The eventual sub-Hinze power-law scaling exponent approaches −3/2, shown by the dashed black line, as d0/dH is increased.
10-20 runs). At early times, the distributions for all d 0 /d H exhibit a peak at d 0 /d H , denoted by the vertical dotted lines. For the two smallest cavities (given in Figure 6 panels (a-b)), for which no break-up was observed during many runs, only a small number of bubbles are formed over time, and the size distribution near the injection scale does not decrease appreciably with time.
Over time, as the larger cavities (given in Figure 6 panels (c-f)) begin to disintegrate, the size distribution for d < d 0 begins to be built up. Even among these larger cavities which produce a considerable number of sub-Hinze bubbles, the increase in the number of sub-Hinze bubbles is much more pronounced for the cavities that are initially larger (evidenced by comparing curves for d 0 /d H = 4.15 and d 0 /d H = 8.30, for example). This suggests that there is a large separation of scales between the sub-Hinze bubbles and the parent bubbles responsible for their creation; more simply, large bubbles are needed for the production of small bubbles. For the largest cavities, the size distribution for sub-Hinze bubbles eventually follows an sketched on all plots as the dashed line for reference. The final curves shown (for t/T int = 4) might constitute an under-estimation for the bubble size distribution for smaller bubbles, since some of the bubbles which may break have risen out of the field of view by this time. Now, we consider the size distributions averaged averaged between 2T int and 4T int . During these times, a significant number of break-ups have occured (for larger d 0 /d H ), but a significant portion of the bubbles have not yet left the field of view, and the bubble size distribution approaches a constant shape.
scaling is approached, indicated by the dashed black line. The size distribution is affected not only by the break-up physics, but is also steepened by the rising dynamics of the bubbles: as small bubbles rise more slowly than larger ones, they tend to linger in the field of view for longer, increasing their concentrations (Garrett et al. 2000).
Integrating the transient size distributions over the bubble diameter, the temporal evolution of the number of resolved bubbles n (with the minimum resolvable size d min = 0.12d H ) is shown in Figure 8 (a). The gray shaded region denotes the times over which the bubble size distributions are averaged in Figure 7 and Figure 8 (b). Figure 8 (b) shows the number of resolved sub-Hinze, super-Hinze, and total bubbles, averaged over the time period considered, with the minimum resolved bubble size d min ≈ 0.12d H . Again, we see an increase in the number of bubbles formed with the initial size of the cavity. Further, the number of sub-Hinze bubbles increases with the parent bubble size more rapidly than the number of super-Hinze bubbles does, making sub-Hinze bubbles constitute a larger portion of the bubble size spectrum for larger d 0 /d H . This is remarkable, since as d 0 /d H is increased, the span of bubble sizes constituting resolvable sub-Hinze bubbles (d min < d < d H ) remains fixed, while the span of potential super-Hinze bubble Taken together, Figures 7 and 8 are congruent with the capillary pinching mechanisms for sub-Hinze bubble production put forward by Rivière et al. (2022). Our figures suggest that their formation relies on the break-up of cavities that power-law scaling in the sub-Hinze bubble size distribution, and the dependence on d 0 of the number of sub-Hinze bubbles produced (shown in Figure 8 (b)) is steeper than that of the number of super-Hinze bubbles produced. We propose in the next section an explanation of the mechanisms leading to this dependence.
3.3. Capillary splitting of ligaments prepared by the turbulence produces small bubbles Visual observations of the large air cavities disintegrating provide clues into the mechanism responsible for the production of sub-Hinze bubbles: child bubbles much smaller than the Hinze scale are seen to originate from a Rayleigh-Plateaulike instability that occurs during the pinching apart of elongated fluid ligaments prepared by the turbulence. However, the turbulence is only able to deform bubbles that are large enough with respect to the Hinze scale that such ligaments might be created, since surface tension is effective at limiting the severity of deformations to smaller bubbles. These experimental observations parallel a recent interpretation of DNSs of bubble break-up (Rivière et al. 2022).
Illustrative examples of bubble break-up are given in Figure 9, which shows the typical break-ups of bubbles of two sizes: one is near the Hinze scale in size (a), and another is seven times larger than the Hinze scale (b). The smaller bubble, with d 0 /d H = 2.1, is initially deformed into two comparably-sized lobes, and the neck separating the two splits at a single point to form two child bubbles, each of a similar scale as the parent. Here, the parent bubble is small enough that surface tension is able to prevent significant deformation during the break-up.
The larger bubble, with d 0 /d H = 7.0, is similarly deformed by the turbulence into two comparably-sized lobes prior to pinch-off. However, the filament of air separating the two just prior to pinch-off has become significantly more elongated inertial deformation inertial deformation capillary instability of the gas ligament pinch-off at a single point Figure 9: Break-up of a bubble initially close to the Hinze scale in size (a, with d0/dH = 2.1) and initially much bigger than the Hinze scale (b, with d0/dH = 7.0). The deformation to the smaller bubble produces two comparably-sized lobes, which split apart to form two comparably-sized child bubbles. The deformation to the larger bubble also produces two comparably-sized lobes, but these are separated by a much more elongated filament of air. The unstable collapse of this filament produces the small "capillary" child bubble between the two larger ones. The small bubbles in the lower left of the image were produced in previous break-ups.
than the neck in the break-up of the smaller bubble. This elongation opens the door to capillary instabilities along the filament during its collapse: in the instance shown in Figure 9 (b), the filament pinches apart at two separate points, leaving a small child bubble (with d d H ) between the two lobes. The two examples of break-up discussed illustrate two mechanisms present in the break-up of bubbles by turbulence. The first is the deformation of the parent bubble by a turbulent structure, likely on the spatial scale of the parent bubble itself. This brings the bubble to an unstable state consisting of two lobes (which will become what we call the "inertial" child bubbles) separated by a neck of air, which begins to pinch apart under capillarity. When the deformation to the bubble is severe enough, this ligament can take on an elongated, deformed shape. Its pinching can become unstable under a Rayleigh-Plateau-like mechanism, leading to the formation of small "capillary" bubbles. Figure 10 shows an additional five instances of deformed ligaments undergoing a capillary instability to produce sub-Hinze bubbles. Many cases, especially those involving large parent bubbles, do not solely involve one ligament separated by two lobes; turbulent deformations cause the bubble shapes to be more irregular. However, in all instances, very small bubbles are produced as an air ligament involved in the turbulent deformation collapses unstably.
This description is clearly a simplified understanding of the bubble pinchoff process, as it does not capture the redistribution of air due to a capillary pressure difference between lobes that may be responsible for the formation of small bubbles (Andersson & Andersson 2006), nor does it describe the "tearing off" of very small bubbles that we observe occurring to large parent bubbles more frequently than 1/T turb (d 0 ). However, the framework serves as a bridge between the inertial deformations to a bubble by the turbulence and the latertime collapse dynamics instigated by capillarity. This understanding mirrors the description of bubble pinch-off in turbulence given in Ruth et al. (2019), in which we showed that turbulence sets an "initial" deformed bubble shape before the collapse dynamics overtake the turbulent dynamics. Once the inertial collapse of the neck becomes fast enough (equivalently, once the neck becomes small enough), however, the turbulence effectively "freezes" in place relative to the accelerating collapse dynamics. The end result is that the final stage of the pinching processin this case, the production of small bubbles through the capillary instability of gas ligaments-is affected by the turbulence only insofar as the turbulence sets the "initial condition" on which the remainder of the process evolves under capillary and, eventually, inertial, dynamics.

Individual break-up event dynamics
So far, we have considered the transient size distributions N (d/d H ) that result from air cavities with d 0 > d H being introduced to turbulence. In this section, we focus on individual break-up events that are tracked individually in three dimensions as described in Section 2.3; these events are the building blocks for the disintegration of larger cavities. We will characterize the break-up events over their typical time scale, which is given by the eddy turn-over time at the parent bubble's scale, T turb (d 0 ) = −1/3 d 4.1. Qualitative discussion of the break-up sequences One break-up producing m = 2 child bubbles is shown in Figure 11, and one producing m = 4 bubbles is shown in Figure 12. In each, images throughout the break-up sequence are shown in (a-c), and the three-dimensional trajectories taken by the bubbles involved are shown in (d). At each point, the bubble's size is computed relative to the local Hinze scale, and d/d H is encoded in the trajectory color. The spatial scale is given in terms of the integral length scale at the breakup location, L int,0 , showing that the bubble trajectories are resolved over multiple integral length scales. Panel (e) shows the dimensional diameters of the bubbles involved over time.
In the case of binary break-up, given in Figure 11, the parent bubble enters the imaged volume from the foreground, and quickly encounters a region of more intense turbulence, where d 0 /d H increases. Eventually, having become deformed, the bubble pinches apart into two child bubbles, each of which are comparable in size to the parent. Both child bubbles persist in the field of view for at least a tenth of a second (∼ an integral-scale turn-over time) without breaking.
In the more complex break-up shown in Figure 12, the parent bubble similarly traverses from a region of less intense turbulence to more intense turbulence, increasing the value of d 0 /d H . Eventually, at t = 0.170 s (shown in panel (a)), the bubble becomes elongated in the vertical direction, and in a sequence of two rapid splitting events produces the three child bubbles that are visible at t = 0.192 s The "family tree" for the single break-up, giving diameters of the bubbles present at each point in time. Fainter lines give the instantaneously-measured diameters, and straight lines give the median for each bubble, which is the quantity we consider in our analysis.
(shown in panel (b)). One is still larger than d H , one is of the order of d H , and the third, left between the two, is smaller than d H . The bubble of the order of the Hinze scale is still significantly deformed, the capillary dynamics involved with the break-up not yet having relaxed. A short time later, by t = 0.201 s (shown in panel (c)), an additional bubble has split from it, leaving a total of four child bubbles.

Identification of break-up events
We identify bubble break-ups like the ones shown in Figures 11 and 12 as being sequences of bubble splitting events not exceeding the eddy turn-over time at the parent bubble scale, T turb (d 0 ) = −1/3 d 2/3 0 . To enforce this temporal constraint, we first construct a "family tree" of all splitting events recorded in one run. Then, if any bubble is present at a time T turb (d 0 ) beyond the initial detected break-up of the first bubble (with diameter d 0 ), we truncate the family tree at that bubble, and start a new family tree with the same bubble (if it later breaks apart). After doing so, we store the sizes of the parent bubble and child bubbles, as well as the turbulence characteristics spatially interpolated at the initial breakup location. To remove spurious break-ups, we discard those for which the sum of The "family tree" for the single break-up, giving diameters of the bubbles present at each point in time. Fainter lines give the instantaneously-measured diameters, and straight lines give the median for each bubble, which is the quantity we consider in our analysis.
the calculated volumes of the m child bubbles is less than 50% of, or more than 200% of, the calculated volume of the parent bubble. In total, we captured 162 bubble break-ups with this dynamical tracking approach that fit the volume conservation criteria. Figure 13 (a) shows the distributions of the break-up conditions (the Hinze scale at the break-up location and the parent bubble size) for the aggregated dataset, which we later break down by the parent bubble's size relative to the Hinze scale. The parent bubble diameter d 0 is typically slightly larger than the Hinze scale, as the distribution of d 0 (the green line) is located just to the right of that of d H (the dashed red line). Thus, the break-ups we capture in this experiment have d 0 /d H ≈ 0.4-3.7. The black curve shows the distribution of the sizes of child bubbles formed during break-ups, integrating to the average number of bubbles formed per break-up event.
To gauge the effect of inhomogeneity in the generated turbulence, we consider how the local turbulence intensity experienced by the bubble (in a Lagrangian sense) varies over timescales relevant to the bubble's break-up. Ideally, a bubble would not be advected through statistically inhomogeneous turbulence during the course of its break-up. Denoting the Hinze scale at the bubble's location at time t as d H (t), Figure 13 (b) shows the Hinze scale at the break-up location scaling, which is approached for large d 0 /d H due to the production of small bubbles by capillary instabilities (Rivière et al. 2022). Qualitatively, the child size distribution for smaller parent bubbles is flatter near the Hinze scale, while that for larger parent bubbles increases more rapidly with decreasing bubble size as a power-law relationship. Note that the child size distribution is defined so that it integrates to the average number of child bubbles formed. This representation of the child size distribution masks the large number of bubbles formed very close to the parent bubble size. To capture these small bubbles, we also compute the volumetric child size distribution, normalized by the volume of the parent bubble V 0 . Since the determination of the volumes of larger bubbles is difficult given their deformations, we approximate the parent bubble volume as the sum of the volumes of the child bubbles, and consider (Vejrazka et al. 2018). The distribution of these dimensionless volumes is shown in Figure 14 (b), exhibiting a ∪ shape that is not strongly dependent on d 0 /d H (though we again see increased small bubble production with larger d 0 /d H ). The large values of this distribution near 1 suggest that in many break-up events, small bubbles are "torn off" of the parent bubble, without inertial deformation producing multiple child bubbles of sizes comparable to that of the parent. We note that the resolution of our experiment (in which the smallest bubble we can detect is approximately 200 µm in diameter) limits the number of bubbles detected.

Small bubble production without significant inertial deformation
In many of the break-ups we observe in the large cavity disintegration and individual break-up experiments, small bubbles were seen to be "torn off" from a parent bubble, without an appreciable large-scale deformation to the parent bubble. These events are reminiscent of tip-streaming (Montanero & Gañán-Calvo 2020). This phenomenon is evidenced by the right side of the ∪-shaped child size distributions shown in Figure 14 (b), as a child bubble that is nearly the size of the parent is the signature of such break-ups. To understand these events, we present in Figure 15 a qualitative discussion of the dynamics of individual splitting events. For each splitting event, we compare the velocity of the parent bubble at break-up v parent (denoted by the gray arrow in panel (a)) to the displacement between the parent bubble's final position x parent (the gray circle) and the initial positions at which the child bubbles are detected x child (the black circles). The child bubble's initial detected position ahead of or behind the parent bubble, κ = v parent · (x child − x parent )/(u L int ), normalized by turbulence quantities, is then computed, and is plotted against the child bubble's size relative to the parent size in panel (b). The color of each marker denotes the size of the splitting event's parent bubble relative to the Hinze scale. The black line shows the expected value of κ given the normalized child bubble size. Smaller child bubbles (with d child /d parent < 0.6, below which the mean value of κ becomes negative) tend to be left in the wake of the parent bubble (κ < 0), while larger child bubbles tend to be produced ahead of the parent bubble (κ > 0). While the conceptual picture for break-up discussed in Section 3.3 describes the role of capillarity during break-ups involving large-scale deformations, it is likely that break-ups solely involving small bubble production are also regulated by capillarity: in these cases, a turbulent motion smaller than the parent bubble may succeed in producing a ligament which extends off of one side of the parent, and this ligament may pinch apart into many small bubbles in a capillary instability as it is retracted back into the bulk of the parent bubble. Specifically, Figure 15 suggests that the bulk of a bubble may often be swept forward by a turbulent eddy, and the trailing ligament may become unstable as it "catches up" with the rest of the parent bubble. Similar to the framework presented in Section 3.3, the process is initiated by a turbulent deformation to the parent, and ends with the capillary instability of a ligament involved in the deformation.

A model for bubble break-up
5.1. Physical ideas The experiments presented in Sections 3 and 4, taken together with the existing literature, point to three important time scales that must be considered in developing a population balance model: the inverse of the break-up frequency, the break-up duration, and the capillary capillary pinching times.
The longest of these is the typical duration until a break-up occurs-that is, the inverse of the break-up frequency, 1/ω(d 0 ). This time scale will control how many break-up events will occur over a given time and will be a function of d 0 /d H . The second timescale is that over which a break-up typically occurs, or the event duration (i.e., lasting from the start of the deformation until the child bubbles have all been formed), and will also be a function of d 0 /d H . The breakups taking the longest time will be those instigated by the largest eddies capable of causing break-up, which are taken to be those at the parent bubble's scale (Luo & Svendsen 1996). Thus, an upper bound and typical scale of the breakup duration is taken to be the eddy turn-over time at the parent bubble's scale, T turb (d 0 ) = −1/3 d 2/3 0 , in agreement with experimental and numerical observations of the time over which bubbles are deformed prior to break-up (Risso & Fabre 1998;Martínez-Bazán et al. 1999b;Rivière et al. 2021). The final timescale we consider is that of the capillary instabilities of gas ligaments that produce a small child bubble of size d, which will occur over the capillary timescale of that child bubble, T cap (d) = (ρ/γ) 1/2 d 3/2 /(2 √ 3) (Rivière et al. 2022). From these three relevant time scales, we define three types of events. At the shortest time, we define the individual binary splitting events. For the production of bubbles with d d H , we have T cap (d) T turb (d 0 ). At the eddy turn-over time, we define a break-up as being a sequence composed of all the splitting events occurring in a time bounded by ∆T break−up = T turb (d 0 ), which permits the production of more than two bubbles in a single event (similar to the definition used for drop break-ups by Solsvik et al. (2016)). Finally, following the nomenclature from Hinze (1955), a disintegration is a longer-duration process involving an arbitrary number of break-ups.
These timescales are sketched in Figure 16, which illustrates two break-up events that stem from a bubble of diameter d A encountering turbulence. The deformation to the parent bubble that instigates the break-up is assumed to happen within a time T turb (d A ) before the first bubble splits from the parent. Then, within an additional time bounded by T turb (d A ), subsequent splitting events occur due to capillary instabilities arising from the deformation. One such instability produces a bubble with diameter d C , and the time over which this instability develops is set by the capillary timescale at the smaller child bubble size, T cap (d C ). Later on, one of the child bubbles produced in the first break-up, with diameter d B , itself breaks up.
Using these ideas, we propose a population balance model that integrates these physical elements and models the evolution of a bubble size distribution with a Boltzmann transport equation using the bubble size as an internal coordinate. The population balance model considers a break-up rate kernel f , constructed from child size distributions computed through a Monte Carlo approach (constrained by results from experiments and DNSs, informing the number of children and the shape of the distribution) and a parent bubble break-up frequency taken from the literature. With the kernel defined, we integrate the model in time to simulate the evolution of the size distribution during a cavity disintegration and compare to our experimental data.

Population balance modeling
In a confined region of homogeneous turbulence, the transient evolution of the absolute dimensionless volumetric bubble size distribution is the absolute dimensional volumetric size distribution,Ṽ = V /V H , andt = t/T int , is described by  Figure 16: Sketch of two bubble break-ups and the associated timescales, with T turb (d) = −1/3 d 2/3 and Tcap(d) = (ρ/γ) 1/2 d 3/2 /(2 √ 3). The gray vertical lines denote the times associated with each of the two break-ups. The shaded region to the left bounds the time over which the deformation to the parent bubble is assumed to occur (the turbulent timescale at the parent bubble size), and the region to the right of the line bounds the time over which the subsequent splitting events are assumed to occur (also taken to be the same turbulent timescale). During the subsequent splitting events, the capillary timescale at the size of the smaller child bubble sets the time over which the splitting event occurs (Rivière et al. 2022). The time between break-ups is set by the inverse of the break-up frequency ω of the bubble which is to break, which we address later in the paper.
where the first term on the RHS gives the rate of consumption of bubbles of volumeṼ due their break-ups, and the second term on the RHS gives the rate of production of bubbles of volumeṼ due to the break-ups of larger bubbles (Martínez-Bazán et al. 2010). The break-up kernelf (δ;∆) = f (δ; ∆)V H T int can be decomposed into a parent break-up frequency and volumetric child size distribution withf (δ;∆) =ω(∆)p(δ;∆), with the dimensionless break-up frequencyω(∆) = ω(d)T int and dimensionless volumetric child size distributioñ p(δ;∆) = p(δ; ∆)V H . Thus, we can moveω(∆) outside the integral in the first term on the RHS and invoke Ṽ 0p (δ;Ṽ )dδ = m (Ṽ ), with m (Ṽ ) the average number of bubbles formed in the break-up of a bubble of volumeṼ , to express the bubble consumption term as simply −N V (Ṽ ,t)ω(Ṽ ). Note that we definep(δ;∆) so that it integrates overδ to the average number of child bubbles formed by the break-up of a bubble of volume∆.

Construction of the child size distributions
We develop a parameterization of the break-up volumetric child size distributioñ p(δ;∆) that accounts both for child bubbles produced by both the slower inertial mechanism (occurring over the eddy turnover time) and the faster capillary pinching mechanism (occurring over the capillary timescale of the small child bubbles) using a Monte Carlo approach. We consider a set of rules constrained by our experimental and numerical observations describing the outcomes of individual break-up events, then aggregate the outcomes of these events into child size distributions.

Statistics on the number of child bubbles formed
A key step in modeling each break-up is to constrain the distribution of the number of bubbles formed in each event. To this end, we first consider the data from our dynamical experiments given in section 4. The average number of child bubbles larger than the experimentally-resolvable minimum size d min /d H ≈ 0.07, m , is shown in Figure 17 (a). As the parent bubble increases in size, more child bubbles are typically produced. Given the steep dependence ofp(δ;∆) oñ δ, we must qualify each observation of m with the minimum resolved bubble size to better enable comparisons between different experiments. For compactness, however, we take all m values to be the number of resolved bubbles larger than 0.07d H unless otherwise noted. Our experimental observations of m , binned by d 0 /d H , are shown in the black squares, and the gray region around them bounds ± one half of a standard deviation around the mean. Next, to consider the number of bubbles produced in the break-ups of larger bubbles, we turn to data from the disintegration of the three largest cavities presented in section 3. With the assumption that the initial splitting event happens nearly instantly after the bubble is released into the turbulence, to apply the same definition of the duration of the break-up, we define m for this dataset as the number of resolved bubbles present after one eddy turnover time T turb (d 0 ) = −1/3 d 2/3 0 has elapsed after the cavity release, which are denoted by the open circles in fig. 8 (a) and fig. 17 (a). We invoke eq. (C 1) to apply a slight adjustment to these numbers in order to extrapolate results to the finer spatial resolution of the tracked break-up experiment, as discussed in Appendix C. The number of bubbles in the extrapolated range constitutes about 30% of the ones in the observable range. These adjusted values are shown as the filled-in light blue circles in Figure 17 (a).
We have additionally re-analyzed the DNSs of bubbles breaking in homogeneous, isotropic turbulence presented in Rivière et al. (2021Rivière et al. ( , 2022, tracking the bubble break-up events in a similar way to what has been done on the experimental data in Section 4. From these DNSs, we can compute the average number of bubbles formed per event as a function of the parent bubble size, included in panel (a) as the red star markers. Open stars give the original observations, for which d min /d H = 0.25, while the filled-in stars give the number adjusted for the spatial resolution. Note that while we consider We c = 1 for the experimental data, the value of d H for the DNS is given by We c = 3 .
Finally, as a comparison, the open gray markers show the (un-adjusted) number of bubbles detected experimentally in break-ups by Vejrazka et al. (2018), in which break-ups varied in and d 0 (which is denoted by the marker style). As shown in their paper, once collapsed to d 0 /d H , the dependence on the dimensional bubble size nearly disappears.
The four datasets (our two experiments, those from Vejrazka et al. (2018), and DNSs from Rivière et al. (2021)) produce a coherent picture regarding the number of bubbles formed. When d 0 /d H is small, break-ups tend to be binary, producing on average 2 child bubbles after T turb (d 0 ). As d 0 /d H increases, the number of child bubbles increases. Surface tension is less effective at preventing the severe deformation of larger bubbles, leading to more complex deformed bubble shapes that yield a greater number of child bubbles. The orange curve in panel (a) shows a fit to the data of the form where m min = 2 and the fit constants are b 1 = 4 and b 2 = 2.3. Figure 17  with m + m min the mean number of children, a function of the parent bubble size.

A stochastic model for each break-up
The Monte Carlo approach involves running many iterations of a stochastic model and developing a statistical representation of the aggregated results. Each discrete simulation of a break-up mirrors the physical processes involved: the bubble, sketched in Figure 18 (a), is first deformed into two lobes, shown in panel (b), and then some number of capillary bubbles are created as the neck separating the lobes collapses to create the two inertial child bubbles. For each iteration (i.e., one simulated breakup) at a given value of∆, we first define the number of bubbles m that will be produced by picking a value of m from the distribution r(m ; d 0 /d H ) given by eq. (5.3), adding m min = 2, and rounding to the nearest integer. We pickδ min = 0.07 3 in order to match the experimental dataset on which the parameterization of m is based. As we will show, once the p.d.f.s have been constructed for this given value ofδ min , it will be straightforward to extend them to lower or higher values ofδ min .
For cases in which m 3, the capillary mechanism produces m = m − 2 bubbles, whose sizes follow a ∝δ α distribution with α = −7/6 (corresponding to the P d (d/d H ; d 0 /d H ) ∝ (d/d H ) −3/2 scaling described by Rivière et al. (2022), as distributions in diameter are related to those in volume by -Bazán et al. 2010;Qi et al. 2020)). As is sketched in Figure 18 (c), the volumeδ cap,i of capillary bubble i is picked from a power-law distribution with slope α, bounded betweenδ min and the maximum allowable volume for a capillary bubble given the previously-produced bubbles,δ cap,max,i . For the production of the first capillary bubble, we setδ cap,max,1 =∆ (noting that the steep slope of P d (d/d H ; d 0 /d H ) with respect to d/d H makes the production of capillary bubbles this large uncommon). For the production of the remaining capillary bubbles, we setδ cap,max,i =∆ − i−1 j=1δ cap,j . At each step of the process, ifδ cap,i is greater thanδ cap,max,i /2, we replace it withδ cap,max,i /2 −δ cap,i , such that for any splitting event, the smaller of the two produced does not further split.
Once the volumes of the m capillary bubbles are specified, we must determine the volumes of the two inertial bubbles. To that end, we first compute the portion of the parent bubble volume that has gone to the capillary bubbles, χ cap = m i=1δ cap,i /∆. The size of the first of the two inertial child bubblesδ inertial,1 is drawn uniformly from the remaining bubble volume, (1 − χ cap )∆, and the second is taken as its complement,δ inertial,2 = (1 − χ cap )∆ −δ inertial,1 . Once this is done, the volumes of all child bubbles produced in this single break-up have been determined.

Aggregation of simulated break-ups into child size distributions
For a given value of d 0 /d H (or the equivalent normalized volume∆ = (d 0 /d H ) 3 ), the process of simulating one break-up stochastically is repeated n MC = 10 5 times.
For each∆, the sizes of the bubbles produced in each of the n MC events are aggregated, and the distribution of all these child bubbles defines the volumetric child size distributionp(δ;∆). The distribution is normalized such that ∆ 0p (δ;∆) = m (∆), with m (∆) the average number of bubbles formed. Since the size distribution is aggregated from geometrically-plausible break-ups, it itself must satisfy any constraints relating to the sizes of the bubbles produced. Figure 19 (a) shows the volumetric child size distributions for five values of∆. When∆ is small, the child size distribution is nearly uniform, as the capillary production mechanism is negligible for small bubbles; for moderate∆, the child size distribution exhibits ap ∝δ α scaling for small bubbles, while remaining close to flat for bubbles near the parent bubble size. For even larger bubbles, for which the capillary production mechanism is the most effective, the entire distribution approaches aδ α scaling. For each∆, we also obtain χ cap (∆), shown in Figure 19 (b), by averaging the portion of the parent bubble volume going to the capillary child bubbles χ cap over the n MC events. When∆ 1, χ cap ≈ 0, and essentially all of the parent bubble volume goes to the two inertial child bubbles. With larger∆, χ cap increases, reaching χ cap = 0.1 at∆ = 60. Even at∆ = 1000, less than half of the parent bubble volume goes to the capillary bubbles.
We then fit each volumetric child size distribution as a sum of two components, each stemming from one of the two mechanisms of child bubble production, with α = −7/6 set by the distribution from which the capillary bubbles are picked and γ(∆) chosen to match the aggregated Monte Carlo simulation data. The two remaining coefficients, a(∆) and b(∆), are constrained by the volume going to bubbles produced by each mechanism, leading to (5.6) The fits to each child size distribution with eq. (5.4) are shown as the faint, thick lines in Figure 19 (a). Figure 19 (c) shows the evolution of the exponent γ(∆) describing the inertial production mechanism. Values of χ cap (∆) and γ(∆), which together contain all the necessary information about the child size distributions, are stored for many values of∆. To implement the child size distributions in a population balance model, we interpolate χ cap (∆) and γ(∆) for a given value of∆. Data for each curve and Python code to construct the child size distributions will be provided online at publication for those wishing to implement the model we have constructed.

Parameterization of the break-up frequency
The next step is to parameterize how often the break-ups will occur. Here, using an approach that has been successfully applied to the break-up of oil droplets in turbulent jets (Aiyer et al. 2019;Aiyer & Meneveau 2020), we integrate the effects of eddies smaller than the parent bubble size (each of dimensional diameter d e ) which contribute to break-up (Prince & Blanch 1990;Tsouris & Tavlarides 1994), (5.7) where K is an order-1 constant we adjust,ũ turb (d e /d H ) = C 1/2 2 1/3 d 1/3 e /u = C 1/3 (d e /d H ) 1/3 (d H /L int ) 1/3 is the dimensionless turbulent velocity scale of the eddy, (d e /L int ) −4 is the approximate dimensionless eddy density (Solsvik et al. 2016), and Ω(d e /d H ; d 0 /d H ) is the break-up efficiency given the eddy and bubble sizes. Neglecting viscous effects (given the low viscosity of air bubbles), the breakup efficiency, which gives the probability that an eddy has sufficient energy to overcome surface tension, is taken as the inverse of the exponential of the ratio between the average change in surface energy associated with the break-up E σ (d 0 ) and the kinetic energy of the eddy E eddy (d e ), exp(−E eddy (d e )/E σ (d 0 ). The average surface energy change is given dimensionally by (5.8) with the proportional change in surface area due to break-up Γ dependent on the form of the child size distribution according to The kinetic energy of the eddy is given by E eddy (d e ) = (π/4)ρd 3 e C 2 ( d e ) 2/3 . Expressed in our non-dimensional units, the break-up efficiency is then with the critical Weber number We c necessary to link the scales of the bubble and the turbulence. With each component specified, eq. (5.7) is evaluated numerically and is shown in Figure 20, using K = 2 picked through a comparison to the experimental data given in Section 3. The break-up rate increases as bubbles approach the Hinze scale and then plateaus due to two competing effects: while larger bubbles are susceptible to a wider range of turbulent scales that may cause break-up, they tend to break into many more bubbles than smaller ones do, leading to a greater surface energy term in eq. (5.7). This means that while more eddies are interacting with the parent bubble, each is less likely to have sufficient energy to cause a break-up.
The thicker gray line in Figure 20 gives the inverse of the turbulent turnover time at the parent bubble scale, which we take to set the duration of each break-up event. The break-up frequency is thus consistent with the breakup duration, sinceω(∆) = gT int being strictly less than T int /T turb (d 0 ) means that the typical duration of a break-up is never longer than the typical time between such break-ups. Finally, the dotted orange line gives the inverse of the (dimensionless) capillary timescale at the parent bubble scale, T int /T cap (d 0 ), showing that capillary effects happen faster than both the break-up duration and time between break-ups (up until the largest bubbles we consider). The capillary pinching events responsible for sub-Hinze bubble creation thus occur over even shorter durations, as the capillary timescales of the small child bubbles formed will be much faster than that of the parent bubble.
5.5. Summary of parameters involved in the model To summarize, Table 2 lists each parameter in the model and explains how each is determined.

Model comparison to transient air cavity disintegration data
Withp(δ;∆) andω(∆) now fully specifyingf (δ;∆), we can simulate the turbulent disintegration of cavities we studied experimentally in Section 3 by picking the appropriate initial condition for each (i.e., N (d/d H ) giving one bubble of size d 0 /d H ) and integrating eq. (5.1) in time. Figure Figure 19 (a)). With larger parent cavities, the sub-Hinze distribution steepens as the capillary mechanism contributes more significantly to the child size distributions for parent bubbles of these larger cavity sizes (as evidenced in the χ cap (∆) curve shown in Figure 19 (b)).   larger-scale turbulent deformations as the mechanism responsible for small bubble production. Crucially, significant small bubble production by this mechanism is limited to parent bubbles with d 0 d H , as only bubbles much larger than the Hinze scale can become deformed to a severe enough extent to produce the ligaments from which the small bubbles originate.

Model element Equation
The first piece of evidence we provide for this role of capillarity is visual: Figures 9 and 10 show a number of instances of small bubbles being left behind after the collapse of gas ligaments. Second, the experimental N (d) ∝ d −3/2 scaling for d < d H with d 0 d H is coherent with the P (d) ∝ d −3/2 scaling for the breakup child size distribution reported by Rivière et al. (2022), who showed that the lifetime of ligaments before their collapse to produce a bubble of size d coincides with the capillary time scale of a bubble of size d, T cap ∝ d −3/2 .
We implemented these physical ideas in a population balance model of turbulent bubble break-up. The child size distributions describing individual break-up events were constructed with a Monte Carlo approach involving simulations of many break-ups. The statistics of each simulated break-up are prescribed by our understanding of the role of capillarity and additional experimental results on individual bubble break-up in which parent and child bubbles were tracked dynamically in three dimensions. The resulting expression for the child size distribution, eq. (5.4), involves two components: one describes the effect of the large-scale deformation to a parent bubble by an energetic turbulent eddy, and the other describes the action of capillarity in producing small bubbles. Finally, the rate at which parent bubbles undergo break-ups was determined by integrating the action of eddies below the bubble's size, which all contribute to break-up. The complete model (consisting of the child size distributions and the parent break-up frequency) yields a good match to our transient experimental data.
Along with the recent analysis of DNSs of bubble break-up in turbulence from Rivière et al. (2021Rivière et al. ( , 2022, this experimental work opens the door to a new understanding role of capillarity in turbulent bubble break-up, in which surface tension not only counteracts the initial turbulent deformation to a bubble but also leads to the formation of sub-Hinze bubbles through capillary instabilities that arise during the final stages of the break-up process. Bubble sizes are detected with image processing of the images of the cavity disintegrations in multiple stages. First, a simple image intensity threshold is applied to binarize each greyscale image, and the bright spots at the interior of each bubble image is filled in. A first pass at extracting the bubble diameter d intensity based on this intensity threshold is then made by computing the equivalent diameter of a circle with the same projected area as the binarized bubble image. As we find that the determined sizes of small bubbles (d intensity < d cutoff , with d cutoff = 1.5 mm) are sensitive to the image intensity threshold chosen, we individually employ a Canny filter (Canny 1986) to the images of each of these small bubbles to find their borders. Their diameter d is then defined as the equivalent diameter of the projected area inside the bubble border. For larger bubbles (with d intensity > d cutoff , for which the Canny edge detection often fails due to the deformed bubble shape), we define the diameter as d = d intensity + σ d , where σ d = −25 µm is the typical value to which d Canny − d intensity asymptotes for bubbles approaching d cutoff .
A.2. Adjusting the size distribution to account for bubble advection and smoothing of the size distributions The time-dependent size distribution N (d/d H , t/T int ) is first computed on a frame-by-frame basis, and then is averaged in time over a window with width τ /T int that increases with t/T int . The width τ is first set as τ /T int = 0.13 + 0.38(t/T int ) and then clipped at τ /T int = 2.67, which it reaches at t/T int = 6.68.
Due to the buoyant rise of the bubbles and their advection by the turbulence in the air cavity disintegration experiments, bubbles leave the measurement region over time. This complicates the analysis of our data: if we were to solely consider bubbles viewed in-frame, we would calculate a rapid loss of bubble volume as bubbles leave the field of view. Results would be further skewed by any size dependence of the bubbles' motions. Indeed, the transient bubble size distributions based on the bubbles viewed in-frame shown in the left of Figure  The adjusted size distributions with the slight smoothing applied, which we consider in the paper. (d) The total number of bubbles detected in-frame (gray) and the total number of bubbles, including the advection adjustment (black). (e) The sum of the volumes of the bubbles detected over time, normalized by the known volume of the air cavity. The gray curve is the volume of bubbles detected in-frame; the black curve is that curve added to the volume of bubbles from the advection adjustment, in green.
the total volume and number of bubbles tracked, shown in the fourth and fifth columns respectively as the gray lines, decreases as bubbles leave the field of view. We address this experimental limitation by tracking the bubbles' motion in two dimensions and making note of when bubbles leave or enter the measurement region near one of its four borders. Then, to compute N (d) for some time t, along with the bubbles of size d imaged at time t, we add counts of all the bubbles of size d that have left the measurement region before time t, and subtract counts of all the bubbles of size d that have entered the measurement region before time t. The resulting adjusted size distributions are shown in the second column of Figure 22.
The final step in the processing is to slightly smooth each N (d/d H , t) curve in bubble size by an amount dependent on the total number of bubbles present at each time, n(t) = ∞ 0 N d(d/d H ). Size distributions are computed with 30 geometrically-spaced bins between d/d H = 0.14 and 17.5. Then, at each time, each curve is smoothed with a Gaussian filter with a standard deviation of σ smooth (t) bins, with σ smooth (t) picked given an empirical function of the number of bubbles present. When fewer than three bubbles are present, σ smooth (t) is set to one bin; when more than thirty are present, σ smooth (t) is 0.1 bins. σ smooth (t) is interpolated between these two limits when a moderate number of bubbles are present. This approach ensures minimal smoothing when a sufficient number of bubbles are present and a moderate amount of smoothing when few bubbles are present. The smoothed size distributions, which we consider in the paper, are shown in the third column of Figure 22.
This advection adjustment approach effectively "freezes" in place the record of bubbles as they leave the measurement volume. The green regions in the third and fourth columns of Figure 22 show the additional number of bubbles and corresponding additional bubble volume added to the bubble record with this method. The black lines show the sum of the in-frame measurements and this adjustment. Any limit obtained by the adjusted n curve is still not especially physically meaningful, as some of the bubbles which have exited the field of view are not much smaller than the Hinze scale, so they would eventually break apart further if left within the turbulence region. However, it is qualitatively closer representation of the "true" behavior than would be obtained through just the in-frame measurements.
The plots of the summed bubble volume (normalized by the cavity volume) shown in the fourth column of Figure 22 reveal a second limitation in our bubble detection method: the total volume of bubbles considered, ∞ 0 V N V (V )dV , is not a constant value equal to the known volume of the air cavity V 0 . At early times, the total bubble volume is under-counted, and becomes over-counted at later times. This is due to the highly-deformed shapes of large bubbles, for which the equivalent diameter determination we employ is only a rough approximation. Further, bubbles whose images overlap can be detected as a single larger bubble.
We note that, aside from some representative images taken at late times, the latest measurement of N (d/d H ) presented in the paper or employed in our analysis is t/T int = 5, at which point the advection adjustment has had only a moderate impact on the size distribution for all cavity sizes.
Appendix B. Comparison of cavity disintegration with and without turbulence The cavity release experiment detailed in Section 2.2 was also run with the turbulence-generating pumps turned off, such that the cavity was released into otherwise still water. In these experiments, the extent of the bubble production is greatly reduced. Figure 23 (a) and (b) shows snapshots of two bubbles of the same size, at the same time after their release, into quiescence and turbulence, respectively. The bubble released into quiescence is deformed by buoyancy but has not broken apart; the bubble released by turbulence has undergone break-ups. Bubbles released into quiescence occaisionally break due to the cup's motion or their large buoyant deformations (Landel et al. 2008), but Figure 23 (c-e), which show the number of bubbles produced over time in turbulence and quiescence for the three largest cavity sizes, indicate that the moderate bubble production from these break-ups is negligible compared to the much greater production in turbulence.
Appendix C. Adjustment of the number of resolved bubbles to account for the minimum resolved child size Since we find a bubble size distribution that scales as N (d) ∝ d α d with α d < −1 for small bubbles, the total number of bubbles above some minimum size will diverge as that minimum size decreases, up until some additional physical limit is encountered. Therefore, to enable a more direct comparison between datasets in which the experimental or numerical resolution differs, we can adjust the total number of bubbles formed in a break-up to account for the different resolved sizes.
Let us denote by m[d > zd H ] the average number of resolved bubbles larger than zd H that are formed in a break-up. Given a known value of m[d > xd H ], we can find the corresponding value of m[d > yd H ], which is the hypothetical number resolved had a minimum spatial resolution of yd H been employed. Following the conceptual model presented in section 3.3, we assume that all but two of the child bubbles produced in each break-up follow a power-law scaling ∝ (d/d H ) α d , with α d = −3/2. With this assumption, we calculate the appropriate prefactor for the sub-Hinze distribution given the observed value of m[d > xd H ], then extend the distribution to d min /d H = y and integrate over all the larger bubble sizes to get the effective number in the range that is resolvable in the hypothetical experiment, yielding