{'fraction_non_alphanumeric': 0.030363157368192884, 'fraction_numerical': 0.015123220062026878, 'mean_word_length': 4.402179787425689, 'pattern_counts': {'":': 0, '<': 4, '<?xml version=': 0, '>': 6, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 2, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Simultaneous measurements of deforming Hinze-scale bubbles with surrounding turbulence

We experimentally investigate the breakup mechanisms and probability of Hinze-scale bubbles in turbulence. The Hinze scale is defined as the critical bubble size based on the critical mean Weber number, across which the bubble breakup probability was believed to have an abrupt transition from being dominated by turbulence stresses to being suppressed completely by the surface tension. In this work, to quantify the breakup probability of bubbles with sizes close to the Hinze scale and to examine different breakup mechanisms, both bubbles and their surrounding tracer particles were simultaneously tracked. From the experimental results, two Weber numbers, one calculated from the slip velocity between the two phases and the other one acquired from local velocity gradients, are separated and fitted with models that can be linked back to turbulence characteristics. Moreover, we also provide an empirical model to link bubble deformation to the two Weber numbers by extending the relationship obtained from potential flow theory. The proposed relationship between bubble aspect ratio and the Weber numbers seems to work consistently well for a range of bubble sizes. Furthermore, the time traces of the bubble aspect ratio and the two Weber numbers are connected using the linear forced oscillator model. Finally, having access to the distributions of these two Weber numbers provides a unique way to extract the breakup probability of bubbles with sizes close to the Hinze scale.


Introduction
The process by which finite-sized gas bubbles and liquid droplets break in a turbulent environment constitutes one of the most fundamental and practically important phenomena in multiphase flows. Details of how this takes place have significant impact in various industrial and natural processes, such as chemical reactors (Jakobsen 2014), bioreactors (Kawase & Moo-Young 1990), air-sea gas transfer (Liss & Merlivat 1986), drag reduction (Lohse 2018;Verschoof et al. 2016), and two-phase heat transfer (Lu & Tryggvason 2008;Lu et al. 2005;Dabiri et al. 2013). Bubbles in strong turbulence can deform, break, and coalesce with each other. The presence of deformation adds to a problem that is already complicated even for the dispersed two-phase flows with rigid, non-deformable particles (Balachandar & Eaton 2010). Moreover, most works on bubble deformation have been limited to simulations (Elghobashi (2019) and the references within) with very few experimental works being able to resolve both phases in 3D simultaneously. It is thus the main objective of this paper to overcome this technical challenge and provide new experimental results to study bubble deformation and breakup in turbulence.
The earliest studies on bubble breakup in turbulence were conducted by Kolmogorov (1949) and Hinze (1955). In particular, Hinze unified the results of numerous preceding investigations. In his seminal work, he argued that only two dimensionless numbers are needed: one is the Weber number W e (also used by Kolmogorov (1949)) and the other one is the viscosity group, N = µ d / ρ d σD/2, in which ρ d and µ d are the density and the dynamic viscosity of the dispersed phase, respectively. σ is the surface tension, and D is the bubble diameter. This is the first time that the idea of the critical Weber number was introduced, and Hinze argued that the critical Weber number must depend only on N following W e crit = c(1 + f (N )), where f (N ) is a function of N . For bubbles with vanishing inner viscosity, the critical Weber number should just be a constant c. A critical Weber number of 0.59 was extrapolated from an earlier experiment conducted by Clay (1940). In this work, the Weber number is defined based on external stresses τ applied on the bubble surface W e = τ D/σ. τ is related to the energy dissipation rate (ǫ) in the form of τ = C 2 (ǫD) 2/3 based on the Kolmogorov theory, in which C 2 ≈ 2.13 is the Kolmogorov constant. This formulation should be, strictly speaking, only applied to homogeneous and isotropic turbulence, yet it has been used in many other flow configurations, including chemical reactors with impellers and jets based on the assumption of local isotropy.
Hinze's framework was constructed primarily for liquid droplets. But he noted that the critical Weber number should not be a universal constant; instead, it depends on the density difference between the two phases. Sevik & Park (1973) extended this framework to bubbles splitting in turbulence, in which a large density difference between the two phases was present. A slightly larger critical Weber number of 1.26 was observed. By assuming bubbles break once they start to resonate with surrounding turbulent eddies, the critical Weber number can be calculated analytically by equating the natural frequency of bubbles (Lamb 1932) with the reciprocal of the eddy turnover time. The predicted value seems to agree with their measured results. Although W e crit has been studied and reported in different types of flows, it should be noted that no consensus on W e crit has been reached so far. For air bubble breaking in different flow configurations, such as linear shear flow, turbulent jets, and homogeneous isotropic turbulence, a wide range of W e crit from 0.59 to 7.8 have been reported to date (Hinze 1955;Sevik & Park 1973;Deane & Stokes 2002;Martínez-Bazán et al. 1999a;Risso & Fabre 1998). Based on this observation, one can only conclude that W e crit is roughly of order unity.
Introducing W e crit also comes with a critical length scale. For a given mean turbulence energy dissipation rate ǫ ( . denotes ensemble average), the critical bubble size is often referred to as the Hinze scale D H , and it is related to W e crit following W e crit = ρC 2 ( ǫ D H ) 2/3 D H /σ. For a given W e crit , it is important to introduce the idea of energyabundant/super-Hinze (W e ≫ W e crit and D ≫ D H ) versus energy-limited/sub-Hinze (W e < W e crit and D < D H ) breakups. The former one has been studied much more extensively than the latter for a simple reason: super-Hinze breakup is much faster and more frequent so it is easier to observe in a finite volume and to collect enough statistics. Super-Hinze breakup is typically studied in several different flow configurations: pipe flow (Hesketh et al. 1991) and turbulent jets (Sevik & Park 1973;Martínez-Bazán et al. 1999b;Vejražka et al. 2018). In these cases, the energy contained in turbulent eddies is so abundant that each bubble is almost guaranteed to break-it is only a matter of time.
Breakup frequency is an important parameter in the population balance equation (Hulburt & Katz 1964;Ramkrishna 2000). However, this framework has one limitationit assumes that all bubbles above the Hinze scale will eventually break and no bubbles below the scale will ever break. This poses an important challenge to numerical simulations to account for sub-Hinze scale microbubbles, which are important to air-sea gas exchange (Deane & Stokes 2002), as well as underwater acoustics as these small bubbles tend to remain in the waterside for an extended period of time.
The breakup mechanisms that have been proposed and evaluated in the literature include (i) persistent stretching by straining flows (parallel flow, plane hyperbolic, axisymmetric hyperbolic, Couette flow, or rotating flow) (Hinze 1955); bubbles tend to exhibit regular affine deformation in these types of flows (lenticular or cigar-shaped). (ii) resonance mechanism that relies on bubble oscillation to siphon energy from the surrounding turbulence until breakup (Sevik & Park 1973;Hesketh et al. 1991;Risso & Fabre 1998). It typically assumes that the surrounding eddies retain a similar frequency with bubbles' natural frequency. (iii) Inertial mechanism relies on bubbles suddenly being exposed to strong flows, which leads to an almost-immediate irregular breakup. This mechanism has been studied in many contexts in addition to turbulence-induced breakup, e.g. raindrop fragmentation (Villermaux & Bossa 2009) and bag breakup in crossflows (Ng et al. 2008). In turbulence, three mechanisms may be all present, so applying only one mean Weber number to account for all breakup mechanisms is questionable.
In addition, as Risso & Fabre (1998) noted, the instantaneous and local Weber number, W e, could be much larger than the mean value W e . They proposed to use the time trace of W e along each bubble trajectory to evaluate its breakup frequency. However, in their experiments, such instantaneous Weber number was not directly accessible. As a result, flow velocity from single-phase turbulence was used as a surrogate. This is a common practice in the community as the simultaneous measurements of both phases, either in 2D or 3D, remain challenging.
To resolve deformation and breakup of the Hinze or sub-Hinze scale bubbles, in this paper, we will introduce an experiment that provides unique simultaneous measurements of both bubble deformation and surrounding flows thanks to the recent advancement of the 3D high-concentration particle shadow tracking (Tan et al. 2019) and 3D virtualcamera visual-hull shape reconstruction (Masuk et al. 2019a). In §2, the experimental setup, i.e. a vertical water tunnel system with a large section of homogeneous and isotropic turbulence, will be introduced. In the same section, the optical system designed to conduct simultaneous measurements of both the phases will also be discussed. In §4, based on the new datasets, we discuss how flow decomposition can be conducted to analyze the relative roles played by different mechanisms. In §4.5, we finally estimate the breakup probability of Hinze-scale bubbles in turbulence.

Experimental Setup
A facility was designed to accomplish two main goals: (i) maintain homogeneous and isotropic turbulence in a large volume to ensure that bubbles within this volume experience similar turbulence characteristics, and (ii) bubble deformation should be driven primarily by turbulence rather than by buoyancy, and bubble sizes remain close to the Hinze scale so that we can investigate the deformation and breakup of the Hinze-scale and sub-Hinze-scale bubbles. Satisfying both criteria is challenging. For example, many systems that feature a large region of homogeneous and isotropic turbulence tend to have Figure 1: Schematic of the V-ONSET vertical water tunnel; two insets show the 3D model of the jet array used to fire high-speed water jets into the test section and a bubble bank to inject bubbles, respectively. Additional details concerning this facility can be found in Masuk et al. (2019b) a low energy dissipation rate (Variano et al. 2004;Mercado et al. 2012) ( ǫ = O(10 −5 -10 −3 ) m 2 /s 3 ), whereas facilities that use water jets to break super-Hinze-scale bubbles can generate a large energy dissipation rate ǫ = O(0.1-10 3 ) m 2 /s 3 at the cost of having strong flow inhomogeneity and anisotropy (Martínez-Bazán et al. 1999b;Vejražka et al. 2018).
The experimental setup used for the current study is essentially a vertical water tunnel capable of generating turbulence with ǫ roughly at 0.16 -0.5 m 2 /s 3 . To extend the residence time of a Hinze-scale bubble in the interrogation volume, the mean flow in the tunnel was configured to move downward in a vertically-oriented test section. The flow speed was adjusted to balance the rise velocity of bubbles with diameters at around 3 mm to increase the residence time of these bubbles in the view area. Combined with ǫ in this region, W e was roughly at 1.19, indicating that most bubbles in the interrogation volume are close to the Hinze scale.
Turbulence in the test section was generated using 88 high-speed water jets (up to 12 m/s), each of which has a diameter d j of 5 mm, firing co-axially downward into the test section along with the mean flow. The firing pattern of these momentum jets was randomized in a way similar to the work by Variano et al. (2004) in order to ensure that no secondary flow structure would develop in the test section (De Silva & Fernando 1994;Srdic et al. 1996;Variano et al. 2004). On average, 12.5% of the jets were kept on at a time as this was found to maximize the turbulence intensity. The test section was set much farther downstream of the jets (about 80d j ) to ensure that the jets were well mixed and turbulence becomes homogeneous and isotropic with very little spatial variation. Additional details concerning this setup and its flow characteristics can be found in Masuk et al. (2019b).
Bubbles were generated at the bottom of the test section using two different sizes of hypodermic needles (Small needles: inner diameter (ID) of 160 µm and outer diameter (OD) of 300 µm; large needles: ID of 260 µm and OD of 500 µm). The range of bubble diameters in the experiment was 2-7 mm, which was set mostly by turbulence generated in our tunnel as large bubbles were broken before entering the interrogation window. Note that the bubble injection was far below the measurement volume to ensure that bubbles entering the measurement volume already lost any memory of the injection.
Both the bubble dynamics and turbulence statistics were collected by using six highspeed cameras each with a 1024×1024 pixel resolution and 4000 frame per second (fps) frame rate. The frame rate was selected to ensure that about 10 frames of images were taken within one Kolmogorov timescale τ η =2.5 ms. These cameras were spatially distributed to cover the entire perimeter of the octagonal test section. Six red LED panels with wavelength at roughly 630 nm were used to provide diffused backlighting to cast shadows of both particles and bubbles onto the imaging planes of all six cameras.

Flow characterization
Before discussing bubble deformation and breakup, single-phase turbulence statistics in the tunnel needs to be characterized to ensure that the flow is close to homogeneous and isotropic, and the bubble size is within the inertial range of turbulence. The details of these statistics can be found in Masuk et al. (2019b). Here, we only show the measured longitudinal second-order structure function i.e. D LL in figure 2. It can be seen that our experiments were able to resolve length scales as small as 2η, with η ≈50 µm being the Kolmogorov length scale η = (ν 3 / ǫ ) 1/4 . ν is the kinematic viscosity of water. We can resolve such a small scale thanks to our in-house high-concentration particle tracking system (Tan et al. 2019) that employs the Shake-The-Box method (Schanz et al. 2016).
The structure function should approach two limits: one in the dissipative range (r ≪ η) and the other in the inertial range (η ≪ r ≪ L). L is the integral scale, which is estimated based on L ≈ u ′3 / ǫ , where u ′ is the fluctuation velocity. The scale separation between η and L is determined by the Taylor-scale Reynolds number Re λ = 15u ′ L/ν, which is estimated to be around 435. In the dissipative range, the structure function follows the relationship of D LL = (ǫ/15ν)r 2 . In the inertial range, the 2/3-scaling law is based on the classical Kolmogorov theory. Although how long the inertial range is and if the Kolmogorov constant C 2 is affected by the finite-Reynolds number effect are subjected to further investigations (Ni et al. 2013), using a standard number C 2 = 2.13 can provide a reasonable estimation of ǫ. The solid line shown in the figure is based on the calculated ǫ = 0.16 m 2 /s 3 . However, if ǫ obtained from the inertial range is used to predict the dissipative range D LL (dashed line), it appears that the dashed line is systematically lower than the experimental results. In sum, the difference of ǫ estimated from either the dissipative or inertial range helps to quantify the experimental uncertainty of the mean energy dissipation rate: ǫ =0.22±0.07 m 2 /s 3 . Moreover, after bubbles getting injected into the system, bubbles can actively modulate turbulence and increase the local energy dissipation rate to around 0.52 m 2 /s 3 , which is calculated not from the structure functions but from the local velocity gradients that will be introduced in §4.2 and figure 4(b).
The shaded area in figure 2(b) marks the size range of bubbles with respect to the Kolmogorov scale η. As one can see, bubbles are well within the inertial range of turbulence, indicating that their deformation and breakup are indeed driven by the velocity fluctuations that can be estimated by the inertial range scaling.

Results and discussions
4.1. Simultaneous bubble and particle tracking As shown in figure 3(a), shadows of both bubbles and particles were projected onto the imaging planes of cameras. It is straightforward to separate their images based on the size difference. An example of segmented images of a bubble and surrounding tracer particles is shown in figure 3(b). The bubble silhouette was then input into a recently-developed virtual-camera visual hull method (Masuk et al. 2019a) for 3D shape reconstruction. Averaging surface points on the reconstructed geometry helps to determine the center of mass, which was then tracked in 3D to obtain a bubble trajectory. This procedure was repeated for all bubbles to acquire both the kinematics (from tracks) as well as geometrical information (from 3D shape reconstruction). On average, there were about 15 bubbles in the interrogation volume at each time instant, and each bubble trajectory roughly lasts about 0.09 seconds (360 frames) before exiting.
Separated images for tracer particles were input into our in-house OpenLPT (Tan et al. 2019) to perform the shake-the-box calculation (OpenLPT has already been open-sourced and is available for the entire community to use @JHU-Ni-Lab on Github). Compared with bubbles, significantly more particles could be found in the interrogation volume. At each time instant, there were about 6,000 tracer particles with the mean track length of about 200 frames. Figure 3(c) shows one example of about 40 tracer trajectories in the vicinity of a bubble with its 3D shape reconstructed from silhouettes segmented from figure 3(a). In this case, trajectories of tracer particles within 4D away from the bubble center were included. These tracks could be used to estimate the flow condition around the bubble. Since a high-concentration of tracer particles were available around almost every bubble, this experiment provided access to almost all key physical quantities-the Weber number,

Flow velocity and velocity gradient
For a bubble at location x 0 , its surrounding flow velocity u p (x 0 + x p ) can be measured at a number of discrete positions x 0 + x p where n tracer particles are located (p = 1, 2, ..., n). These tracer particles are sought within a radius of D s /2 from the bubble center with D s being the diameter of a spherical search volume. The flow field within this range can be decomposed into leading terms by applying the Taylor expansion: where u i = n p=1 u p i (x 0 + x p )/n represents the local mean flow estimated by averaging the velocity of n tracer particles. A ij and H jik indicate the velocity gradient tensor and the Hessian matrix, respectively, and the tilde denotes coarse-graining at the bubble size.
x p is the separation vector directed from the bubble center at x 0 to the p th tracer particle location. For small micro-bubbles with sizes in the dissipative range (D ≪ η), the flow is linear so the velocity Hessian is negligibly small. This higher-order term grows as a function of bubble size and eventually becomes important for bubbles with sizes in the inertial range (η ≪ D ≪ L). Although the velocity gradient around each bubble can be measured accurately (Ni et al. 2015), the velocity Hessian, on the other hand, requires measuring the gradient of the velocity gradient (three 3 × 3 matrices). Even though it is possible to calculate the velocity Hessian given sufficient number of tracer particles, the uncertainty becomes large due to the second-order spatial derivative. As a result, we limit only to the first two orders, i.e. the mean flow velocity u i and the velocity gradient A ij (x 0 ), to capture the key mechanisms of deformation.
The velocity gradient tensor A ij can be uniquely solved if we have four particles around a bubble. In practice, on average, 30-40 particles were used to perform least-squares fit by seeking the minimum value of the squared residuals Pumir et al. 2013;Ni et al. 2015). Although finite-sized bubbles typically come with a large search radius and abundant nearby tracer particles thanks to our tracking method (Tan et al. 2020), particles in the vicinity of a bubble are still randomly distributed in space. If nearby particles stay primarily within a quasi-2D plane, the estimation of the out-ofplane velocity gradient will have large uncertainty. Similar to previous studies (Xu et al. 2011;Ni et al. 2015), an inertia tensor I = p x p i x p j /tr( p x p i x p j ) was adopted to evaluate the shape factor of the particle cloud. If particles are uniformly distributed in 3D, three eigenvalues of this inertia tensor (γ i ) all equal to 1/3. For a quasi-2D distribution, the smallest eigenvalue (γ 3 ) will be very close to zero, and the gradient along that direction cannot be calculated. In practice, events with γ 3 /γ 1 smaller than 0.15 was therefore removed from the statistics.
Based on A ij , the coarse-grained rate-of-strain tensor, S ij , and rotation tensor, Ω ij can be directly obtained: Figure 4(a) shows the probability density function (PDF) of two eigenvalues of S ij (the largest λ 1 and the smallest λ 3 ) based on different D s . The PDFs of | λ 1 | and | λ 3 | overlap with each other for three D s considered, indicating that the magnitude of flow stretching and compression near a bubble on average is similar. The PDF progressively shifts leftward as D s becomes larger because coarse-graining velocity gradients at a larger D s works effectively as enlarging a low-pass filter, which will continue to reduce the gradient as D s increases. As a result, the calculated velocity gradient using particles within a search diameter of D s should always underestimate A ij at the bubble scale D because D s > D. Fortunately, both D s and D are in the inertial range, and the eigenvalues of A ij can be related to the local energy dissipation rate in the form of C 2 ( ǫd) 2/3 = ( λ 3 d) 2 , where C 2 = 2.13 is the Kolmogorov constant (Batchelor 1953;Sreenivasan 1995;Ni et al. 2013) and ǫ is the coarse-grained energy dissipation rate for a range of length scales d considered.
To check if the distribution of ǫ is indeed the same for d varying between D to D s , in figure 4(b), the PDFs of the estimated local ǫ using three different D s are shown. Although the distribution of the calculated λ 3 are sensitive to D s , once converted to ǫ, three curves from all three D s fall right on top of each other, indicating that the local ǫ is roughly the same for the range of D s considered. Therefore, although D s > D is needed to include enough tracer particles for calculating velocity gradients, the statistics reported are insensitive to D s thanks to the universal inertial range scaling in homogeneous and isotropic turbulence.
The coarse-grained energy dissipation rate can be described by the log-normal distribution based on the Kolmogorov refined theory in 1962 (Kolmogorov 1962) and multi-fractal spectrum (Meneveau & Sreenivasan 1991).
where ǫ r is the energy dissipation rate coarse-grained at a scale r. µ ≈ 0.25 is the intermittency exponent. L = 3.2-6 cm is the integral length scale. A is a parameter that needs to be fitted to the experimental data to determine the variance of ǫ when r = L, which was found to be around one. Based on the definition, ǫ measured from our experiments is equivalent to ǫ r | r=D , which is shown as the black solid line in figure 4. The nice agreement between the experimental data and the log-normal distribution (equation 4.2) shows that the measured coarse-grained energy dissipation rate is consistent with the classical Kolmogorov theory.

Bubble deformation by the velocity gradient W e vg
In turbulence, the difference of dynamic pressure across a bubble acts to push the bubble interface inward to drive bubble deformation. Based on this argument, λ 3 , which is associated with the direction that compresses the most, should be the more relevant eigenvalue of A ij . Following the argument, the Weber number can be defined as This Weber number definition is based on the local coarse-grained A ij and ǫ, which is different from the mean Weber number defined by Kolmogorov (1949) and Hinze (1955). log-normal distribution of the local ǫ explains the observed shape of the PDF of W e vg , from which bubble breakup probability can be determined.

Slip-velocity induced deformation W e slip
As Hinze stated in his original seminal work (Hinze 1955), employing the velocity gradient to evaluate the deformation and breakup of droplets should only be applied if there is no large density difference between the dispersed phase and the carrier phase. For bubbles in water, such a large density difference does exist, and it is not surprising that W e vg may not capture the total stress acted on bubbles by turbulence. For example, the instantaneous velocity mismatch between the two phases could also lead to significant dynamic pressure that needs to be evaluated. This effect can be captured by the so-called slip velocity, u slip = u b − u f . As its name suggests, u slip quantifies the drift of a bubble of velocity u b away from the instantaneous local flow velocity u f . u f represents the continuous-phase fluid velocity at the center of a bubble if the bubble was not there. In practice, u f needs to be estimated from the continuous-phase velocities measured in the vicinity of the bubble. Therefore, we assume u f to be the same as u i in equation 4.1, which can be estimated by averaging the tracer velocities around the bubble. Figure 6(a) shows the PDF of only one horizontal component of u f normalized by its own standard deviation. u f can be calculated around bubbles of different sizes, which are shown by different colored symbols. The solid line indicates the standard normal distribution, which seems to agree well with the horizontal velocity distribution of bubbles of all sizes, at least for the range of bubble sizes considered. Since the PDFs of bubbles of all sizes are nearly the same, they can be combined and the results are shown in figure 6(b). Furthermore, to rule out the possible D s effect, the same procedure was repeated for three different D s =2-4D to 6-8D. As shown in figure 6(b), no discernible difference is observed for the flow velocity PDF at three D s , which suggests that u f is not sensitive to D s either. u b denotes the bubble velocity with one of its horizontal components along the x-axis being u b,x . Figure 7(a) shows the distribution of u b,x , normalized by its own standard deviation, for a wide range of bubble sizes, and the distribution for all bubble sizes seem to agree with a Gaussian distribution (solid line) very well. The standard deviation of u b for both horizontal directions are shown as dashed lines in figure 7(b), and they exhibit a weak, if at all, dependence on D. Note that, in the other limit for bubbles rising in a quiescent medium with no turbulence, since the horizontal velocity is coupled with the size-dependent rise velocity ( will eventually cross at around 200 mm/s for bubble size close to zero, which gives the right limit as extremely-small bubbles should behave similarly to tracers u 2 b,x Although both u f and u b along the horizontal directions appear to follow the Gaussian distribution, the slip velocity u slip = u f − u b does not. As shown in figure 8(a), for bubbles of all sizes, the tails of the slip-velocity PDF (u slip,x ) are systematically higher than that of the Gaussian function (black solid line), indicating that the slip velocity is more intermittent than the velocity of either phase alone. For the distribution of the normalized slip velocity, similar to the PDFs of u f,x and u b,x , no obvious bubble-size dependence is observed. Note that the PDF of u slip resembles that of the velocity increment between two points in single-phase turbulence (Kailasnath et al. 1992;Sreenivasan 1999;Li & Meneveau 2005). The PDF of the velocity increment has been fitted with a stretched exponential function (Kailasnath et al. 1992), which is adopted here to describe the observed PDF of the slip velocity.
where C is the normalization factor, and Q and m are fitting parameters in the stretched exponential function. For single-phase turbulence, the degree to which the tail of the PDF is stretched depends on the scale separation. If the velocity separation is close to the integral length scale, the PDF recovers the Gaussian distribution (m=2). As the separation becomes smaller and smaller, the PDF becomes more and more intermittent; at m = 1, the PDF follows an exponential function. If we take the bubble size 0.03L to 0.12L as the scale separations to calculate the velocity increment in single-phase turbulence, the scaling exponent m should vary between 0.8 to 1.05 based on the work by Kailasnath et al. (1992). In our case, although the slip velocity distribution also follows the stretched exponential, the PDF preserves its shape for all bubble sizes considered in this work with no obvious scale dependence, and all symbols in figure 8(a) collapse with one another. Therefore, the distributions of the normalized slip velocity for different sizes of bubbles were fitted together with one stretched exponential function and one set of constants, i.e. Q, and m. In particular, m is found to be a constant close to 6/5, which is slightly larger than the range of m from 0.8 to 1.05 in single-phase turbulence. This observation suggests that the slip velocity between the two phases is less intermittent compared with the velocity increment between two points in single-phase turbulence under the same scale separation, which is not surprising since bubbles are capable of filtering out intermittent small-scale fluctuations.
Furthermore, the fluctuation slip velocity ( u 2 slip 1/2 ) increases as a function of bubble size D, suggesting that larger bubbles with a larger inertia tend to deviate further away from the surrounding fluid velocity. At the same time, the typical velocity scale of an eddy of the bubble size D also increases with D, following ( ǫ D) 1/3 . After assuming that these two velocity scales are related, the measured u 2 slip 1/2 along all three directions are fitted with γ( ǫ D) 1/3 by performing the least-square regression to obtain the fitting coefficient γ, which turns out to be 0.62. The fitted result is shown in figure 8(b) as black dash-dotted line. It is clear that the fit reproduces the growth of the measured standard deviation of the slip velocity as a function of D, but the agreement between the fitted and the measured results is not perfect. Nevertheless, for simplicity and without any other alternative velocity scales, this fit using the eddy velocity is used to estimate u 2 slip 1/2 for bubbles with size in the inertial range. With this relationship and two coefficients, i.e. Q = 3/4, m = 6/5, the distribution of u slip can be estimated from equation 4.4.
Finally, the distribution of W e slip , calculated based on equation 4.4 and 4.5, is shown in figure 5. The blue solid line indicates the predicted W e slip based on the stretched exponential fit to the horizontal slip velocity u slip,x (equation 4.4). The distribution also peaks at around W e ≈ 1, which is slightly smaller than the most probable value of W e vg . The right tails of both PDFs (W e vg and W e slip,x ), corresponding to the range of W e that is important for deformation and breakup, are very close to each other. This may suggest that, for bubble deformation, slip velocity and velocity gradient may be equally important; completely relying on the velocity gradient may not account for all stresses that bubbles experience in turbulence.

Buoyancy-induced deformation
Although the turbulence energy dissipation rate has been set as high as possible in our facility, the buoyancy effect is not negligible. In figure 5, the PDF of W e slip in the vertical direction based on the z-axis slip velocity, i.e. W e slip,z is also shown. This PDF has a bump near W e slip,z ≈3-4 because of the buoyancy effect, but both the left and right tails seem to agree with those of W e slip,x . This suggests that both the turbulence effect and the buoyancy effect are present in W e slip,z , but the effect of buoyancy is rather limited to a comparatively narrower region near the peak of the PDF. Nevertheless, the exact functional form of the PDF close to the peak is unknown and requires further investigations to understand the coupling between the local bubble rise velocity and the surrounding turbulence.
Note that W e slip,z is similar to the Eötvös number: Eo = ρgD 2 /σ, as the terminal vertical slip velocity u slip,z driven primarily by buoyancy should be proportional to √ gD. Note that this relationship is approximate, as the buoyancy-driven terminal rise velocity is also sensitive to the bubble geometry, orientation, and the drag coefficient. In intense turbulence, these parameters could also be functions of ǫ. In a recent paper (Salibindla et al. 2020), the drag coefficient of bubble with different sizes in intense turbulence was reported, and it follows C D = max(24/Re b (1 + 0.15Re 0.687 b ), min(f (Eo), f (Eo)/W e 1/3 )) where f (Eo) = 8Eo/3(Eo+ 4). Based on this equation, the most probable slip velocity in the vertical direction can be calculated following: u 2 slip,z = 2V b (ρ − ρ b )g/ρAC D where V b is the volume of a bubble, A is the projected area of volume-equivalent spherical bubble, and ρ b is the density of bubble. For the bubble size range considered, W e slip,z calculated based on C D is about 4, which is consistent with the bump of W e slip,z observed in the PDF. This agreement confirms that the observed bump in the distribution of W e slip,z is indeed driven by buoyancy.
All together, it seems that the bump in the distribution of W e slip,z is limited to a narrow range, and the right tail of W e slip,z seems to be close to that of W e slip,x and W e vg . This suggests that, at least for our parameters when ǫ ≈ 0.2-0.5 m 2 /s 3 , the buoyancy-induced deformation is limited. If we keep increasing ǫ , the buoyancy effect will become even weaker.

Bubble aspect ratio vs. Weber numbers
So far, we have focused primarily on discussing the distribution of different definitions of Weber numbers and understanding the connection between these Weber numbers and their associated turbulence characteristics. In this section, the instantaneous Weber numbers along bubble trajectories will be used to study the mechanisms of bubble deformation and breakup in turbulence. Figure 9 shows two examples of the simultaneous measurements of bubble aspect ratio (α) and the Weber numbers (W e vg and W e slip ). α was calculated as α = r 1 /r 3 where semi-major axis r 1 and semi-minor axis r 3 are the instantaneous maximum and minimum radius of a bubble, respectively. These two axes were determined by finding the maximum and minimum vertex-center distances from the 3D-reconstructed bubble geometry, respectively. As shown in figure 9(a), the reconstruction seems to successfully capture the oscillation of a bubble that undergoes small-amplitude deformation. For this particular case, both of the Weber numbers for almost the entire duration are smaller than 5. The time trace of α is not similar to that of either W e at first glance. The only evident correlation is probably at t = 10-20 ms, when a small peak observed in the time trace of α seems to be driven by a large W e slip maybe 5 ms earlier. For t = 50-70 ms, despite W e slip drops close to zero, α continues to rise thanks to a relatively large value of W e vg during this period. This indicates that bubbles probably respond to both Weber numbers, likely to be the maximum instantaneous Weber number, i.e. max(W e vg ,W e slip ).

Simultaneous measurements of bubble geometry and W e
The simultaneous measurements also provide a framework to test models for bubble deformation and breakup. One such a model has been proposed before by Risso & Fabre (1998) and Lalanne et al. (2019), which is essentially a forced oscillator model that connects bubble deformation directly to local W e through a linear differential equation. This model is designed to follow the interaction between a bubble with surrounding turbulent eddies along its Lagrangian trajectory, exactly how our experiments were performed. The dimensionless form of the equation is: where ξ = 1/2πτ d f 2 is the damping coefficient, τ d = D 2 /80ν is the damping time scale (Risso & Fabre 1998), and f 2 = 96σ/ρD 3 is the mode 2 natural frequency of bubble oscillation (Lamb 1932).
We recognize that this model is linear, and the deformation and breakup process of bubbles in turbulence is nonlinear, especially when bubbles exhibit non-affine deformation and subsequently break. Linearizing this problem relies on two assumptions: (i) The Weber numbers are not very large ( O(1)); and (ii) Bubble deformation is driven primarily by eddies of the bubble size. For (i), similar to the roles played by the Reynolds number in laminar-turbulence transition in single-phase pipe flows, the Weber number here should determine when the process becomes nonlinear. Since there is no consensus on the range of Weber number where the nonlinearity becomes important, we can only argue that the process is linear or nonlinear if the Weber number is much smaller or much larger than O(1). Note that if the flow is turbulent or if the flow Reynolds number is large is irrelevant here. For example, for a small bubble with size smaller than η in turbulence, the bubble only senses linear flows around itself even though the flow Reynolds number is large, so the deformation process of this bubble is always linear as long as W e 1. In this work, the Weber number has a wide distribution. For most Hinze-scale bubbles in our experiments with the Weber numbers close to one, the linear model should capture some of the key dynamics. But for highly-deformed bubbles with Weber numbers in the order of O(10) to O(10 2 ), the linear model is not expected to work, but the discrepancy between the model prediction and measured results may still shed new light on the problem of bubble deformation and breakup in turbulence. For (ii), the Weber numbers defined based on flows of the bubble size essentially ignore the contributions from subbubble-scale eddies. Statistically, this assumption probably holds because smaller eddies tend to be weaker, even though they could occasionally become exceedingly strong due to turbulent intermittency. But it would require further investigations to understand their contributions.
In equation 4.6, the amplitude of the instantaneous W e(t) is controlled by the prefactor K ′ . W e(t) was not available before in other experiments, and it had to be estimated based on two-point velocity measurements from single-phase turbulence (Risso & Fabre 1998). In this work, in addition to measuring W e(t) directly, the method also allows us to distinguish between W e vg and W e slip . But since the model did not explicitly account for individual W e, here we apply three different inputs: W e(t) = W e vg , W e(t) = W e slip , and W e(t) = max(W e vg , W e slip ) to equation 4.6 to obtain model predictions, which are shown in figure 9 and 10 as blue dashed line, red dashed line, and black solid line, respectively.
Note thatâ in equation 4.6 is defined as the ratio of the deformed radius (a = R − D/2) to D, where R is the major axis of the deformed bubble and D is the diameter of an volume-equivalent sphere. It has to be converted to α * = 2(âD + D/2)/D for comparisons with the measured α. However, despite our best efforts, α = α * because α * does not contain information about the minor axis, which has to be replaced with D/2. Nevertheless, α and α * should share the same trend as bubbles deform.
Indeed, similarities can be observed between α and α * in figure 9(a). The magnitude of α * is affected by the prefactor K ′ in equation 4.6, which is fixed at 0.1 in this work. In figure 9(a), the model successfully captures roughly four oscillation periods, which can also be observed in the experimental results. This agreement suggests that the model can explain the deformation and shape oscillations for some bubbles undergoing small-amplitude deformation with small Weber numbers. To be more precise, the Weber numbers, including both W e vg and W e slip , are around 1-3, and the resulting bubble aspect ratio is about 2, which is considered as small-amplitude linear deformation. Note that this definition of small-amplitude deformation and small Weber number are completely based on our observation, and the transition between linear and nonlinear deformation is likely to be smooth over a large range of Weber numbers without having a clean demarcation.
For figure 9(a), there is a small gap in the time trace with no data at around 100 ms because, during this time period, the velocity gradient calculation does not meet the requirement mentioned in §4.2. In addition, despite the overall agreement, the phase lag of each period of α keeps changing in the experimental data, e.g. the second peak is much closer to the first one and further away from the third one. This varying phase lag is a feature that cannot be reproduced from the model. Figure 9(b) shows another example to compare α with α * . A large α observed at 32 ms seems to correlate with an event of large W e vg occurred at 20 ms, whereas a small bump of α at 60 ms seems to correspond to a sudden increase of W e slip at 50 ms. This observation is still consistent with the argument that bubble deformation tends to be driven by both Weber numbers. In this case, three model-predicted time traces of α * differ from one another; W e vg is systematically larger than W e slip for the entire duration. Nevertheless, for this case, although the model-predicted time trace still embraces some oscillatory features, the measured results do not. Over a similar period of time compared with figure 9(a), only one distinct peak is observed in figure 9(b). This suggests that the linearoscillator model provided by Risso & Fabre (1998) may capture the dynamics of bubbles undergoing small-amplitude deformation (α ≈1-3) with small W e (W e 3), like in figure 9(a), but not for all bubbles, particularly not for bubbles with large W e undergoing large aspect ratio changes. In addition, sometimes the surrounding flow maintains its strength, and the bubble is not allowed to oscillate freely. For these cases, the oscillation amplitude becomes smaller, and the model tends to overpredict the bubble aspect ratio.
Since the model provided by Risso & Fabre (1998) is primarily designed to characterize the breakup process, figure 10 shows two examples of bubbles that eventually break. W e for breaking bubbles is clearly much larger: one reaches close to 40 and the other one climbs up to almost 20, nearly a factor of 4-8 larger than the cases for smallamplitude deformation in figure 9(a). In figure 10(a), W e slip dominates, but a local event at W e slip ≈ 40 did not break the bubble, even though it did successfully deform the bubble to a large α at around 8. Following a peak of W e slip at 85 ms, this bubble eventually broke at t = 90 ms, with the instantaneous W e about 10 and local aspect ratio close to 2.
This example represents many bubbles we observed that do not break at the moment when α reaches its peak; instead, they split at a later time during the process of retraction. As bubbles retract towards a spherical shape, the excess surface energy stored on the bubble interface is transferred back to the surrounding flows in the form of turbulent kinetic energy (Dodd & Ferrante 2016). However, this process is unstable because of the large density difference between the two phases, and it eventually leads to breakup before bubbles return back to a spherical shape.
For the first example ( figure 10(a)), W e slip is intermittent with a large variation of magnitude in a short period of time, which seems to be consistent with the notion of eddy-bubble collision (Prince & Blanch 1990;Risso & Fabre 1998) that this bubble keeps encountering different eddies with varying intensity. For the second example shown in figure 10(b), both Weber numbers slowly increase with time until the bubble breaks. The aspect ratio does not vary much throughout the entire time trace. For the last 20 ms, the aspect ratio of this bubble is close to a constant. Rather than eddy-bubble collision, the results seem to suggest an alternative mechanism: bubbles entrained in an eddy slowly get pulled apart by this eddy as it grows in strength over time.
The model predictions are also shown alongside with these two examples of breakup. In both cases, oscillations clearly observed in α * from the model calculation are not so obvious in experimental results. As explained before, it is not surprising as large deformation is expected to be highly nonlinear and should deviate from the linear equation 4.6. In addition, this disagreement also implies that the bubble oscillation may not be the right mechanism for bubble breakup at large Weber numbers. Bubbles might just be stretched and deformed by the local strains and slip velocity until the surface tension could no longer hold it together. Such a mechanism that needs to account for 3D couplings between bubbles and surrounding flows is clearly missing in the current model framework. This calls for future investigations into improving the model for large Weber numbers and other breakup mechanisms.

Distribution of bubble aspect ratio
In addition to the response of individual bubbles to different Weber numbers, the distribution of α could also be connected to that of W e to examine if the bubble aspect ratio can be solely determined by W e in a statistical sense. This relationship between α and W e was first derived by Moore (1965) for bubbles rising in a quiescent medium: where W e was defined to account for the dynamic pressure driven by the rising motion of bubbles, not by turbulence. In addition, the key assumption in this model is that W e ≪ 1 so that the departure from a spherical shape is so small that any high-order terms associated with W e 2 can be ignored.
For W e ≈ 1 or above, the potential flow theory applied to oblate ellipsoids with fore-aft symmetry yields: This equation has a maximum aspect ratio of 6 when the Weber number is close to 3.745, above which the symmetric shape is impossible to attain for a bubble. Although this formulation provides a better framework for our studies of bubbles with W e ≈ 1, it cannot predict the relationship between α and instantaneous W e for W e > 3.745, which is about 21.5% of the total events in our experiments. Figure 11(a) shows the PDF of α for bubble size D = 4.5 mm. The PDF peaks at α = 1.7 and has a long tail that skews towards larger values of α. On top of the experimental results, the PDF of α calculated from equation 4.7 using W e = W e vg as the input is also plotted as the red dashed line. Although the peak location is slightly different from the experimental results, the overall trend is close. Since the long tail in the PDF of W e vg comes from the log-normal distribution of ǫ, a large probability of strong deformation α > 2 is likely to be contributed by intermittent events with a large ǫ. If α versus W e vg follows a linear relationship, the two PDFs should overlap with each other. To make the red dashed line closer to the experimental results, we tried to add a second-order correction to equation 4.7, which did not provide satisfactory results (not shown here). Finally, after adjusting the parameters in equation 4.7, we settle down to a simple new equation equation 4.9. Similar trend of α can be seen even with this Weber number. But compared with W e vg , the results based on W e slip tend to underpredict α, which is consistent with the observation in figure 5 that the peak of W e slip PDF is on the left side of W e vg PDF. Nevertheless, the right tail of α can be reproduced by the calculations using both W e slip and W e vg . This may imply that the dynamic stresses contributed by both velocity gradients and the slip velocity are equally important, but turbulent velocity gradients seem to work better and thus more important for mild deformation. Figure 11(b) compiles the PDFs of α for five different bubble sizes, from 2.5 mm to 6.5 mm with an interval of 1 mm. The peaks of these PDFs progressively shift rightward towards a larger α as D grows, which is consistent with our intuition that large bubbles are more deformable and thus have a larger α on average. Moreover, the PDF becomes wider (the right tail of the PDF rises) as D increases, which implies that the probability of bubbles with α much larger than the mean also increases. To explain this, equation 4.9 is applied to all these cases with different D, and results are shown as solid lines with corresponding colors to compare with the PDF of measured α. The modeled PDF of α agrees with the measured results really well for most bubble sizes except for the largest bubbles where the buoyancy effect may deform bubbles even further. This agreement suggests that the observed change of the PDFs of α as a function of D is driven mostly by the change of W e, but the relationship between α and W e may not be linear for bubbles deformed by turbulence. An alternative method to predict the relationship between α and W e is to use the model provided by equation 4.6. As discussed before, although the model-predicted time trace of α does not match with the measured one exactly, for small Weber numbers, the model is still able to capture some bubble oscillatory deformation, as shown in figure  9(a). Here, we want to extend the test beyond single time traces. The comparison is shown in figure 12 for only two sizes of 2.5 mm and 5.5 mm for simplicity. The model seems to reproduce the overall trend of the PDF. But the PDF of α * calculated from the model is more flat than that of the measured α, indicating a much higher probability of strongly-deformed bubbles compared with the measured results. This observation is consistent with the time traces shown in figure 10 that the linear-oscillator model seems to overpredict the number of large-deformation events.
This observed difference can also be attributed to other possible reasons. For example, equation 4.6 is a linear one-dimensional model with both W e and α * being scalars. In our experiments, W e vg has an implicit direction that follows the largest compression direction of the rate-of-strain tensor, and W e slip should be aligned with the slip velocity direction. Their combined effect to bubble deformation may not follow a simple relationship of max(W e vg , W e slip ). In certain circumstances, they could potentially work against each other, which is not accounted for in the linear-oscillator model.

Breakup probability
One condition that is implicitly assumed in many breakup models (Martínez-Bazán et al. 1999b) is that all bubbles will eventually break, just a matter of time, as long as W e > W e crit . This implies a breakup probability (p b ) close to 100%, which should be valid for large W e ≫ 1 and D ≫ D H . However, for Hinze-scale or sub-Hinze-scale bubbles with W e 1 and D D H , p b could be anywhere from 0 to 100%. This number has not been reported before from previous experiments as it is challenging to estimate p b given the limited residence time of bubbles within the view area. The data that has been generated in this work provides a unique way to evaluate p b indirectly based on two assumptions: (i) turbulence remains close to homogeneous and isotropic so that the statistics collected in the entire view area can be compiled together to predict the breakup probability; and (ii) the local Weber number is the sole parameter that determines the status of bubble deformation and breakup. Introducing and measuring local W e is one step further from the Hinze's seminal work, in which the ensembleaveraged W e was adopted to quantify the breakup probability. The limitation of using W e is that, based on W e being larger or smaller than the critical W e crit , p b is close to a step function (p b = 1 if W e > W e crit ; p b = 0 if W e < W e crit ). However, in turbulence, local flows could be orders of magnitude stronger than the mean; bubbles could break in response to the local W e instead of the mean Weber number. To transfer this intuition to quantitative results, the main objective of this section is to determine p b by linking local W e distribution to bubble breakup probability.
Before estimating p b , we would like to extend the PDF of W e vg and W e slip beyond our experiments to other turbulent flows with different ǫ and bubble size D. Based on equations 4.2 and 4.4, the distribution of local W e for three ǫ from 0.1 m 2 /s 3 to 10 m 2 /s 3 are shown in figure 13(a). Both W e slip and W e vg shift rightward systematically with a seemingly-identical shape. It is important to note that the shape remains the same on the logarithmic scale, indicating that the distribution actually widens on the linear scale and the standard deviation of local W e increases as ǫ grows.
From the distribution, we can estimate p b based on: Here p b is written as a function of the mean Weber number W e to help with other experiments or simulations that have access only to W e from the mean energy dissipation rate. It can be seen that, as W e grows, either due to a larger D or larger ǫ, the bubble breakup probability p b will increase. Figure 13(b) shows p b as a function of W e . In most previous works assuming that the breakup probability has a sharp transition at W e crit , p b would behave like a step function: p b = 1 for W e W e crit and p b = 0 for W e < W e crit . Although these two limits still apply in figure 13(b), the transition is much smoother, spanning over a few orders of magnitude of W e . Note that the choice of W e crit does not affect the shape of the curve. As shown in figure 13(b), when we change W e crit from 1 to 4, it just shifts the transitional W e towards the new W e crit without affecting the overall trend.
This framework applies to W e vg and W e slip , both of which contribute to bubble breakup. Since the right tail of their respective distribution is very close to each other, an equal contribution from the two breakup mechanisms was assumed. As the result, the curves of p b for either W e alone approach 0.5 for W e much larger than W e crit to ensure that the sum of the two p b equals to one. Furthermore, even for the same W e crit , p b of W e slip (blue lines) stays mostly above that of W e vg (red lines) until they cross at a location very close to the plateau near p b =0.5. This difference can be ascribed to the difference of the PDFs: compared with W e vg , the PDF of W e slip has a larger probability for small W e. This also suggests that bubbles close to the Hinze scale may be deformed more often by the slip velocity.
Finally, the total breakup probability p b by summing the contribution from W e vg and W e slip and using W e crit = 1 to 4 are shown as three cyan lines in figure 13. These three lines are fitted with the same switch function: p b = [1 + (2.8 W e /W e crit ) −1.7 ] −1 (4.11) which includes W e crit as the input. The fitted results are shown in figure 13(b) as three black solid lines. One may not see the cyan lines at all because the fit overlap perfectly with the calculated results over the entire W e range for three W e crit considered. Equation 4.11 provides a method to estimate bubble breakup probability in turbulence, particularly for bubbles close to the Hinze scale and W e ≈ 1.

Conclusion
Bubble deformation and breakup in intense turbulence is ubiquitous in many applications, but details of how this takes place for a bubble close to the Hinze scale remain elusive because of the lack of data to probe the interaction between finitesized bubbles and surrounding turbulence. In this study, both 3D bubble geometry and nearby 3D particle tracks were acquired simultaneously using our in-house virtual camera reconstruction and particle tracking algorithm. The experiments were performed in a system that can reach a high turbulent energy dissipation rate that can significantly deform and even break bubbles, while maintaining homogeneous and isotropic turbulence throughout the entire measurement volume. Since the 3D information of both phases is available, this unique dataset allows us to interrogate the couplings between the two phases, in particular the key mechanisms that drive bubble deformation and breakup. The flow velocity was decomposed into two components, the local flow velocity and velocity gradient, both coarse-grained at the bubble scale. Each component can be used to define its own Weber number as a way to quantify their relative contributions to bubble deformation. In this study, in addition to directly measuring the Weber numbers, bubble deformation is also connected to the log-normal distribution of the local coarse-grained energy dissipation rate ǫ. The modeled distribution of both ǫ and W e vg based on turbulence characteristics agree well with the measured results. Moreover, because of the density mismatch between the two phases and the finite bubble size effect, the slip velocity also plays an important role. Based on this observation, a different Weber number is defined to measure deformation driven by the slip velocity, whose distribution can be fitted with a stretched exponential function inspired by the distribution of two-point velocity increments in single-phase turbulence. Based on this function, the distribution of the slip-velocity-based Weber number can be connected to ǫ and bubble size.
The distribution of the Weber number was also connected to the reconstructed bubble geometry. It has been shown that the relationship that was developed for describing bubbles rising in a quiescent medium does not work well for the turbulent case. A new non-linear model was proposed to improve the fit, and it seems to work well for a range of bubble sizes considered. In addition, the results were tested against a linear forced oscillator model that was proposed before. Although the model does seem to reproduce some key features of a few example time traces qualitatively, the distribution of the predicted aspect ratio does not match with the directly-measured results quantitatively.
Finally, the Weber number distribution is generalized for different bubble sizes and energy dissipation rates in order to evaluate breakup probability, which was estimated based on the mean energy dissipation rate in many other works. In contrast to what has been proposed before that the bubble breakup probability experiences a precipitous drop as bubble size decreases below the Hinze scale, accounting for the distribution of local Weber number helps to smooth the curve near the Hinze scale. The final calculated relationship between breakup probability and the mean Weber numbers was fitted with a simple function that can help future works to estimate bubble breakup probability based on the mean Weber number.
We acknowledge the financial support from the National Science Foundation award numbers: 1854475 and CAREER-1905103. We would also like to acknowledge Charles Meneveau for suggestions. Declaration of Interests: The authors report no conflict of interest.