Fragmentation in turbulence by small eddies

From air-sea gas exchange, oil pollution, to bioreactors, the ubiquitous fragmentation of bubbles/drops in turbulence has been modeled by relying on the classical Kolmogorov-Hinze paradigm since the 1950s. This framework hypothesizes that bubbles/drops are broken solely by eddies of the same size, even though turbulence is well known for its wide spectrum of scales. Here, by designing an experiment that can physically and cleanly disentangle eddies of various sizes, we report the experimental evidence to challenge this hypothesis and show that bubbles are preferentially broken by the sub-bubble-scale eddies. Our work also highlights that fragmentation cannot be quantified solely by the stress criterion or the Weber number; The competition between different time scales is equally important. Instead of being elongated slowly and persistently by flows at their own scales, bubbles are fragmented in turbulence by small eddies via a burst of intense local deformation within a short time.

P erhaps no other area of fluid dynamics has borne a twin problem more than bubble breakup 1 and turbulence cascade 2 both by Andrey N. Kolmogorov, based on a key idea of elementary entities, i.e., bubbles and eddies, being fragmented into smaller and smaller sizes, following a universal mechanism. In 1955, Hinze 3 extended Kolmogorov's original idea 1 , and this Kolmogorov-Hinze (KH) framework has since posed deep and lasting impacts on modeling turbulent bubble/ drop fragmentation in various flow configurations 4-6 and applications, including emulsion 7 , spray formation 8 , and raindrop dynamics 9 .
The key hypothesis in the KH framework is that, in turbulence, bubbles/drops with diameter D are broken by eddies of the same size and the contribution from sub-bubble scale eddies is negligible. The most important dimensionless number based on D is thus the Weber number. The fundamental challenge associated with this key hypothesis is not about its correctness but its falsifiability. For fully-developed turbulence, eddies of many length scales are present at the same time. In these situations, bubbles always encounter eddies of various sizes, so it is extremely difficult to disentangle them cleanly 10 , not to mention establishing their roles in bubble breakup. Therefore, there has been no direct experimental evidence so far to either support or refute this hypothesis.
To overcome the aforementioned problem and to directly test this hypothesis, we seek a "magic knife" to cleanly separate the sub-bubble-scale eddies from the bubble-sized ones, and expose bubbles to only one type at a time. A flow configuration that features two identical vortex rings colliding symmetrically head-on was identified as an excellent candidate because it possesses two unique stages. In the early stage, both rings expand but remain intact (Fig. 1a, the red ring), and the vortical structure only contains one large scale. As it progresses, two rings become unstable, and the large-scale flow breaks down into a turbulent cloud, comprised mainly of small eddies (Fig. 1a, the blue ring). Recently, it has been shown that the generated turbulent cloud in this later stage shares similar statistics with other types of fully-developed turbulence 11,12 . The power spectrum follows the same −5/3 power law in the inertial range, and it lasts for an extended period of time before turbulence starts to decay. In this flow, a bubble tends to experience either only the bubble-sized vortices or turbulence consisting of the sub-bubble-scale eddies. This allows us to distinguish contributions from eddies with various sizes and directly examine the key hypothesis in the classical KH framework.

Results
Turbulence production. Figure 1a shows a schematic of the experimental apparatus that features a vortex collision sub-system (Fig. 1b) and a bubble injection sub-system. The dashed box indicates the measurement volume close to the bottom of the rings. Additional details can be found in Methods. Two distinct stages of the developed flows are highlighted in red and blue colors. The early stage was dominated by smooth and intact vortex rings, and the later stage was filled with many small eddies. Careful system control was designed to ensure that a bubble always rises to the same height when the two rings just touch each other. As shown in Fig. 1c, bubbles (indicated by the green blobs) that got entrained into one of the vortex rings were carried downward and experienced two different types of flows. The smooth red ring and rough blue ring illustrate the intact vortex ring at the early stage and a turbulent cloud at the late stage, respectively. b Picture of the linear actuator and piston-cylinder assembly for generating the vortex rings. c Breaking bubbles (3D green blobs) at two different stages within the view volume. d The longitudinal structure function D LL ; The blue and green shaded areas represent the range of the bubble size and vortex size, respectively. e The averaged ratio between the z-component vorticity ω z and the total vorticity magnitude ω in the vortex ring as a function of time t. The two dashed lines indicate the two extreme limits from the uni-directional vorticity in the early stage to the fully-isotropic case later.
To quantify the statistics of the flow structures, the longitudinal second-order structure functions D LL (for details, see Supplementary Information) at different times are shown in Fig. 1d. D LL indicates the turbulence kinetic energy distributed across different length scales. At early times (t < 0.1 s) after the ring collision, D LL exhibits a clear peak at around 10-15 mm, close to the vortex size. Most of the kinetic energy is kept within this scale. This peak decays over time as the two vortex rings become unstable and break. As the kinetic energy cascades down to smaller scales, D LL at these scales rise up until t = 0.12 s when D LL reaches a 2/3 scaling law. This inertial-range scaling law is well known in the fully-developed turbulence based on the Kolmogorov theory 2 , i.e. D LL (r) = C 2 (〈ϵ〉r) 2/3 , which is indicated by solid lines in Fig. 1d. The prefactor C 2 ≈ 2 is the Kolmogorov constant, and 〈ϵ〉 is the mean turbulence energy dissipation rate. From t = 0.12 s to 0.18 s, although the kinetic energy at the vortex size continues to drop rapidly, it does not decay as much for the inertial-range scales where D LL exhibits the 2/3 scaling. It implies that the kinetic energy generated via the breakdown of large vortices compensates the energy dissipated at the small scales, leaving a much slower energy decay for the intermediate inertial-range scales. This inertial-range scaling lasts for roughly one integral timescale (56 ms), which is much longer than the typical bubble breakup time that will be introduced later (see Fig. 2a bottom). This 2/3 scaling law also covers the range of bubble sizes used in the current experiments, marked by the blue shaded area. From this scaling law, the range of turbulence energy dissipation rate 〈ϵ〉 can be estimated to be roughly 0.20-0.25 m 2 /s 3 .
One may expect that, as the vortex rings break down to a turbulent cloud, the flow should become more isotropic. To quantify the flow isotropy, the ratio between the z-component vorticity ω z and the total vorticity magnitude ω (Supplementary Information) is shown in Fig. 1e. Two dashed lines mark the two limits of 〈ω z /ω〉: 〈ω z /ω〉 = 1 if the original vortex rings remain intact and hω z =ωi ¼ 1= ffiffi ffi 3 p if the flow becomes fully isotropic. In Fig. 1e, 〈ω z /ω〉 drops gradually with time, indicating that the flow indeed approaches the isotropic turbulence as the cascade process continues.
Bubble breakup modes. Once the bubble was entrained into one of the vortices at the collision point, it was carried downwards by the flow. During this process, we observed two distinct bubble breakup modes, the examples of which are shown in Fig. 2a. For the first case, a bubble was deformed consistently along the z-axis until the moment of breakup. This process is relatively slow, and the bubble's interface seems to be smooth throughout the entire process, similar to what was observed in the linear-extensional flows 13,14 . For the rest of the paper, this type of breakup is referred to as the primary breakup as it occurs first and always before the moment when the two vortex rings break down to a turbulent cloud.
After the primary breakup, based on the KH framework, the daughter bubbles should become harder to break because their sizes are smaller and the bubble-scale eddies have weakened, yet it is surprising to find that the daughter bubble experiences a more violent breakup, as shown in the second case of Fig. 2a. This more violent breakup is referred to as the secondary breakup hereafter. The secondary breakups have three features: (i) a rough bubble interface with large local curvatures; (ii) complicated deformation along non-persistent directions; and (iii) short breakup time. The secondary breakup occurs within 5.1 ms, which is much smaller than 32.1 ms for the primary breakup. The two breakup modes are always correlated with the bubble breakup locations. In practice, a critical height at y c = −51 mm (corresponding to the vortex ring bottom location at t = 0.10 s after their collision) was used to separate the two breakup modes (primary y > y c ; secondary y < y c ). More discussions of this separation criterion can be found in Supplementary Information.
To quantitatively compare the two breakup modes, several key statistics of the bubble geometry, orientation, and breakup time, obtained from the 3D shape reconstruction, are provided. In Fig. 2b, the probability density functions (PDFs) of the bubble aspect ratio α, obtained from the 3D reconstructed bubble geometries from 6 ms before to the moment of breakup, for both breakup modes are illustrated. It is evident that the primary breakups typically feature a larger α compared with the secondary breakups. Furthermore, Fig. 2c shows the PDF of the bubble NATURE COMMUNICATIONS | https://doi.org/10.1038/s41467-022-28092-3 ARTICLE orientation, indicated by the angle between the bubble semimajor axis and the z-axis (θ), suggesting that bubbles have preferential alignment with the z-axis during the primary breakup, while the distribution of θ for the secondary breakup is wider due to the disturbances from the surrounding turbulence. The third statistics that can be used to distinguish the two breakup modes is the breakup timescale t c , which is defined as the time delay between the start time to the breakup instant. Note that the start time is not chosen immediately after the previous breakup, but at the minimum bubble aspect ratio closest to the breakup moment, when the bubble begins to be deformed by an eddy that will eventually break it. Figure 2d shows the PDF of t c for the two breakup modes. The secondary breakup skews significantly more towards a smaller t c compared with the primary breakup. These three statistical quantities show a consistent picture as the two examples in Fig. 2a.
The difference in the two breakup modes can be linked to the distinct breakup mechanisms involved. For the primary breakup, the large-scale vortex entrains the bubble towards its center-a local pressure minimum. For a bubble with a size close to the vortex diameter, as it reaches the vortex center, it experiences a pressure gradient that tends to compress it along the radial direction and extend it along the z-axis. Secondary breakup is more irregular, driven by a turbulent cloud filled with sub-bubble scale eddies. Bubbles break without significant elongation because the process is disrupted and accelerated by these eddies with smaller timescales, which also explains the observed difference in breakup time t c .
We emphasize that the primary breakup follows the key hypothesis made in the classical KH framework, in which a bubble is assumed to be broken by a clean and isolated vortex filament with a size close to the bubble diameter. However, most bubble breakups observed in fully developed turbulence are closer to the secondary case, where the contribution from a cloud of smaller eddies cannot be ignored.
Bubble breakup mechanism. To understand these two breakup modes, we measured the 3D flow in the vicinity of a breaking bubble along with its 3D geometry to show two examples of the primary and secondary breakups in Fig. 3a. In particular, the 3D isosurfaces of the vorticity magnitude at t = 0.06 s for the primary breakup and t = 0.14 s for the secondary breakup are used to illustrate the flow structures. For the primary breakup, the bubble is trapped within the bottom portion of one vortex ring and the flow surrounding the bubble is smooth, whereas the secondary breakup takes place in a chaotic turbulent cloud, indicated by the rough isosurface of vorticity.
Based on the 3D flows, the first quantity that can be acquired is the Weber number. At the bubble scale, the flow can be decomposed into the straining and rotational components. In practice, the velocity gradient coarse-grained at a scale close to the bubble size, i.e., e A ij , is extracted based on the velocity of multiple tracer particles around the bubble (for details, see Supplementary  Information). From e A ij , the coarse-grained strain rate tensor e S ij and rotation tensor e Ω ij can be determined: , which lead to the definition of two Weber numbers: where λ 3 (the largest compression rate) is the smallest eigenvalue of e S ij , and ω is the vorticity magnitude. The new definition of the two Weber numbers extends the original one-dimensional version in the KH framework to emphasize the contributions from the 3D straining and rotational flows. Nevertheless, the key assumption in the KH framework that the only relevant length scale is the bubble size is still applied here. Figure 3b, c show the PDFs of We S and We Ω over 6 ms duration before the moment of breakup for both the primary and secondary modes. For the primary breakup, We Ω appears to be systematically larger than We S because the bubble is compressed more by the radial pressure gradient due to the flow rotation than by the straining flow. For the secondary breakup, the peaks of the PDFs of We S and We Ω both locate at values smaller than one, and more importantly, smaller than their primary breakup counterparts.
The KH framework implies that bubbles with larger Weber numbers tend to break more easily. If it were right, we should expect a more violent primary breakup. However, the observations suggested otherwise, which clearly refute the key hypothesis in the KH framework. For the secondary breakup, although the eddy of the bubble size is much weaker, many sub-bubble-scale eddies begin to emerge. To demonstrate their appearance, we apply a high-pass rolling-average spatial filter with a filter length l = 3 mm (which is selected to be close to the bubble mean diameter) to the velocity field. The residual fluctuation velocity u < and its variance hu 2 < i only contain the contribution from small eddies (D e < l). Figure 3d shows the PDF of non-dimensionalized residue fluctuation velocity around bubbles for both primary and secondary breakups. The secondary breakup has systematically larger hu 2 < i compared with the primary counterpart, which implies that the contribution of sub-bubble scale eddies that were missed in the KH framework may be the key to understand the bubble breakup in turbulence.
Model. To extend the KH framework, we argue that both the bubble diameter and eddy size are relevant length scales. To include both, let us consider a simple scenario: a bubble with an equivalent spherical diameter of D encounters and gets deformed by a small eddy of size D e < D, as shown in Fig. 4a. Consistent with the features observed in the secondary breakup in Fig. 2a, the concave interface has a large local curvature that is presumably linked to the eddy size D e , so the interfacial stress can be determined by using σ/D e , instead of σ/D in the KH framework. Furthermore, the eddy inertia ρu 2 e has to overcome the interfacial stress to break the bubble. Therefore, a new Weber number based on the eddy velocity scale u e and length scale D e can be defined, and this Weber number (We e based on the eddy size D e ) has to be larger than one for the bubble to break, In addition to this stress constraint, another key missing piece in the KH framework is time. The eddy timescale has to be shorter than the bubble's relaxation time, otherwise, the bubble will return back to the sphere. The eddy lifetime or turnover time scales with D e /u e , whereas the bubble relaxation time can be estimated by using the natural frequency of the bubble, i.e. ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ffi 96σ=ðρD 3 Þ p =ð2πÞ (Lamb mode 2 15 ). So the second criterion can be set as D e =u e <2π ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ffi ρD 3 =ð96σÞ p . Rearranging this relationship yields another new dimensionless number (Ti, which represents the time ratio) that involves both D and D e : Eqs. (2) and (3) together provide two new dimensionless numbers, We e and Ti, which set two constraints on the eddy velocity u e , whose dependence on the eddy size D e are shown as blue and red solid lines in Fig. 4b. Based on these two constraints, the minimum eddy velocity to break a bubble (u e,min (D e , D)) can be estimated following: This equation divides Fig. 4b into two regions: breakup (u e > u e,min (D e , D)) and no breakup (u e ≤ u e,min (D e , D)). It is evident that u e,min has a non-monotonic dependence on D e /D, reaching its own minimum at the intersection of the two lines: , which results in D e /D = 0.74. To estimate the breakup probability, we ignore the bubbleinduced flow modulation and assume the surrounding turbulence follows the same statistics as their single-phase counterpart, the distribution of the eddy velocity can be calculated based on the log-normal distribution of the energy dissipation rate 16,17 follows: where σ 2 ln ϵ ¼ A þ μ ln ðL=D e Þ is the variance; A represents a large-scale variability, which is set at A = 0 for convenience; μ ≈ 0.25 is the intermittency exponent; and L is the integral length scale of turbulence. In this work, the original vortex diameter is used to estimate the integral scale, which is L ≈ 15 mm; we use 〈ϵ〉 = 0.20 m 2 /s 3 based on the previous estimation. The instantaneous eddy velocity u e can be directly related to the eddy-size based dissipation rate, following u e ¼ ffiffi ffi 2 p ðϵ e D e Þ 1=3 . As a result, the PDF of the eddy velocity, P(u e |D e ), for any given eddy size D e can be expressed as Based on this eddy velocity distribution, the probability of a bubble of size D being broken by a single eddy of size D e , p(D e , D), can be calculated by following pðD e ; DÞ ¼ R 1 u e;min Pðu e jD e Þdu e , where eddies with a velocity larger than u e,min (Eq. (4)) are integrated. Figure 4c (the red-curve group) shows p(D e , D) versus the non-dimensionalized eddy size D e /D for 〈ϵ〉 = 0.2 m 2 /s 3 . In contrast to the hypothesis made in the KH framework, a wide range of eddies smaller than the bubble size can drive bubble breakup; in fact, for the selected energy dissipation rate, eddies with D e close to 0.74D are actually much more efficient in breaking a bubble than bubble-sized eddies. Furthermore, p(D e , D) of all D e increases systematically as D grows.
So far, the discussions were limited to the breakup probability driven by a single eddy without considering that smaller eddies are more abundant. To account for this effect, the collision rate ω c (for details, see Supplementary Information) is used to define a weighted breakup probability (p*(D e , D)) following the expression of: p Ã ðD e ; DÞ ¼ ω c pðD e ; DÞ= R D 10η ω c pðD e ; DÞdD e , which is shown in Fig. 4d as red curves for different bubble sizes at 〈ϵ〉 = 0.2 m 2 /s 3 . It appears that bubbles of all sizes are still preferentially broken by eddies of size D e = 0.74D. But the contributions by even smaller eddies are growing as D increases.
In addition to the size dependence, p(D e , D) and p * (D e , D) for the same range of D but using a higher 〈ϵ〉 at 5.0 m 2 /s 3 is shown as green curves in Fig. 4c, d, respectively. In Fig. 4c, as 〈ϵ〉 increases, the probability (p(D e , D)) of one bubble being broken by one eddy of size D e increases for all D e , while maintaining its peak at D e = 0.74D, until p(D e , D) saturates for a range of D e . But once accounting for the number density difference of eddies of different sizes, the peak locations of p * (D e , D) for all bubble sizes shift to a much smaller D e , suggesting that bubbles are preferably fragmented by much smaller eddies in stronger turbulence. This result implies that the deviation from the hypothesis in the KH framework becomes more significant for either larger bubbles or stronger turbulence.
To validate the proposed model, we quantify the breakup time of the secondary breakups since these events are driven by small eddies. The expected bubble breakup time for different bubble sizes D is shown as open circles in Fig. 4e, and their respective distributions are indicated by the error bars with the upper and lower bounds marking the 30th and 70th percentiles, respectively.
In the model, if we assume that the breakup time scales with D/u e , the mean breakup time (t c (D e , D)) for a bubble of size D being broken by an eddy of size D e can thus be expressed as: Pðu e jD e Þdu e . By considering the contribution of all sub-bubble scale eddies within the inertial range (10η < D e < D), the expected bubble breakup time can finally be estimated following: where 〈. . .〉 e represents the integration over all eddy sizes. The model predicted breakup time ht c ðDÞi e is shown in Fig. 4e in comparison with the experimental results of the secondary breakups. A nice overall agreement can be observed without using any fitting parameters.
To ensure that the generalization of the model to other types of turbulent flow configuration is appropriate, in addition to our experiments, we also compare the model prediction against another dataset in fully-developed turbulence driven by a jet array 18 . In this experiment, the breakup frequency g(D, 〈ϵ〉) is measured for different bubble sizes (D) over a wide range of 〈ϵ〉. The results are shown in Fig. 4f as open symbols. The dashed line indicates the scaling predicted based on the KH framework by assuming that the breakup frequency scales with the reciprocal of the turnover time of the bubble-sized eddies, i.e. (ϵD) 1/3 /D ∝ ϵ 1/3 . It is evident that the measured breakup frequency exhibits a much steeper scaling than the 1/3 power law suggested by the KH framework. Based on our model, g(D, 〈ϵ〉) can be estimated by including the contributions from all small eddies, following gðD; hϵiÞ ¼ R D 10η ω c pðD e ; DÞdD e . The model predictions for different 〈ϵ〉 and D are shown in Fig. 4f as solid lines. It can be seen that the modeled breakup frequency agrees with the experimental results well, capturing the steeper scaling that the data suggests. This scaling originates from the growing contributions from small eddies with a larger frequency as 〈ϵ〉 increases. Furthermore, our model suggests that there is no simple power law relationship between the breakup frequency and 〈ϵ〉 because of the non-trivial dependence of p * (D e , D) on D and 〈ϵ〉.

Discussion
The classical KH framework hypothesized that bubbles tend to be broken by eddies of similar sizes in turbulence. So the only length scale used in the definition of the Weber number is the bubble diameter (D), and the contributions by the sub-bubble-scale eddies are ignored. To directly test this hypothesis, one would need to conceive an experiment that can cleanly disentangle eddies of different sizes. To that end, we developed an experiment by injecting a bubble into the flow driven by the head-on collision between the two vortex rings. As the rings break down into turbulence, the flow experiences two distinct stages with the early stage consisting of only the large vortical structures of the bubble size and the late fully-developed turbulence filled with a spectrum of the sub-bubble-scale eddies. In this flow, two distinct bubble breakup modes were observed. Against our intuition, the more violent and rapid secondary breakup occurred when the bubble was smaller, the flow kinetic energy was lower, and the resulting Weber number was much weaker. This observation clearly shows the inadequacies of employing only the Weber number and the bubble length scale, disapproving the key hypothesis in the KH framework. A model for bubble breakup in turbulence based on two new dimensionless numbers that account for two scales, the bubble and eddy sizes, is therefore proposed. Our framework supplement the classical framework that relies only on the stress criterion by adding a key criterion-the time scale. Our model is validated by the excellent agreement between the model prediction and experiments from two seemingly different flow configurations. Our work emphasizes the importance of the sub-bubble scale eddies for bubble fragmentation and adds a dimension to the existing framework on the bubble/drop dynamics in turbulence [19][20][21] .

Methods
Vortex ring. Vortex rings were generated using two identical piston-cylinder assemblies in a water tank with a size of 60 × 60 × 150 cm 3 . The cylinders and pistons were carefully adjusted to make sure they are aligned. Both cylinders have an inner diameter D 0 of 2.54 cm with the separation between their exits about 27.9 cm. Within each cylinder, a piston is driven by a stainless steel shaft connected to a pneumatic linear actuator. The piston can be moved forward and backward at different stroke times T s controlled by the air pressure. The stroke length L s was kept at 10.2 cm, and the stroke ratio was, therefore, SR = L s /D 0 = 4. Each assembly can generate vortex rings with different Reynolds numbers Re = V p D 0 /ν based on the mean velocity of the piston V p = L s /T s and the kinematic viscosity of the water ν. A hypodermic needle with a 1.7 mm inner diameter connected to a syringe filled with air was used to generate bubbles. By altering the injection rate through a computer-controlled syringe pump and needle size, bubbles with diameters ranging roughly from 2 to 5 mm can be generated. In this study, the initial bubble size is fixed at around 5 mm. Time delays and system control. One critical element of our experiments is the timing. Bubbles must be entrained into the vortices exactly at the moment when the two rings collide with each other. Therefore, three sub-systems, including two vortex-ring generators, one syringe pump for bubble injection, and four cameras were all controlled by a digital data acquisition (DAQ) system. Since t = 0 s is the reference time when two vortex rings began to touch each other, the syringe pump was activated at t 1 = −1.34 s for bubbles to naturally rise into the collision point; the vortex-ring generators were actuated at t 2 = −0.63 s for vortex rings to travel from the generators to the collision point; and the cameras were triggered at t 3 = −0.18 s to capture some frames before the collision. Even though the three time delays can be accurately controlled by our DAQ system, bubbles or vortex rings may not arrive at the exact same moment from one run to another. In order to quantify the uncertainty of the arrival time, we used high-speed cameras (7500 fps with 0.13 ms uncertainty) to capture the arrival time of rings and bubbles. The measured uncertainty of the arrival time difference between rings and bubbles is about ± 38 frames, corresponding to about ±5.1 ms. This number is much smaller than the time that it takes for the two vortex rings to break down into turbulence, i.e., about 100 ms, or the integral timescale of the turbulence generated, i.e., about 60 ms, which suggests that the synchronization is sufficiently accurate for reproducible experiments.
3D measurements. In each run, the shadows of bubbles and many surrounding tracer particles were projected onto all four cameras by the back LED panels. The images of individual phases were segmented based on the contrast and size differences. For the liquid phase, the 3D tracer trajectories were determined by performing the Lagrangian particle tracking with the in-house OpenLPT code 22 that implemented the parallelized Shake-The-Box algorithm 23 . These trajectories were then smoothed by convoluting them with the Gaussian kernels 24,25 , from which the particle velocity and acceleration could be obtained along their trajectories. For the gas phase, the bubble 3D geometry was reconstructed using the virtual-camera visual hull method 26 that calculates the intersection of the cone-like volume extruded from the silhouette on each camera. Once the 3D bubble geometry was obtained, the bubble center of mass was linked frame by frame to get the bubble trajectory.
Statistics. In total, we collected 83 datasets and observed about 150 breakup events for vortex ring Reynolds number ranging from 7.5 × 10 4 to 1.1 × 10 5 . The majority of the data were taken at the largest Reynolds number, in which almost each experiment yielded both primary and secondary breakups. For smaller Reynolds number, the turbulence generated is not sufficiently intense to fragment bubbles. So, these datasets were not included in the statistics. The breakups that were used in this manuscript are plotted as a point cloud in the Supplementary Fig. S7 to show the breakup locations relative to the flow.

Data availability
All the data supporting this work are available from the corresponding author upon reasonable request.