1. 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 Reference Jakobsen2014), bioreactors (Kawase & Moo-Young Reference Kawase and Moo-Young1990), air–sea gas transfer (Liss & Merlivat Reference Liss and Merlivat1986), drag reduction (Verschoof et al. Reference Verschoof, Van Der Veen, Sun and Lohse2016; Lohse Reference Lohse2018) and two-phase heat transfer (Lu, Fernandez & Tryggvason Reference Lu, Fernandez and Tryggvason2005; Lu & Tryggvason Reference Lu and Tryggvason2008; Dabiri, Lu & Tryggvason Reference Dabiri, Lu and Tryggvason2013). 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 Reference Balachandar and Eaton2010). Moreover, most works on bubble deformation have been limited to simulations (see Elghobashi Reference Elghobashi2019, and references therein) with very few experimental works being able to resolve both phases in three dimensions 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 (Reference Kolmogorov1949) and Hinze (Reference Hinze1955). 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 
$We$ (also used by Kolmogorov Reference Kolmogorov1949) and the other is the viscosity group, 
$N=\mu _d/\sqrt {\rho _d\sigma D/2}$, in which 
$\rho _d$ and 
$\mu _d$ are the density and the dynamic viscosity of the dispersed phase, respectively. 
$\sigma$ 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 
$We_{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 (Reference Clay1940). In this work, the Weber number is defined based on external stresses 
$\tau$ applied on the bubble surface 
$We=\tau D/\sigma$. Here 
$\tau$ is related to the energy dissipation rate (
$\epsilon$) in the form of 
$\tau =C_2(\epsilon D)^{2/3}$ based on the Kolmogorov theory, in which 
$C_2\approx 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. 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 (Reference Sevik and Park1973) 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 Reference Lamb1932) with the reciprocal of the eddy turnover time. The predicted value seems to agree with their measured results. Although 
$We_{crit}$ has been studied and reported in different types of flows, it should be noted that no consensus on 
$We_{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 
$We_{crit}$ from 0.59 to 7.8 have been reported to date (Hinze Reference Hinze1955; Sevik & Park Reference Sevik and Park1973; Risso & Fabre Reference Risso and Fabre1998; Martínez-Bazán, Montanes & Lasheras Reference Martínez-Bazán, Montanes and Lasheras1999a; Deane & Stokes Reference Deane and Stokes2002). Based on this observation, one can only conclude that 
$We_{crit}$ is roughly of order unity.
Introducing 
$We_{crit}$ also comes with a critical length scale. For a given mean turbulence energy dissipation rate 
$\langle \epsilon \rangle$ (
$\langle \cdot \rangle$ denotes ensemble average), the critical bubble size is often referred to as the Hinze scale 
$D_H$, and it is related to 
$We_{crit}$ following 
$We_{crit}=\rho C_2(\langle \epsilon \rangle D_H)^{2/3}D_H/\sigma$. For a given 
$We_{crit}$, it is important to introduce the idea of energy-abundant/super-Hinze (
$We\gg We_{crit}$ and 
$D\gg D_H$) versus energy-limited/sub-Hinze (
$We< We_{crit}$ and 
$D< D_H$) breakups. The former has been studied much more extensively than the latter for a simple reason: breakup of super-Hinze bubbles is much faster and more frequent so it is easier to observe in a finite volume and to collect enough statistics. Breakup of super-Hinze bubbles is typically studied in several different flow configurations: pipe flow (Hesketh, Etchells & Russell Reference Hesketh, Etchells and Russell1991) and turbulent jets (Sevik & Park Reference Sevik and Park1973; Martínez-Bazán, Montanes & Lasheras Reference Martínez-Bazán, Montanes and Lasheras1999b; Vejražka, Zedníková & Stanovskỳ Reference Vejražka, Zedníková and Stanovskỳ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 Reference Hulburt and Katz1964; Ramkrishna Reference Ramkrishna2000). However, this framework has one limitation: it 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 Reference Deane and Stokes2002), 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 Reference Hinze1955); bubbles tend to exhibit regular affine deformation in these types of flows (lenticular or cigar-shaped); (ii) a resonance mechanism that relies on bubble oscillation to siphon energy from the surrounding turbulence until breakup (Sevik & Park Reference Sevik and Park1973; Hesketh et al. Reference Hesketh, Etchells and Russell1991; Risso & Fabre Reference Risso and Fabre1998); it typically assumes that the surrounding eddies retain a similar frequency with bubbles’ natural frequency; (iii) an inertial mechanism that 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 Reference Villermaux and Bossa2009) and bag breakup in crossflows (Ng, Sankarakrishnan & Sallam Reference Ng, Sankarakrishnan and Sallam2008). In turbulence, the three mechanisms may all be present, so applying only one mean Weber number to account for all breakup mechanisms is questionable.
In addition, as Risso & Fabre (Reference Risso and Fabre1998) noted, the instantaneous and local Weber number, 
$We$, could be much larger than the mean value 
$\langle We\rangle$. They proposed to use the time trace of 
$We$ 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 two or three dimensions, remain challenging.
To resolve deformation and breakup of the Hinze- or sub-Hinze-scale bubbles, in this paper, we introduce an experiment that provides unique simultaneous measurements of both bubble deformation and surrounding flows thanks to the recent advancement of the three-dimensional (3-D) high-concentration particle shadow tracking (Tan et al. Reference Tan, Salibindla, Masuk and Ni2019) and 3-D virtual-camera visual-hull shape reconstruction (Masuk, Salibindla & Ni Reference Masuk, Salibindla and Ni2019a). In § 2, the experimental set-up, i.e. a vertical water tunnel system with a large section of homogeneous and isotropic turbulence, is introduced. In the same section, the optical system designed to conduct simultaneous measurements of both the phases is also discussed. In § 4, based on the new datasets, we discuss how flow decomposition can be conducted to analyse the relative roles played by different mechanisms. In § 4.5, we finally estimate the breakup probability of Hinze-scale bubbles in turbulence.
2. Experimental set-up
As shown in figure 1, a facility named V-ONSET 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 a low energy dissipation rate (Variano, Bodenschatz & Cowen Reference Variano, Bodenschatz and Cowen2004; Mercado et al. Reference Mercado, Prakash, Tagawa, Sun and Lohse2012) (
$\langle \epsilon \rangle =O(10^{-5}\text {--}10^{-3})\ \textrm {m}^{2}\ \textrm {s}^{-3}$), whereas facilities that use water jets to break super-Hinze-scale bubbles can generate a large energy dissipation rate 
$\langle \epsilon \rangle =O(0.1\text {--}10^{3})\ \textrm {m}^{2}\ \textrm {s}^{-3}$ at the cost of having strong flow inhomogeneity and anisotropy (Martínez-Bazán et al. Reference Martínez-Bazán, Montanes and Lasheras1999b; Vejražka et al. Reference Vejražka, Zedníková and Stanovskỳ2018).
Figure 1. Schematic of the V-ONSET vertical water tunnel; two insets show the 3-D 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. (Reference Masuk, Salibindla, Tan and Ni2019b).
The experimental set-up used for the current study is essentially a vertical water tunnel capable of generating turbulence with 
$\langle \epsilon \rangle$ roughly at 
$0.16\text {--}0.5\ \textrm {m}^{2}\ \textrm {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 
$\langle \epsilon \rangle$ in this region, 
$\langle We \rangle$ 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 (with jet velocity up to 
$12\ \textrm {m}\ \textrm {s}^{-1}$), each of which has a diameter 
$d_j$ of 5 mm, firing coaxially downward into the test section along with the mean flow. The firing pattern of these momentum jets was randomised in a way similar to the work by Variano et al. (Reference Variano, Bodenschatz and Cowen2004) in order to ensure that no secondary flow structure would develop in the test section (De Silva & Fernando Reference De Silva and Fernando1994; Srdic, Fernando & Montenegro Reference Srdic, Fernando and Montenegro1996; Variano et al. Reference Variano, Bodenschatz and Cowen2004). On average, 12.5 % of the jets were kept on at a time as this was found to maximise the turbulence intensity. The test section was set much farther downstream of the jets (about 80
$d_j$) to ensure that the jets were well mixed and turbulence becomes homogeneous and isotropic with very little spatial variation. Additional details concerning this set-up and its flow characteristics can be found in Masuk et al. (Reference Masuk, Salibindla, Tan and Ni2019b).
Bubbles were generated at the bottom of the test section using two different sizes of hypodermic needles (small needles with inner diameter (ID) of 
$160\ \mathrm {\mu }\textrm {m}$ and outer diameter (OD) of 
$300\ \mathrm {\mu }\textrm {m}$; large needles with ID of 
$260\ \mathrm {\mu }\textrm {m}$ and OD of 
$500\ \mathrm {\mu }\textrm {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 high-speed cameras each with a 
$1024\times 1024$ pixel resolution and 4000 frame per second (fps) frame rate. The frame rate was selected to ensure that each particle could be imaged about 10 times within one Kolmogorov timescale 
$\tau _\eta =2.5\ \textrm {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.
3. Flow characterisation
Before discussing bubble deformation and breakup, single-phase turbulence statistics in the tunnel needs to be characterised 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. (Reference Masuk, Salibindla, Tan and Ni2019b). 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\eta$, with 
$\eta \approx 50\ \mathrm {\mu }\textrm {m}$ being the Kolmogorov length scale 
$\eta =(\nu ^{3}/\langle \epsilon \rangle )^{1/4}$. Here 
$\nu$ 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. Reference Tan, Salibindla, Masuk and Ni2020) that employs the Shake-the-Box method (Schanz, Gesemann & Schröder Reference Schanz, Gesemann and Schröder2016).
Figure 2. The longitudinal structure function 
$D_{LL}$ as a function of the scale separation 
$r$ normalised by the Kolmogorov length scale 
$\eta$. The dashed and solid lines indicate the dissipative and inertial range scalings based on the Kolmogorov theory, respectively.
The structure function should approach two limits: one in the dissipative range (
$r\ll \eta$) and the other in the inertial range (
$\eta \ll r\ll L$). 
$L$ is the integral scale, which is estimated based on 
$L\approx u'^{3}/\langle \epsilon \rangle$, where 
$u'$ is the fluctuation velocity. The scale separation between 
$\eta$ and 
$L$ is determined by the Taylor-scale Reynolds number 
$Re_{\lambda }=\sqrt {15 u'L/\nu }$, which is estimated to be around 435. In the dissipative range, the structure function follows the relationship of 
$D_{LL}=(\epsilon /15\nu )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 & Xia Reference Ni and Xia2013), using a standard number 
$C_2=2.13$ can provide a reasonable estimation of 
$\epsilon$. The solid line shown in the figure is based on the calculated 
$\langle \epsilon \rangle =0.16\ \textrm {m}^{2}\ \textrm {s}^{-3}$. However, if 
$\langle \epsilon \rangle$ 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 
$\langle \epsilon \rangle$ estimated from either the dissipative or inertial range helps to quantify the experimental uncertainty of the mean energy dissipation rate: 
$\langle \epsilon \rangle =0.22\pm 0.07\ \textrm {m}^{2}\ \textrm {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\ \textrm {m}^{2}\ \textrm {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 
$\eta$. 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.
4. Results and discussion
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. Reference Masuk, Salibindla and Ni2019a) for 3-D shape reconstruction. Averaging surface points on the reconstructed geometry helps to determine the centre of mass, which was then tracked in three dimensions 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 3-D 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 s (360 frames) before exiting.
Figure 3. (a) Raw images of one highly deformed bubble observed by six high-speed cameras simultaneously. (b) The outline and silhouette of the same bubble extracted from camera 4. (c) 3-D tracks of about 40 tracer particles within 
$4D$ (
$D$ is the volume-equivalent sphere diameter) from the centre of the bubble that is shown as a 3-D reconstructed geometry.
Separated images for tracer particles were input into our in-house OpenLPT (Tan et al. Reference Tan, Salibindla, Masuk and Ni2019) 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 6000 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 3-D shape reconstructed from silhouettes segmented from figure 3(a). In this case, trajectories of tracer particles within 
$4D$ away from the bubble centre were included. These tracks could be used to estimate the flow condition around the bubble. As 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, turbulence energy dissipation rate and the full coarse-grained velocity gradients, locally, instantaneously and along each bubble trajectory. Additional information concerning the set-up and measurement techniques can be found in works by Masuk et al. (Reference Masuk, Salibindla and Ni2019a,b) and Tan et al. (Reference Tan, Salibindla, Masuk and Ni2019).
4.2. Flow velocity and velocity gradient
For a bubble at location 
$\boldsymbol {x}_{\boldsymbol {0}}$, its surrounding flow velocity 
$\boldsymbol {u}^{p}(\boldsymbol {x}_{\boldsymbol {0}}+\boldsymbol {x}^{p})$ can be measured at a number of discrete positions 
$\boldsymbol {x}_{\boldsymbol {0}}+\boldsymbol {x}^{p}$ where 
$n$ tracer particles are located (
$p=1,2,\ldots ,n$). These tracer particles are sought within a radius of 
$D_s/2$ from the bubble centre 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:
\begin{equation} \left.\begin{gathered} u^{p}_i(\boldsymbol{x}_{\boldsymbol{0}}+\boldsymbol{x}^{p}) \approx \overline{u_i}(\boldsymbol{x}_{\boldsymbol{0}}) + \skew4\tilde{{\mathsf{A}}}_{ij}(\boldsymbol{x}_{\boldsymbol{0}})x^{p}_j + O\left(x^{p}_j \skew4\tilde{{\mathsf{H}}}_{jik} (\boldsymbol{x}_{\boldsymbol{0}})x^{p}_k\right),\\ \skew4\tilde{{\mathsf{A}}}_{ij}(\boldsymbol{x}_{\boldsymbol{0}})=\frac{\partial u^{p}_i}{\partial x^{p}_j}\quad \text{and}\quad \skew4\tilde{{\mathsf{H}}}_{jik}(\boldsymbol{x}_{\boldsymbol{0}})=\frac{\partial^{2} u^{p}_i}{\partial x^{p}_j\partial x^{p}_k}, \end{gathered}\right\} \end{equation}
where 
$\overline {u_i}=\sum ^{n}_{p=1} u^{p}_i(\boldsymbol {x}_{\boldsymbol {0}}+\boldsymbol {x}^{p})/n$ represents the local mean flow estimated by averaging the velocity of 
$n$ tracer particles. Here 
$\skew4\tilde{{\mathsf{A}}}_{ij}$ and 
$\skew4\tilde{{\mathsf{H}}}_{jik}$ indicate the velocity gradient tensor and the Hessian matrix, respectively, and the tilde denotes coarse graining at the bubble size. We use 
$\boldsymbol {x}^{p}$ to denote the separation vector directed from the bubble centre at 
$\boldsymbol {x}_{\boldsymbol {0}}$ to the 
$p$th tracer particle location. For small micro-bubbles with sizes in the dissipative range (
$D\ll \eta$), 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 (
$\eta \ll D\ll L$).
Although the velocity gradient around each bubble can be measured accurately (Ni et al. Reference Ni, Kramel, Ouellette and Voth2015), the velocity Hessian, on the other hand, requires measuring the gradient of the velocity gradient (three 
$3\times 3$ matrices). Even though it is possible to calculate the velocity Hessian given sufficient number of tracer particles, the uncertainty becomes large owing to the second-order spatial derivative. As a result, we limit only to the first two orders, i.e. the mean flow velocity 
$\overline {u_i}$ and the velocity gradient 
$\skew4\tilde{{\mathsf{A}}}_{ij}(\boldsymbol {x}_{\boldsymbol {0}})$, to capture the key mechanisms of deformation.
The velocity gradient tensor 
$\skew4\tilde{{\mathsf{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 
$\sum _p[u^{p}_i-\skew4\tilde{{\mathsf{A}}}_{ij}x^{p}_j]^{2}$ (Pumir, Bodenschatz & Xu Reference Pumir, Bodenschatz and Xu2013; Ni et al. Reference Ni, Kramel, Ouellette and Voth2015). Although finite-sized bubbles typically come with a large search radius and abundant nearby tracer particles thanks to our tracking method (Tan et al. Reference Tan, Salibindla, Masuk and Ni2020), particles in the vicinity of a bubble are still distributed randomly in space. If nearby particles stay primarily within a quasi-2-D plane, the estimation of the out-of-plane velocity gradient will have large uncertainty. Similar to previous studies (Xu, Pumir & Bodenschatz Reference Xu, Pumir and Bodenschatz2011; Ni et al. Reference Ni, Kramel, Ouellette and Voth2015), an inertia tensor 
$I=\sum _px^{p}_ix^{p}_j/\textrm {tr}(\sum _px^{p}_ix^{p}_j)$ was adopted to evaluate the shape factor of the particle cloud. If particles are uniformly distributed in three dimensions, three eigenvalues of this inertia tensor (
$\gamma _i$) all equal to 1/3. For a quasi-2-D distribution, the smallest eigenvalue (
$\gamma _3$) will be very close to zero, and the gradient along that direction cannot be calculated. In practice, events with 
$\gamma _3/\gamma _1$ smaller than 0.15 was therefore removed from the statistics.
Based on 
$\skew4\tilde{{\mathsf{A}}}_{ij}$, the coarse-grained rate-of-strain tensor, 
$\skew4\tilde{{\mathsf{S}}}_{ij}$, and rotation tensor, 
$\tilde{\Omega}_{ij}$, can be obtained directly: 
$\skew4\tilde{{\mathsf{S}}}_{ij}=(\skew4\tilde{{\mathsf{A}}}_{ij}+\skew4\tilde {{\mathsf{A}}}_{ji})/{2}$, 
$\tilde{\Omega}_{ij}=(\skew4\tilde{{\mathsf{A}}}_{ij}-\skew4\tilde {{\mathsf{A}}}_{ji})/{2}$. Figure 4(a) shows the probability density function (PDF) of two eigenvalues of 
$\skew4\tilde{{\mathsf{S}}}_{ij}$ (the largest 
$\widetilde {\lambda _1}$ and the smallest 
$\widetilde {\lambda _3}$) based on different 
$D_s$. The PDFs of 
$|\widetilde {\lambda _1}|$ and 
$|\widetilde {\lambda _3}|$ overlap 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 
$\skew4\tilde{{\mathsf{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 
$\skew4\tilde{{\mathsf{A}}}_{ij}$ can be related to the local energy dissipation rate in the form of 
$C_2(\tilde {\epsilon } d)^{2/3}=(\widetilde {\lambda _3}d)^{2}$, where 
$C_2=2.13$ is the Kolmogorov constant (Batchelor Reference Batchelor1953; Sreenivasan Reference Sreenivasan1995; Ni & Xia Reference Ni and Xia2013) and 
$\tilde {\epsilon }$ is the coarse-grained energy dissipation rate for a range of length scales 
$d$ considered. To check whether the distribution of 
$\tilde {\epsilon }$ is indeed the same for 
$d$ varying between 
$D$ to 
$D_s$, in figure 4(b), the PDFs of the estimated local 
$\tilde {\epsilon }$ using three different 
$D_s$ are shown. Although the distribution of the calculated 
$\widetilde {\lambda _3}$ are sensitive to 
$D_s$, once converted to 
$\tilde {\epsilon }$, three curves from all three 
$D_s$ fall right on top of each other, indicating that the local 
$\tilde {\epsilon }$ 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.
Figure 4. (a) The distribution of the two eigenvalues (
$\widetilde {\lambda _1}$ and 
$\widetilde {\lambda _3}$) of the local rate-of-strain tensor coarse grained at the bubble scale 
$D$ (
$|\tilde {\lambda }|$ is used here because 
$\widetilde {\lambda _3}<0$). Three search diameters ranging from 2–4
$D$ to 6–8
$D$ are denoted by different colours. (b) The distribution of the local coarse-grained energy dissipation rate 
$\tilde {\epsilon }$. The log-normal distribution from (4.2) is shown as the black solid line.
The coarse-grained energy dissipation rate can be described by the log-normal distribution based on the Kolmogorov refined theory in 1962 (Kolmogorov Reference Kolmogorov1962) and multi-fractal spectrum (Meneveau & Sreenivasan Reference Meneveau and Sreenivasan1991):
\begin{align} P\left(\frac{\epsilon_r}{\langle\epsilon\rangle}\right) = \frac{1}{\epsilon_r/\langle\epsilon\rangle}\frac{1}{\sqrt{2{\rm \pi} (A+\mu \ln(L/r))}} \exp\left[-\frac{(\ln(\epsilon_r/\langle\epsilon\rangle) + 1/2(A+\mu \ln(L/r)))^{2}}{2(A+\mu \ln(L/r))}\right], \end{align}
where 
$\epsilon _r$ is the coarse-grained energy dissipation rate at a scale 
$r$. Here 
$\mu \approx 0.25$ is the intermittency exponent, 
$L=3.2\text {--}6\ \textrm {cm}$ is the integral length scale and 
$A$ is a parameter that needs to be fitted to the experimental data to determine the variance of 
$\tilde {\epsilon }$ when 
$r=L$, which was found to be around one. Based on the definition, 
$\tilde {\epsilon }$ measured from our experiments is equivalent to 
$\epsilon _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 (4.2) shows that the measured coarse-grained energy dissipation rate is consistent with the classical Kolmogorov theory.
4.3. Different types of deformation
4.3.1. Bubble deformation by the velocity gradient 
$We_{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, 
$\widetilde {\lambda _3}$, which is associated with the direction that compresses the most, should be the more relevant eigenvalue of 
$\skew4\tilde{{\mathsf{A}}}_{ij}$. Following the argument, the Weber number can be defined as
\begin{equation} We_{vg}=\frac{\rho(\widetilde{\lambda_3}D)^{2}D}{\sigma}\sim\frac{C_2(\tilde{\epsilon} D)^{2/3}}{\sigma}. \end{equation} This Weber number definition is based on the local coarse-grained 
$\skew4\tilde{{\mathsf{A}}}_{ij}$ and 
$\tilde {\epsilon }$, which is different from the mean Weber number defined by Kolmogorov (Reference Kolmogorov1949) and Hinze (Reference Hinze1955). Figure 5 shows the distribution of local 
$We_{vg}$ based on the measurements of 
$\skew4\tilde{{\mathsf{A}}}_{ij}$ along each bubble track. The distribution peaks at around one, suggesting that those bubbles are indeed Hinze-scale bubbles. All data points on the right-hand side of the peak represent bubbles deforming under strong velocity gradients. On top of the experimental results, the model of 
$We_{vg}$ based on (4.2) and (4.3) is also shown. Similar to figure 4(b), the log-normal distribution of the local 
$\tilde {\epsilon }$ explains the observed shape of the PDF of 
$We_{vg}$, from which bubble breakup probability can be determined.
Figure 5. The distribution of the measured Weber numbers, based on the slip velocity, i.e. 
$We_{slip,x}$ (blue circle) and 
$We_{slip,z}$ (black plus), and the velocity gradient, 
$We_{vg}$ (red triangle). Two solid lines represent the modelled Weber number distributions based on the log-normal distribution of 
$\tilde {\epsilon }$ ((4.2) and (4.3), red line) and the stretched-exponential fit of the slip velocity ((4.4) and (4.5), blue line), respectively.
4.3.2. Slip-velocity induced deformation $We_{slip}$
As Hinze stated in his original seminal work (Hinze Reference Hinze1955), 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 
$We_{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, 
$\boldsymbol {u}_{slip} = \boldsymbol {u}_b - \boldsymbol {u}_f$. As its name suggests, 
$\boldsymbol {u}_{slip}$ quantifies the drift of the bubble velocity 
$\boldsymbol {u}_b$ away from the instantaneous local flow velocity 
$\boldsymbol {u}_f$.
Here 
$\boldsymbol {u}_f$ represents the continuous-phase fluid velocity at the centre of a bubble if the bubble was not there. In practice, 
$\boldsymbol {u}_f$ needs to be estimated from the continuous-phase velocities measured in the vicinity of the bubble. Therefore, we assume 
$\boldsymbol {u}_f$ to be the same as 
$\overline {u_i}$ in (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 
$\boldsymbol {u}_f$ normalised by its own standard deviation. 
$\boldsymbol {u}_f$ can be calculated around bubbles of different sizes, which are shown by different coloured 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. As 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\text {--}4D$ to 6–8
$D$. As shown in figure 6(b), no discernible difference is observed for the flow velocity PDF at three 
$D_s$, which suggests that 
$\boldsymbol {u}_f$ is not sensitive to 
$D_s$.
Figure 6. The distribution of the horizontal flow velocity 
$u_{f,x}$ (normalised by its own standard deviation) nearby bubbles of (a) different sizes 
$D$ and (b) different search diameters 
$D_s$. The black solid lines in (a,b) show the standard normal distribution for reference.
We use 
$\boldsymbol {u}_b$ to denote 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}$, normalised 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 
$\boldsymbol {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 (Ern et al. Reference Ern, Risso, Fabre and Magnaudet2012), 
$\langle u_{b,x}^{2}\rangle ^{1/2}$ should depend on the bubble size. Therefore, the observed nearly constant 
$\langle u_{b,x}^{2}\rangle ^{1/2}$ clearly indicates that the buoyancy effect is negligible in the horizontal directions owing to the background intense turbulence. Figure 7(b) also displays the standard deviation of 
$\boldsymbol {u}_f$ along two horizontal directions. In contrast to 
$\langle u_{b,x}^{2}\rangle ^{1/2}$ , 
$\langle u_{f,x}^{2}\rangle ^{1/2}$ seems to decrease as 
$D$ increases because a finite-sized bubble effectively serves as a filter that reduces the local flow fluctuations. By extrapolating both 
$\langle u_{b,x}^{2}\rangle ^{1/2}$ and 
$\langle u_{f,x}^{2}\rangle ^{1/2}$ to small bubble sizes, 
$\langle u_{b,x}^{2}\rangle ^{1/2}$ and 
$\langle u_{f,x}^{2}\rangle ^{1/2}$ will eventually cross at around 
$200\ \textrm{mm s}^{-1}$ for bubble size close to zero, which gives the right limit as extremely small bubbles should behave similarly to tracers 
$\langle u_{b,x}^{2}\rangle ^{1/2}\approx \langle u_{f,x}^{2}\rangle ^{1/2}$.
Figure 7. (a) The distribution of the horizontal velocity 
$u_{b,x}$ of bubbles (normalised by its own standard deviation) with different diameter 
$D$. The black solid line indicates the standard normal distribution for comparison. (b) The fluctuation of mean flow velocity (solid) and bubble velocity (dashed lines) along two different directions versus bubble size 
$D$.
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 normalised 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, Sreenivasan & Stolovitzky Reference Kailasnath, Sreenivasan and Stolovitzky1992; Sreenivasan Reference Sreenivasan1999; Li & Meneveau Reference Li and Meneveau2005). The PDF of the velocity increment has been fitted with a stretched exponential function (Kailasnath et al. Reference Kailasnath, Sreenivasan and Stolovitzky1992), which is adopted here to describe the observed PDF of the slip velocity.
\begin{gather} P(u_{slip,x})=C \exp\left[-Q\left( \frac{u_{slip,x}}{\langle u^{2}_{slip,x} \rangle^{1/2}}\right)^{m}\right] , \end{gather}
\begin{gather}We_{slip} = \frac{\rho u_{slip}^{2} D}{\sigma}, \end{gather}
where 
$C$ is the normalisation 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.03
$L$ to 0.12
$L$ 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. (Reference Kailasnath, Sreenivasan and Stolovitzky1992). 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 normalised 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.
Figure 8. (a) The distribution of the normalised horizontal slip velocity between the two phases. Symbols denote bubbles of different sizes and the black solid line indicates the standard normal distribution. The red solid line shows the stretched exponential (4.4) fit to the data. (b) The fluctuation slip velocity of all three components versus the bubble diameter 
$D$; The black dash-dotted line indicates the estimation from the second-order structure function. The prefactor 
$4/9$ is chosen to minimise the offset between the solid line and the data.
Furthermore, the fluctuation slip velocity (
$\langle u_{slip}^{2}\rangle ^{1/2}$) increases as a function of bubble size 
$D$, suggesting that larger bubbles with 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 
$(\langle \epsilon \rangle D)^{1/3}$. After assuming that these two velocity scales are related, the measured 
$\langle u_{slip}^{2}\rangle ^{1/2}$ along all three directions are fitted with 
$\gamma (\langle \epsilon \rangle D)^{1/3}$ by performing the least-square regression to obtain the fitting coefficient 
$\gamma$, 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 
$\langle u_{slip}^{2}\rangle ^{1/2}$ for bubbles with size in the inertial range. With this relationship and two coefficients, i.e. 
$Q=3/4$ and 
$m=6/5$, the distribution of 
$u_{slip}$ can be estimated from (4.4).
Finally, the distribution of 
$We_{slip}$, calculated based on (4.4) and (4.5), is shown in figure 5. The blue solid line indicates the predicted 
$We_{slip}$ based on the stretched exponential fit to the horizontal slip velocity 
$u_{slip,x}$ (4.4). The distribution also peaks at around 
$We\approx 1$, which is slightly smaller than the most probable value of 
$We_{vg}$. The right tails of both PDFs (
$We_{vg}$ and 
$We_{slip,x}$), corresponding to the range of 
$We$ 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.
4.3.3. 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 
$We_{slip}$ in the vertical direction based on the 
$z$-axis slip velocity, i.e. 
$We_{slip,z}$ is also shown. This PDF has a bump near 
$We_{slip,z}\approx 3\text {--}4$ because of the buoyancy effect, but both the left and right tails seem to agree with those of 
$We_{slip,x}$. This suggests that both the turbulence effect and the buoyancy effect are present in 
$We_{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 
$We_{slip,z}$ is similar to the Eötvös number, 
$Eo=\rho gD^{2}/\sigma$, as the terminal vertical slip velocity 
$u_{slip,z}$ driven primarily by buoyancy should be proportional to 
$\sqrt {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 
$\epsilon$. In a recent paper (Salibindla et al. Reference Salibindla, Masuk, Tan and Ni2020), 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_b^{0.687}),{{\min}}(\,f(Eo), f(Eo)/We^{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} = 2 V_b (\rho - \rho _b) g / \rho A C_D$, where 
$V_b$ is the volume of a bubble. 
$A$ is the projected area of volume-equivalent spherical bubble and 
$\rho _b$ is the density of bubbles. For the bubble size range considered, 
$We_{slip,z}$ calculated based on 
$C_D$ is about 4, which is consistent with the bump of 
$We_{slip,z}$ observed in the PDF. This agreement confirms that the observed bump in the distribution of 
$We_{slip,z}$ indeed comes from the buoyancy-induced bubble rise velocity.
All together, it seems that the bump in the distribution of 
$We_{slip,z}$ is limited to a narrow range, and the right tail of 
$We_{slip,z}$ seems to be close to that of 
$We_{slip,x}$ and 
$We_{vg}$. This suggests that, at least for our parameters when 
$\langle \epsilon \rangle \approx 0.2\text {--}0.5 \ \textrm {m}^{2}\ \textrm {s}^{-3}$, the buoyancy-induced deformation is limited. If we keep increasing 
$\langle \epsilon \rangle$, the buoyancy effect will become even weaker.
4.4. Bubble aspect ratio versus 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 are used to study the mechanisms of bubble deformation and breakup in turbulence.
4.4.1. Simultaneous measurements of bubble geometry and $We$
Figure 9 shows two examples of the simultaneous measurements of bubble aspect ratio (
$\alpha$) and the Weber numbers (
$We_{vg}$ and 
$We_{slip}$). Here 
$\alpha$ was calculated as 
$\alpha =r_1/r_3$ where the 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-centre distances from the 3-D-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 five. The time trace of 
$\alpha$ is not similar to that of either 
$We$ at first glance. The only evident correlation is probably at 
$t=10\text {--}20\ \textrm {ms}$, when a small peak observed in the time trace of 
$\alpha$ seems to be driven by a large 
$We_{slip}$ maybe 5 ms earlier. For 
$t=50\text {--}70\ \textrm {ms}$, despite 
$We_{slip}$ drops close to zero, 
$\alpha$ continues to rise thanks to a relatively large value of 
$We_{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 (We_{vg},We_{slip})$.
Figure 9. (a,b) Two example time traces of bubble aspect ratio (black circles) and their corresponding Weber numbers (c,d) during deformation. Each measured Weber number time trace (
$We_{vg}$ (blue squares) and 
$We_{slip}$ (red triangles)) in panels (c) and (d) can be input into (4.6) to obtain a calculated time trace of 
$\alpha ^{*}$, which are shown as the dashed lines of the same colour in panels (a) and (b). In addition, the solid black line represents the calculated results by using 
$\max$(
$We_{vg}$,
$We_{slip}$) as an input. (Note that the blue dashed line in (b) is not visible because it overlaps perfectly with the black line.)
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 (Reference Risso and Fabre1998) and Lalanne, Masbernat & Risso (Reference Lalanne, Masbernat and Risso2019), which is essentially a forced oscillator model that connects bubble deformation directly to local 
$We$ 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
\begin{equation} \frac{\textrm{d}^{2}\hat{a}}{\textrm{d}\hat{t}^{2}}+2\xi\frac{\textrm{d}\hat{a}}{\textrm{d}\hat{t}}+\hat{a}=K'We(t), \end{equation}
where 
$\xi =1/2{\rm \pi} \tau _d\, f_2$ is the damping coefficient, 
$\tau _d = D^{2} / 80 \nu$ is the damping time scale (Risso & Fabre Reference Risso and Fabre1998) and 
$f_2 = \sqrt {96\sigma / \rho D^{3}}$ is the mode 2 natural frequency of bubble oscillation (Lamb Reference Lamb1932).
We recognise 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. Linearising this problem relies on two assumptions: (i) the Weber numbers are not very large (
${\lesssim }O(1$)); and (ii) bubble deformation is driven primarily by eddies of the bubble size. For assumption (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. As there is no consensus on the value of the critical Weber number, 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 
$\eta$ 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 
$We\lesssim 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. However, 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 assumption (ii), the Weber numbers defined based on flows of the bubble size essentially ignore the contributions from sub-bubble-scale eddies. Statistically, this assumption probably holds because smaller eddies tend to be weaker, even though they could occasionally become exceedingly strong owing to turbulent intermittency. However, it would require further investigations to understand their contributions.
In (4.6), the amplitude of the instantaneous 
$We(t)$ is controlled by the prefactor 
$K'$. Note that 
$We(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 Reference Risso and Fabre1998). In this work, in addition to measuring 
$We(t)$ directly, the method also allows us to distinguish between 
$We_{vg}$ and 
$We_{slip}$. However, because the model did not explicitly account for individual 
$We$, here we apply three different inputs, 
$We(t)=We_{vg}$, 
$We(t)=We_{slip}$ and 
$We(t)={\max}(We_{vg},We_{slip})$, to (4.6) to obtain model predictions, which are shown in figures 9 and 10 as a blue dashed line, red dashed line and black solid line, respectively.
Figure 10. (a,b) Example time traces of two breakup events. Symbols and lines are the same as those in figure 9. The bubble breaks at 89.5 ms for (a) and at 56.5 ms for (b).
Note that 
$\hat {a}$ in (4.6) is defined as the ratio of the deformed radius (
$a=r_1-D/2$) to 
$D$, where 
$D$ is the diameter of an volume-equivalent sphere. It has to be converted into 
$\alpha ^{*}=2(\hat {a}D+D/2)/D$ for comparisons with the measured 
$\alpha$. However, despite our best efforts, 
$\alpha \neq \alpha ^{*}$ because 
$\alpha ^{*}$ does not contain information about the minor axis, which has to be replaced with 
$D/2$. Nevertheless, 
$\alpha$ and 
$\alpha ^{*}$ should share the same trend as bubbles deform.
Indeed, similarities can be observed between 
$\alpha$ and 
$\alpha ^{*}$ in figure 9(a). The magnitude of 
$\alpha ^{*}$ is affected by the prefactor 
$K'$ in (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 
$We_{vg}$ and 
$We_{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 keeps changing in the experimental data, e.g.the second peak is much closer to the first and further away from the third. This varying phase lag is a feature that cannot be reproduced from the model.
Figure 9(b) shows another example to compare 
$\alpha$ with 
$\alpha ^{*}$. A large 
$\alpha$ observed at 32 ms seems to correlate with an event of large 
$We_{vg}$ occurred at 20 ms, whereas a small bump of 
$\alpha$ at 60 ms seems to correspond to a sudden increase of 
$We_{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, the three model-predicted time traces of 
$\alpha ^{*}$ differ; 
$We_{vg}$ is systematically larger than 
$We_{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 linear oscillator model provided by Risso & Fabre (Reference Risso and Fabre1998) may capture the dynamics of bubbles undergoing small-amplitude deformation (
$\alpha \approx 1\text {--}3$) with small 
$We$ (
$We\lesssim 3$), as in figure 9(a), but not for all bubbles, particularly not for bubbles with large 
$We$ 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.
As the model provided by Risso & Fabre (Reference Risso and Fabre1998) is primarily designed to characterise the breakup process, figure 10 shows two examples of bubbles that eventually break. The Weber number 
$We$ for breaking bubbles is clearly much larger: one reaches close to 40 and the other climbs up to almost 20, nearly a factor of 4–8 larger than the cases for small-amplitude deformation in figure 9(a). In figure 10(a), 
$We_{slip}$ dominates, but a local event at 
$We_{slip}\approx 40$ did not break the bubble, even though it did successfully deform the bubble to a large 
$\alpha$ at around 8. Following a peak of 
$We_{slip}$ at 85 ms, this bubble eventually broke at 
$t=90$ ms, with the instantaneous 
$We$ about 10 and local aspect ratio close to 2.
This example represents many bubbles we observed that do not break at the moment when 
$\alpha$ reaches its peak; instead, they split at a later time during the process of retraction. As a bubble retracts back to 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 Reference Dodd and Ferrante2016). However, this process is unstable because of the large density difference between the two phases, and it eventually leads to breakup.
For the first example (figure 10a), 
$We_{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 Reference Prince and Blanch1990; Risso & Fabre Reference Risso and Fabre1998) 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 the 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 
$\alpha ^{*}$ 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 3-D 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.
4.4.2. Distribution of bubble aspect ratio
In addition to the response of individual bubbles to different Weber numbers, the distribution of 
$\alpha$ could also be connected to that of 
$We$ to examine whether the bubble aspect ratio can be solely determined by 
$We$ in a statistical sense. This relationship between 
$\alpha$ and 
$We$ was first derived by Moore (Reference Moore1965) for bubbles rising in a quiescent medium:
\begin{equation} \alpha=1+\tfrac{9}{64}{We}+O(We^{2}), \end{equation}
where 
$We$ 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 
$We\ll 1$ so that the departure from a spherical shape is so small that any high-order terms associated with 
$We^{2}$ can be ignored.
For 
$We\approx 1$ or above, the potential flow theory applied to oblate ellipsoids with fore–aft symmetry yields
\begin{equation} We(\alpha) = 4\alpha^{-4/3}(\alpha^{3}+\alpha-2)[\alpha^{2}\sec^{-1}\alpha-(\alpha^{2}-1)^{1/2}]^{2}(\alpha^{2}-1)^{-3}. \end{equation} 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 
$\langle We\rangle \approx 1$, it cannot predict the relationship between 
$\alpha$ and instantaneous 
$We$ for 
$We>3.745$, which is about 21.5 % of the total events in our experiments.
Figure 11(a) shows the PDF of 
$\alpha$ for bubble size 
$D=4.5$ mm. The PDF peaks at 
$\alpha =1.7$ and has a long tail that skews towards larger values of 
$\alpha$. On top of the experimental results, the PDF of 
$\alpha$ calculated from (4.7) using 
$We=We_{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. As the long tail in the PDF of 
$We_{vg}$ comes from the log-normal distribution of 
$\tilde {\epsilon }$, a large probability of strong deformation 
$\alpha >2$ is likely to be contributed by intermittent events with a large 
$\tilde {\epsilon }$. If 
$\alpha$ versus 
$We_{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 (4.7), which did not provide satisfactory results (not shown here). Finally, after adjusting the parameters in (4.7), we settle down to a simple new equation
\begin{equation} \alpha = \tfrac{2}{5} We^{{2}/{3}} + 1.2 \end{equation}
to fit the data, which is shown as the red solid line in figure 11(a). This new fit shows an excellent agreement with the measured PDF of 
$\alpha$. Moreover. the solid blue line in figure 11(a) shows the PDF of 
$\alpha$ by implementing a different Weber number 
$We=We_{slip}$ in (4.9). Similar trend of 
$\alpha$ can be seen even with this Weber number. However, compared with 
$We_{vg}$, the results based on 
$We_{slip}$ tend to underpredict 
$\alpha$, which is consistent with the observation in figure 5 that the peak of the PDF of 
$We_{slip}$ is on the left-hand side of the PDF of 
$We_{vg}$. Nevertheless, the right tail of 
$\alpha$ can be reproduced by the calculations using both 
$We_{slip}$ and 
$We_{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, are more important for mild deformation.
Figure 11. The distribution of bubble aspect ratio 
$\alpha$ for (a) one size 
$D=4.5$ mm to test against different 
$We$ versus 
$\alpha$ relationship listed in two different equations (4.7) and (4.8) by using either 
$We_{slip}$ (blue) and 
$We_{vg}$ (red) and for (b) a range of sizes from 2.5 to 6.5 mm; solid lines are calculated from (4.9).
Figure 11(b) compiles the PDFs of 
$\alpha$ for five different bubble sizes, from 2.5 to 6.5 mm with an interval of 1 mm. The peaks of these PDFs progressively shift rightward towards larger 
$\alpha$ as 
$D$ grows, which is consistent with our intuition that large bubbles are more deformable and, thus, have a larger 
$\alpha$ 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 
$\alpha$ much larger than the mean also increases. To explain this observation, (4.9) is applied to all these cases with different 
$D$, and results are shown as solid lines with corresponding colours to compare with the PDF of measured 
$\alpha$. The modeled PDF of 
$\alpha$ 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 
$\alpha$ as a function of 
$D$ is driven mostly by the change of 
$We$, but the relationship between 
$\alpha$ and 
$We$ may not be linear for bubbles deformed by turbulence.
An alternative method to predict the relationship between 
$\alpha$ and 
$We$ is to use the model provided by (4.6). As discussed before, although the model-predicted time trace of 
$\alpha$ 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 and 5.5 mm for simplicity. The model seems to reproduce the overall trend of the PDF. However, the PDF of 
$\alpha ^{*}$ calculated from the model is flatter than that of the measured 
$\alpha$, 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.
Figure 12. The distribution of bubble aspect ratio 
$\alpha$ for two different sizes (
$D=2\text {--}3$ and 
$D=5\text {--}6\ \textrm {mm}$) from experiments (circle) and linear forced oscillation model (line, (4.6)).
This observed difference can also be attributed to other possible reasons. For example, (4.6) is a linear one-dimensional model with both 
$We$ and 
$\alpha ^{*}$ being scalars. In our experiments, 
$We_{vg}$ has an implicit direction that follows the largest compression direction of the rate-of-strain tensor, and 
$We_{slip}$ should be aligned with the slip velocity direction. Their combined effect to bubble deformation may not follow a simple relationship of 
$\max(We_{vg},We_{slip})$. In certain circumstances, they could potentially work against each other, which is not accounted for in the linear oscillator model.
4.5. Breakup probability
One condition that is implicitly assumed in many breakup models (Martínez-Bazán et al. Reference Martínez-Bazán, Montanes and Lasheras1999b) is that all bubbles will eventually break, just a matter of time, as long as 
$We>We_{crit}$. This implies a breakup probability (
${p_b}$) close to 100 %, which should be valid for large 
$We\gg 1$ and 
$D\gg D_H$. However, for Hinze-scale or sub-Hinze-scale bubbles with 
$We\lesssim 1$ and 
$D\lesssim 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 
$We$ is one step further from the Hinze's seminal work, in which the ensemble-averaged 
$\langle We\rangle$ was adopted to quantify the breakup probability. The limitation of using 
$\langle We\rangle$ is that, based on 
$\langle We\rangle$ being larger or smaller than the critical 
$We_{crit}$, 
${p_b}$ is close to a step function (
${p_b}=1$ if 
$\langle We\rangle >We_{crit}$; 
${p_b}=0$ if 
$\langle We\rangle <We_{crit}$). However, in turbulence, local flows could be orders of magnitude stronger than the mean; bubbles could break in response to the local 
$We$ 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 
$We$ distribution to bubble breakup probability.
Before estimating 
$p_b$, we would like to extend the PDF of 
$We_{vg}$ and 
$We_{slip}$ beyond our experiments to other turbulent flows with different 
$\langle \epsilon \rangle$ and bubble size 
$D$. Based on (4.2) and (4.4), the distribution of local 
$We$ for three 
$\langle \epsilon \rangle$ from 
$0.1\ \textrm {m}^{2}\ \textrm {s}^{-3}$ to 
$10\ \textrm {m}^{2}\ \textrm {s}^{-3}$ are shown in figure 13(a). Both 
$We_{slip}$ and 
$We_{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 
$We$ increases as 
$\langle \epsilon \rangle$ grows.
Figure 13. (a) The predicted distribution of both 
$We_{slip,x}$ and 
$We_{vg}$ for different energy dissipation rates from 
$\epsilon =0.1$ to 
$10 \ \textrm {m}^{2}\ \textrm {s}^{-3}$. (b) Breakup probability 
$p_b$ calculated based on different mean Weber number 
$\langle We\rangle$. Three sets of lines indicate three different 
$We_{crit}$ from 1 to 4. Within each set, 
$p_b$ based on the total Weber number (cyan), or either 
$We_{vg}$ (red) or 
$We_{slip}$ (blue) alone, are shown. The curves from the total 
$We$ were fitted with (4.11) to predict 
$p_b$ for any mean 
$We$ and any selected 
$We_{crit}$.
From the distribution, we can estimate 
$p_b$ based on
\begin{equation} p_b(\langle We\rangle)=\int^{+\infty}_{We_{crit}}p(We) \,\textrm{d}(We), \end{equation}
where the local 
$We$ could be either 
$We_{vg}$ or 
$We_{slip}$, and 
$\langle We\rangle =\int ^{+\infty }_{-\infty }[We\times p(We)] \,\textrm {d}(We)$. Here 
$p_b$ is written as a function of the mean Weber number 
$\langle We\rangle$ to help with other experiments or simulations that have access only to 
$\langle We\rangle$ from the mean energy dissipation rate. It can be seen that, as 
$\langle We\rangle$ grows, either due to a larger 
$D$ or larger 
$\epsilon$, the bubble breakup probability 
$p_b$ will increase.
Figure 13(b) shows 
$p_b$ as a function of 
$\langle We\rangle$. In most previous works assuming that the breakup probability has a sharp transition at 
$We_{crit}$, 
$p_b$ would behave like a step function: 
$p_b=1$ for 
$\langle We \rangle \geq We_{crit}$ and 
$p_b=0$ for 
$\langle We \rangle < We_{crit}$. Although these two limits still apply in figure 13(b), the transition is much smoother, spanning over a few orders of magnitude of 
$\langle We \rangle$. Note that the choice of 
$We_{crit}$ does not affect the shape of the curve. As shown in figure 13(b), when we change 
$We_{crit}$ from 1 to 4, it just shifts the transitional 
$\langle We\rangle$ towards the new 
$We_{crit}$ without affecting the overall trend.
This framework applies to 
$We_{vg}$ and 
$We_{slip}$, both of which contribute to bubble breakup. As 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 
$We$ alone approach 0.5 for 
$\langle We\rangle$ much larger than 
$We_{crit}$ to ensure that the sum of the two 
$p_b$ equals to one. Furthermore, even for the same 
$We_{crit}$, 
$p_b$ of 
$We_{slip}$ (blue lines) stays mostly above that of 
$We_{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 
$We_{vg}$, the PDF of 
$We_{slip}$ has a larger probability for small 
$We$. 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 
$We_{vg}$ and 
$We_{slip}$ and using 
$We_{crit}=1$ to 4 are shown as three cyan lines in figure 13. These three lines are fitted with the same switch function:
\begin{equation} p_b=[1 + (2.8\langle We\rangle/We_{crit})^{-1.7}]^{-1} \end{equation}
which includes 
$We_{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 range of 
$\langle We\rangle$ for three 
$We_{crit}$ considered. Equation (4.11) provides a method to estimate bubble breakup probability in turbulence, particularly for bubbles close to the Hinze scale and 
$\langle We \rangle \approx 1$.
5. 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 finite-sized bubbles and surrounding turbulence. In this study, both 3-D bubble geometry and nearby 3-D 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. As the 3-D 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 
$\tilde {\epsilon }$. The modeled distribution of both 
$\tilde {\epsilon }$ and 
$We_{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 
$\langle \epsilon \rangle$ 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 nonlinear 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 generalised 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.
Acknowledgements
We acknowledge the financial support from the National Science Foundation award numbers 1854475 and CAREER-1905103. We would also like to acknowledge C. Meneveau for suggestions.
Declaration of interests
The authors report no conflict of interest.
