%~Mouliné par MaN_auto v.0.37.1 (6604bb9b) 2025-11-28 13:20:22
\documentclass[CRMECA,Unicode,biblatex,published]{cedram}

\addbibresource{CRMECA_Riviere_20250747.bib}

\usepackage{siunitx}

\DeclareMathOperator{\erf}{erf}
\DeclareMathOperator{\We}{We}
\DeclareMathOperator{\Rey}{Re}
\DeclareMathOperator{\Cov}{Cov}
\DeclareMathOperator{\Var}{Var}

\newcommand\D{d}

\newcommand{\boldif}{\mathop{}\!\mathbf{d}}

\newcommand\R{R_0}

\newcommand\force{\mathcal{T}}

%%%---------------------------------------------------------------------------

%%% DELIMITERS

\usepackage{mathtools}
\DeclarePairedDelimiter{\parens}{\lparen}{\rparen}
\DeclarePairedDelimiter{\abs}{\lvert}{\rvert}
\DeclarePairedDelimiter{\bracks}{[}{]}

%%%---------------------------------------------------------------------------

\graphicspath{{./figures/}}

\newcommand*{\mk}{\mkern -1mu}
\newcommand*{\Mk}{\mkern -2mu}
\newcommand*{\mK}{\mkern 1mu}
\newcommand*{\MK}{\mkern 2mu}

\hypersetup{urlcolor=purple, linkcolor=blue, citecolor=red}

\let\oldforall\forall
\renewcommand*{\forall}{\mathrel{\oldforall}}


\title{Bubble breakup probability in turbulent flows}
\alttitle{Probabilité de cassure pour des bulles en turbulence}

\author{\firstname{Ali\'enor} \lastname{Rivi\`ere}\CDRorcid{0000-0002-5658-0759}\IsCorresp}
\address{LFMI, School of Engineering, EPFL, Lausanne, 1000, Switzerland}
\address{PMMH, CNRS, ESPCI Paris, Universit\'e PSL, Sorbonne Universit\'e, Universit\'e de Paris, 75005 Paris, France}
\email{alienor.riviere@epfl.ch}
\thanks{A.~Rivière is funded by the SNSF Swiss Postdoctoral Fellowship (project number~233920)}
\CDRGrant[SNSF]{233920}

\author{\firstname{St\'ephane} \lastname{Perrard}\CDRorcid{0000-0002-9315-2892}}
\address[2]{PMMH, CNRS, ESPCI Paris, Universit\'e PSL, Sorbonne Universit\'e, Universit\'e de Paris, 75005 Paris, France}
\email{stephane.perrard@espci.fr}
\thanks{S.\ Perrard aknowledges financial support from PSL University through a Junior Fellowship 2022-305 and from the National Research Agency (reference LASCATURB ANR-23-CE30-0043-03)}
\CDRGrant[National Research Agency]{ANR-23-CE30-0043-03}

\keywords{\kwd{Bubble} \kwd{turbulence} \kwd{breakup} \kwd{stochastic model}}
\altkeywords{\kwd{Bulle} \kwd{turbulence} \kwd{cassure} \kwd{modèle stochastique}}

\begin{abstract}
Bubbles drive gas and chemical transfers in various industrial and geophysical contexts, in which flows are typically turbulent. As gas and chemical transfers are bubble size dependent, a robust quantification requires a prediction of bubble breakup. The most common idea, introduced independently by Kolmogorov in 1949 and Hinze in 1955, is to consider a sharp limit between breaking and non breaking bubbles, given by $\operatorname{We}_c\approx 1$, where the Weber number $\operatorname{We}$ is the ratio between inertial and capillary forces at the bubble scale. In the statistically stationary state, $\operatorname{We}_c$ sets the maximum bubble size. Yet, due to the inherent stochasticity of the flow every bubble might in reality break. In this work, we use a stochastic linear model previously developed to infer the breakup probability of bubbles in turbulence as function of both $\operatorname{We}$ and the residence time. This allows us to introduce a definition of the critical Weber number accounting for the time spent by bubbles within a turbulent region, extending the stationary description. We show that bubble breakup is a memoryless process, whose breakup rate varies exponentially with $\operatorname{We}^{-1}$. The linear model successfully reproduces experimental breakup rates from the literature. We show that the stochastic nature of bubble breakup is central when the residence time of bubbles is smaller than ten correlation times of turbulence at the bubble scale: the transition between breaking and non breaking bubbles is smooth and most bubbles can break. For large residence times, the original vision of Kolmogorov and Hinze is recovered.
\end{abstract}

\begin{altabstract}
Les bulles contrôlent les réactions chimiques et les transferts de gaz pour un grand nombre de situations industrielles et naturelles caractérisées, typiquement, par des écoulements turbulents. Les transferts de gaz aux interfaces et les réactions chimiques étant fonction de la taille des bulles, pour les décrire précisément il est nécessaire de connaître leur distribution de taille et en particulier de prédire leur cassure. L'idée la plus répandue, introduite indépendamment par Kolmogorov en 1949, puis Hinze en 1951, consiste à considérer une limite stricte entre bulles stables et instables, donnée par $\operatorname{We}_c \approx 1$, où le nombre de Weber, $\operatorname{We}$, est le ratio entre forces inertielles et capillaires à l'échelle de la bulle. En régime statistiquement stationnaire, $\operatorname{We}_c$ fixe une taille maximum. Or, de part le caractère intrinsèquement fluctuant de la turbulence, toutes les bulles pourraient en réalité casser. Dans ce travail, nous utilisons un modèle linéaire stochastique développé précédemment pour prédire la probabilité de cassure en fonction du $\operatorname{We}$ mais aussi du temps passé par les bulles dans l'écoulement turbulent, nommé le temps de résidence. Ceci nous permet d'introduire un $\operatorname{We}_c$ qui prend maintenant en compte le temps de résidence, étendant ainsi la vision stationnaire. Nous montrons que le processus de cassure est un processus sans mémoire, associé à des taux de cassure variant exponentiellement en $\operatorname{We}^{-1}$. Ce modèle linéaire reproduit les taux de cassure expérimentaux présents dans la littérature. Nous montrons que la nature stochastique de la cassure est centrale lorsque les temps de résidence sont inférieurs à 10~temps de corrélation de la turbulence à l'échelle de la bulle\tralicstex{}{~}: la transition entre bulles qui cassent et qui ne cassent pas est continue et la majorité des bulles peuvent casser. En revanche, pour des temps longs, nous retrouvons la limite stricte introduite par Kolmogorov et Hinze.
\end{altabstract}

\dateposted{2025-12-02}
\begin{document}

\maketitle

\section{Introduction}

Bubbles play a crucial role in mass transport across interfaces. By increasing the surface of exchange they lead chemical and gas transfers in various industrial processes, as homogenizers~\cite{haakansson2019}, bubble column reactors~\cite{kantarci2005,han2007,risso2018}, and geophysical situations, such as rivers, waterfalls~\cite{beaulieu2012,demars2013}, and oceans~\cite{melville1996,woolf2007,deike2022}. Since gas exchanges are bubble size dependent, predicting breakups in turbulent flows is central to understand the bubble size distribution and then quantify the transfers in a robust manner. In inertial flows, bubble fate is controlled by the ratio between inertial and capillary forces, namely, the Weber number, $\We = \rho U^2 \D/\gamma$ where $\rho$ is the fluid density, $\D$ the bubble equivalent diameter, $U$ a typical velocity, and $\gamma$ the surface tension. Low Weber number bubbles are stable, while they break above a critical Weber number $\We_c$ which depends on the flow geometry~\cite{muller2008,chu2019,rodgar2023,riviere2023}. It is worth mentioning that in the case of dense turbulent emulsions~\cite{crialesi2023} defined the critical Weber number as the one for which the net energy flux from capillarity is zero. In turbulence, according to the theory developed independently by Kolmogorov~\cite{kolmogorov1949bubble} and Hinze~\cite{hinze1955}, bubbles are mainly deformed by eddies of their size so that $U$ can be estimated as the average velocity increment at the bubble scale $\bigl\langle \delta u(\D)^2\bigr\rangle$. In homogeneous and isotropic turbulence (HIT), for size lying within the inertial range, the velocity increment relates to the energy dissipation rate~$\epsilon$~\cite{kolmogorov1941}, which quantifies the energy transfers across scales, through $\bigl\langle \delta u(\D)^2\bigr\rangle = 2 (\epsilon \D)^{2/3} $. The Weber number then reads $\We = 2\rho \epsilon^{2/3}\D^{5/3}/\gamma$. Note that this definition relies on average quantities. It does not account for the inherent fluctuations of the flow. As a consequence, since any bubble might encounter a large enough pressure fluctuation which breaks it, in turbulence $\We_c$ is only defined in a statistical sense~\cite{vela2022,ni2023}. In addition, in practical situations, flows are inhomogeneous (bubble columns) or unsteady (breaking waves) or both, and the time spent by a bubble within a homogeneous turbulent region, called the residence time, affects the critical Weber number $\We_c$. To the best of our knowledge, the residence time has not yet been taken into account in theoretical models.

In this article, we aim to quantify the bubble breakup probability as a function of both $\We$ and the residence time in order to incorporate the effect of finite residence times in the definition of $\We_c$.
To do so, we quantify the statistics of bubble lifetime using a linear stochastic model for bubble deformations, whose parameters were inferred using direct numerical simulations (DNS)~\cite{riviere2024}. Our approach is based on the reminiscent idea that non linear effects may be negligible until the critical deformation for breakup is reached~\cite{risso1998}, as introduced by several authors for bubbles~\cite{risso1998,galinat2007,lalanne2019,masuk2021simultaneous} and drops~\cite{roa2023}. Conversely to our approach, all previous studies used models with parameters derived in the quiescent case, in particular for the damping rate~\cite{rayleigh1879,lamb1932}, which turns out to be significantly enhanced by the presence of turbulent fluctuations, as observed experimentally by Ravelet et~al.~\cite{ravelet2011}.

This article is organized as follows. We first introduce the stochastic linear model. Then, running our reduced dynamics, we measure the probability for a bubble to break as a function of We and the residence time. We deduce the associated breakup rate and compare it to experimental datasets. We eventually introduce a definition of the critical Weber number accounting for the residence time and discuss practical applications in inhomogeneous and unsteady flows.

\section{Linear stochastic oscillator}

For small amplitude deformations, the bubble surface can be described by the local radius $R(\theta, \phi)$ with $\theta$ the colatitude and $\phi$ the longitude. This parametrization holds as long as the local radius is mono-valued in $(\theta, \phi)$, corresponding to limited bubble deformations, and it will fail close to breakup. In a co-moving frame of reference, the local bubble radius $R$ decomposed onto the real spherical harmonics base $Y_\ell^ m(\theta, \phi)$ reads
\begin{equation}
R(\theta, \phi) = \R\bracks[\Bigg]{1 + \sum_{\ell=2}^\infty \sum_{m=-\ell}^\ell x_{\ell, m}(t)Y_\ell^m(\theta, \phi)}
\end{equation}
where $\R$ is the volume-equivalent bubble radius. Since the frame of reference moves with the bubble, the three harmonics $(\ell=1)$, encoding bubble translation, are null at all time. Note that the spherical harmonics decomposition is not invariant under frame rotation: rotating the reference frame redistributes the weight associated to a mode $(\ell, m)$ on all the modes $(\ell, m^\prime)_{m^\prime \in [-\ell, \ell]}$. In turbulence, bubble shape continuously reorients due to random fluctuations and torques. Since we do not prescribe any special frame orientation relative to the bubble shape, we cannot distinguish the dynamics of $x_{\ell, m}$ from the ones of $x_{\ell, m^\prime}$ with $m^\prime \neq m$. Hence the $2\ell +1 $ modes $\ell $ will have equivalent dynamics. Each mode amplitude $x_{\ell, m}$ evolves over time under the stochastic forcing imposed by the surrounding turbulent flow. In turbulence, bubbles are known to mainly deform and break following oblate-prolate oscillations, which correspond to the five modes $(\ell = 2)$~\cite{risso1998,perrard2021} illustrated on Figure~\ref{fig:evolution}a. In the following, we will therefore focus on the dynamics of these five modes. We denote by $x$ the amplitude of one of them. As an illustration, we show on Figure~\ref{fig:evolution}b the temporal evolution of one mode $\ell=2$ at $\We = 0.71$ measured in a DNS following the procedure described in~\cite{riviere2024}.

\begin{figure}[h!]
\includegraphics{fig1-shape-evolve.pdf}
\caption{(a)~Shape of deformation of the five modes $\ell =2$. (b)~Typical temporal evolution of a mode $\ell=2$, at $\We = 0.71$ measured in a DNS~\cite{riviere2024}. (c)~Typical temporal evolution given by the linear stochastic model defined by equation~\eqref{eq:model}, at the same Weber number $\We = 0.71$.}\label{fig:evolution}
\end{figure}

We recently showed~\cite{riviere2024} that each mode $\ell = 2$ dynamics can be described in a statistical sense by a linear oscillator randomly forced by the surrounding flow,
\begin{equation}\label{eq:model}
\ddot{x} + \Lambda \dot{x} + \Omega^2 x = \force(t),
\end{equation}
where lengths are normalized using the bubble equivalent radius $\R$, and times are normalized using the eddy turnover time at the bubble scale $t_c=\epsilon^{-1/3}d^{2/3}$, which encodes the correlation time of turbulent structures of bubble's size. The values of the damping coefficient $\Lambda$ and frequency $\Omega$ as well as the statistical properties of the effective turbulence forcing $\force$ were previously inferred from a DNS~\cite{riviere2024}. We summarize here the main results and assumptions. We found that the dimensionless damping coefficient $\Lambda=12$ (in units of $t_c$) is independent of the Weber number. At $\Rey(d)= \sqrt{2}\epsilon^{1/3}d^{4/3}/\nu = 124$, with $\nu$ the liquid kinematic viscosity, this value is about 7~times higher than the viscous dissipation in a quiescent fluid. Indeed, turbulence induces an effective turbulent damping, as observed for drop shape oscillations in the presence of an internal turbulent flow~\cite{berry2005,xiao2021}. As the bubble lies within the inertial range, where the flow is scale invariant, we assume that $\Lambda$, given in units of $t_c$, is also scale invariant. Note that, in physical units, this implies that the effective damping only varies with the velocity increments at the bubble scale. Experimentally, in a fluid of fixed viscosity, changing the bubble Reynolds number will change the eddy turnover time at the bubble scale and the damping parameter in physical units. The scale invariance hypothesis implies that $\Lambda$ is independent of the Taylor Reynolds number. Although there is no proof for bubbles that deformations are independent of $\Rey_\lambda$, there are some hints that it is the case for drops, as long as the drop $\Rey$ is large enough and its size lies within the inertial range. For drop breakup Vela \& Avila (2022)~\cite{vela2022} varied the Taylor Reynolds number by a factor~3 without seeing any effects on the drop breakup rate. The natural frequency remains unchanged compared to the quiescent value~\cite{rayleigh1879,lamb1932} so that $\Omega^2 = 192/\We $. The forcing term $\force$ models the erratic forcing from turbulence by a time-correlated random variable independent of the bubble properties, and therefore of $\We$. The statistically stationary forcing~$\force$ is fully characterized by its probability distribution function (pdf), and its auto-correlation function. We found that the pdf of $\force$ can be modeled by a hyperbolic secant distribution $1/(2\sigma_\force)/\cosh{\parens[\big]{\pi x/(2\sigma_\force)}}$ with a standard deviation $\sigma_\force = 20$ independent of $\We$. This pdf coincides with the pdf of the pressure mode $\ell =2$ integrated over a sphere of radius~$\R$ in the absence of bubble~\cite{riviere2024}. The important features of this pdf, visible on Figure~\ref{fig:comp}b, are the exponential tails: large amplitude fluctuations are frequent. We also found, by inverting the forcing spectrum, that the auto-correlation function of $\force$ decreases on a characteristic timescale given by the eddy-turnover time at the mode scale. More precisely, the auto-correlation function writes $\exp(-2\pi. 2^{2/3}t)(1 + 2\pi.2^{2/3}t)$, where the $2^{2/3}$ factor originates from the eddy turnover time at the scale $d/\ell$ with $\ell=2$. Note that the typical correlation time of the effective forcing is $0.2 t_c$, which is significantly shorter than the eddy turn-over time $t_c$. Eventually, equation~\eqref{eq:model} models adequately the dynamics of one mode $\ell=2$ of amplitude $x$.

\begin{figure}[h!]
\includegraphics{fig2-comp-model.pdf}
\caption{(a)~Standard deviation of modes $\ell=2$ as a function of the Weber number from the DNS~\cite{riviere2024} (black dots) and given by the model~\eqref{eq:model} (solid line). The dotted line shows $x_c$, the critical value used in Section~\ref{sec:breakup} to model breakup. The standard deviation crosses~$x_c$ at \mbox{$\We \approx 18$.} (b)~Normalized distributions of the mode amplitude obtained by running the model~\eqref{eq:model} for various Weber numbers. The black dotted line is the hyperbolic secant distribution, which is the probability distribution of the stochastic forcing $\force$. For comparison, the black solid line shows a Gaussian distribution.}\label{fig:comp}
\end{figure}

In this paper, we use the linear model to generate a large sample of mode dynamics. We first generate 5000 independent temporal sequences of $\force$ using the procedure described in Appendix~\ref{appendix:noise}. Each temporal sequence lasts 50$t_c$ and has a temporal resolution $\Delta t_1 = 0.01t_c$. This choice ensures to have 10 points per correlation time of the forcing. Starting with an initial condition $\parens[\big]{x(0),\dot x(0)}$, we integrate over small time steps $\Delta t_2/\Delta t_1 = 0.1$, the exact solution of equation~\eqref{eq:model} for a $\force$ piecewise constant. The dynamics become independent of the initial conditions on a timescale given by the maximum between the three following time scales: the bubble period $t_\gamma = \Omega/(2\pi)$, the eddy turnover time or correlation time of the forcing, and the dissipation time~$t_\nu$. The ratio between the bubble period and the eddy turnover time is a function of the Weber number: $t_\gamma /t_c \sim \We^2$. The ratio between the effective dissipation time and $t_c$ is a function of~$\Lambda$: $t_\nu/t_c \sim 1/\Lambda=0.08$. Hence the dynamics become independent of the initial conditions on a timescale given by $\max\parens[\big]{1, \sqrt{\We}, 0.08}$, in units of $t_c$. In the following, we focus on bubbles for which $\We\leq 200$. Hence, we use the algorithm to generate 5000 temporal evolutions of $x$ of length $50 t_c$, with initial conditions $x=0$, $\dot x = 0$, corresponding to an initially spherical bubble at rest, but only analyze times such that $t\geq 10 t_c$. The new initial conditions are then the mode amplitude and velocity at $t_i=10t_c$ which are both randomly distributed. Doing so, we directly access statistically stationary dynamics. A typical dynamics of the reduced model~\eqref{eq:model} is illustrated on Figure~\ref{fig:evolution}c for $\We = 0.71$. The evolution of the standard deviation $\sigma_x$ with the Weber number is captured by the linear model, as shown in Figure~\ref{fig:comp}a, which compares the DNS results from~\cite{riviere2024} with the linear model. The pdf of $x$ is shown on Figure~\ref{fig:comp}b for different values of the Weber number. For small $\We$, the pdf is close to the forcing pdf, the hyperbolic secant distribution. However, when $\We$ increases, the pdf becomes more and more Gaussian. Indeed, even if the dynamical equation is linear, it effectively acts as a filter in Fourier space of the forcing signal. The pdf of $x$ results from a convolution between the forcing statistics and the response function of the system.

\section{A memoryless process}\label{sec:breakup}

In order to model breakup, we introduce a critical linear deformation at which a bubble breaks. We choose $x_c = 0.74$ for consistency with our previous work~\cite{perrard2021}. Accordingly, we define the breakup time of one oscillator $x$ as the first time at which the amplitude reaches $x_c$. The cumulative probability, $p_b(t)$, is the probability that one mode reaches $x_c$ before the observation duration $t$.
As shown in Figure~\ref{fig:comp}a, for $\We\leq 18$, $\sigma_x \leq x_c$ and stochasticity will matter while for $\We \geq 18$, $\sigma_x \geq x_c$ and the dynamics will be mostly controlled by initial conditions. In particular, as shown in~\cite{perrard2021}, at large We, breakup results from the advection of the two sides of the interface toward one another, and happens in one eddy turnover time. We therefore focus on $\We \leq 18$. Figure~\ref{fig:proba} shows the survival probability associated to one mode $\ell =2 $, $1-p_b(t)$, as a function of the observation duration, for increasing Weber numbers. The survival probability decays exponentially in time, as expected for a memoryless process, with a breakup rate $\kappa$ independent of time, associated to each mode of deformation. For a breakup time much larger than the correlation time of both the velocity fluctuations at the bubble scale and the bubble deformations, the dynamics can indeed be modeled as a succession of independent events. Eventually, the breakup probability associated to one mode $\ell = 2$ reads $p_b(t) = 1 - \exp\bracks[\big]{-\kappa(\We) t}$. Note that, thanks to our choice of initial time, $t_i\geq 10 t_c$, the initial conditions are randomly sampled, and the time distributions are exponentially distributed from $t=0$. There is no equilibrium time in contrast to what has been observed in DNS of drops~\cite{vela2022} for which initial conditions are always a spherical drop at rest.

We now connect the breakup associated to one mode $\ell=2$ with bubble breakup. We assume that a bubble breaks when one of the five modes $\ell=2$ reaches the critical deformation. These modes, in the linear limit, are uncoupled and follow the dynamics given by equation~\eqref{eq:model}. As a consequence, the probability that a bubble breaks before $t$ is given by $P_b(t) = 1-\parens[\big]{1-p_b(t)}^5$ and can be written as
\begin{equation}
P_b(t) = 1 - \exp\bracks[\big]{-5\kappa(\We) t},\label{eq:proba}
\end{equation}
which extends to all Weber numbers the previous experimental observation of memoryless bubble breakups at large Weber number by Ravelet et~al.~\cite{ravelet2011} also observed numerically for drop breakups~\cite{vela2022}. Hence the bubble breakup rate $\kappa_b$ is five times larger than that of one single oscillator: $\kappa_b = 5 \kappa$. A similar argument has been used by Brouzet et~al.~\cite{brouzet2021} to model fiber breakup in turbulence. In their case, since a fiber is described as a collection of $N$ elementary rigid elements, the number of breakup modes can be approximated by the number $N$ of elements.

\begin{figure}
\includegraphics{fig3-1mprobabreakup-ttc.pdf}
\caption{Survival probability associated to one mode, as a function of time for several $\We$ (solid color lines) in the model~\eqref{eq:model}. Black dotted lines are exponential fits for each Weber number. Inset plot: associated breakup rate $\kappa$ as a function of $\We$. The black dotted line is an exponential fit, see equation~\eqref{eq:kappawe}. The dashed vertical line shows the limit below which $\sigma_x \geq x_c$, and the initial conditions control breakup.}\label{fig:proba}
\end{figure}

The breakup rate associated with one mode $\ell=2$, $\kappa$, varies exponentially with $\We^{-1}$ (see inset Figure~\ref{fig:proba}). It follows that the bubble breakup rate $\kappa_b$ reads
\begin{equation}\label{eq:kappawe}
\kappa_b = 5 \kappa = \kappa^\infty \exp\bracks[\Bigg]{-\frac{\We^t}{ \We}}
\end{equation}
with $\kappa^\infty = 1.57$ and $\We^t =7.20 $, two numerical factors obtained by a least-square fit. $\kappa^\infty$ is the breakup rate in the limit of infinitely large $\We$. Its inverse, $1/\kappa^\infty=0.64 $, of the order of the eddy turnover time, represents the associated bubble lifetime in the absence of surface tension and is controlled by the typical advection time of the two bubble sides toward each other~\cite{perrard2021}. The transitional Weber number $\We^t$ sets the transition between fast breaking and unlikely slow breakups. This law, which suggests a mechanism of random activation process, was first proposed by Coulaloglou and Tavlarides (1977)~\cite{coulaloglou1977} in the context of emulsions, based on the idea of drop-eddy collisions. Note that when bubble size lies within the inertial range of the turbulent cascade, bubble dynamics are inertial and we expect equation~\eqref{eq:kappawe} to hold for a wide range of Reynolds number, as similarly shown recently for drops~\cite{vela2022}. Figure~\ref{fig:compdata}a compares $\kappa_b$ with the breakup rates measured in various experimental conditions~\cite{martinez1999a,solsvik2015,vejravzka2018}, $\kappa_b^{\exp}$ in~\unit{s^{-1}}, made dimensionless by the eddy turnover time $t_c$ at the bubble scale, together with the DNS data from Rivi\`ere et~al.~\cite{riviere2021}. The experimental data points were originally gathered together in Zhong et~al.~\cite{zhong2024}, and are reproduced here with courtesy of Pr.~Rui Ni. Note that experimental breakup rates are estimated using different expressions, which induces an additional scattering of the data~\cite{haakansson2020}. Even though the data sets were obtained in different experimental facilities with various turbulence intensities corresponding to different bubble Reynolds number, our linear model originally derived for small amplitude perturbations, and extrapolated to $\We \geq 18$, follows the experimental points. This result confirms the intuition of Risso \& Fabre~\cite{risso1998} that bubble lifetime distribution can be predicted using a linear model. Note that the present model only describes independent break-up events. The lifetime distribution of rapid correlated events such as filament breaking~\cite{riviere2022} do not fall under the scope of model~\eqref{eq:proba}.

\begin{figure}[h!]
\includegraphics{fig4-breakuprates-map.pdf}
\caption{(a)~Comparison between our predicted breakup rate and several experimental datasets: MB99 = Mart\'inez-Baz\'an et~al.~(1999) \cite{martinez1999a}; S\&V15 = Solsvik (2015)~\cite{solsvik2015}; V18 = Vejravzka et~al.~(2018) \cite{vejravzka2018}; R21 = Rivi\`ere et~al.~(2021) \cite{riviere2021}. (b)~Probability for a bubble to break depending on its We and the residence time $t_R$ normalized by the eddy turnover time. Black lines are isoprobabilities separated by $0.1$. The red line corresponds to a breakup probability of $1/2$ and defines $\We_c(t_R)$.}\label{fig:compdata}
\end{figure}

\section{Breakup probability and critical Weber number}

In real flows, homogeneous turbulent regions have a finite extent in both space and time. As a consequence, bubbles remain only a finite time, the residence time, $t_R$ in~\unit{s}, within the turbulent region with a probability to break before $t_R$ of
\begin{equation}\label{unifying}
P_b(t_R, \We) = 1 - \exp \bracks[\Bigg]{\frac{- \kappa^\infty t_R}{ \epsilon^{-1/3}d^{2/3}} \exp \parens[\Bigg]{-\frac{\We^t}{\We} } }.
\end{equation}
This equation predicts that any bubble whose size lies within the inertial range can break, provided the residence time is long enough. However, in practice, residence times are finite and the breakup of bubbles at sufficiently small $\We$ is unlikely to occur. Indeed, one could expect that the velocity fluctuations are in reality bounded. A feature which is not captured by the forcing pdf used in the present model. In addition, other physical mechanisms such as dissolution might dominate the dynamics of small bubbles. The curve $P_b(t_R, \We_c) = 0.5$ defines the critical Weber number which now depends on the residence time. Figure~\ref{fig:compdata}b shows the probability to break as a function of both $\We$ and the residence time in unit of the eddy turnover time. For short residence times, $P_b$ gently varies with $\We$. As a consequence, there is a wide range of $\We$ for which a bubble may break ($0.1<P_b<0.9$). The probabilistic nature of the breakup process smooths the transition around $P_b=0.5$. At long residence time compared to the eddy turnover time ($t_R>40 t_c$), the limit between breaking and non breaking bubbles sharpens. We conclude that the probabilistic nature of bubble breakup in turbulence is essential at short resident times while the deterministic vision of Kolmogorov and Hinze is recovered for large resident times ($t>O(10 t_c)$).

Even though other mechanisms might dominate over turbulent fluctuations in some practical situations, such as shear and vorticity~\cite{joseph1998,muller2008,chu2019,riviere2023}, it is tantalizing to compare the results from our linear stochastic model to unsteady and inhomogeneous flows. In what follows, we discuss three different situations --- plunging jets, oceanic bubbles and industrial bubble columns --- which correspond to three different regions of the diagram of Figure~\ref{fig:compdata}b. For example, for a typical plunging jet, with speed~\qty{5}{\meter\per\second}~\cite{kiger2012} at the interface, the bubble jet depth is about~\qty{50}{\milli\meter}, leading to a residence time of~\qty{0.01}{s}. The energy dissipation rate can be estimated using the velocity fluctuations $u^\prime$, which is typically~\qty{0.1}{\percent} of the entrance velocity~\cite{kiger2012}, and the integral lengthscale which is the nozzle radius $D \sim \qty{10}{mm}$, giving $\epsilon = (u^\prime)^3/D\sim \qty{10}{\square\meter\per\cubic\second}$. As a consequence, the residence time is of the order of two eddy turnover times. At this residence time merely all bubbles can break (see Figure~\ref{fig:compdata}b). We expect the coexistence of large bubbles with much smaller ones, leading to a broad bubble size distribution, which is indeed the case below plunging jets and waterfalls.

In the oceans, underneath breaking waves, turbulence is sustained for approximately $1/3$ of the wave period, which is of the order of one to few seconds, with a turbulent intensity of~\qty{1}{\square\meter\per\cubic\second}. For centimetric to millimetric bubbles the residence time ranges from $6$ to $30$ eddy turnover time, corresponding to a critical Weber number between 3 and 5, consistent with the experimental measurements~\cite{deane2002}. In industrial bubble columns, bubbles rise with a velocity of about~\qty{0.1}{\meter\per\second} on meter columns, inducing a weak turbulence of around~\qty{0.1}{\square\meter\per\cubic\second}~\cite{kantarci2005}. The residence time for millimetric bubbles is then of the order of 1000 eddy turnover time, where the breaking transition is extremely sharp. In this regime, all bubbles will have $\We<\We_c\approx 1$.

\section{Conclusion}

In this article, we use a stochastic linear model for bubble deformations to quantify the bubble breakup probability. We show that the breakup statistics can be adequately modeled by considering a critical deformation for breakup for any of the five modes of deformation. This result validates the idea of Risso and Fabre~\cite{risso1998} that bubble lifetime can be predicted from a linear deformation model. We provide an explicit expression to compute the breakup probability as a function of the bubble Weber number and the resident time. We introduce a definition of the critical Weber number which accounts for the resident time, unifying the probabilistic vision of bubble breakup with
the sharp transition between breaking and non breaking bubbles from Kolmogorov and Hinze which holds in statistically steady state. One limit of the current model is the behaviour at a vanishing Weber number. As we did not introduce a cut-off mechanism to prevent a bubble of arbitrary small size to break, our model does not exhibit statistically steady state: for an infinite resident time, all bubbles will eventually break. One solution would be to introduce a cut-off in the pdf of the forcing, which would imply a vanishing break-up rate at a finite Weber number. Provided that turbulent fluctuations dominate over other breaking mechanisms, the model can be applied to practical situations with finite resident time such as unsteady and inhomogeneous flows.

\section*{Acknowledgments}

We thank Rui Ni and Shijie Zhong for sharing the data used to generate Figure~\ref{fig:compdata}a. We also thank Luc Deike, François P\'etr\'elis, S\'ebastien Gom\'e, Louis-Pierre Chaintron and Emmanuel Villermaux for fruitful discussions. This work was performed using HPC resources from GENCI-IDRIS (Grant 2023-AD012B14107). This work was also granted access to the HPC resources of MesoPSL financed by the Region Île-de-France and the project Equip@Meso (reference ANR-10-EQPX-29-01) of the programme Investissements d'Avenir supervised by the Agence Nationale pour la Recherche.

\section*{Declaration of interests}

The authors do not work for, advise, own shares in, or receive funds from any organization that could benefit from this article, and have declared no affiliations other than their research organizations.

\appendix

\section{Generation of a sequence of random forcing using Gaussian copulas}\label{appendix:noise}

In this section, we explain how to generate a realisation of $\force(t)$, the effective forcing from turbulence on bubble shape, using Gaussian copula.

\subsection{Principle and notations}

\subsubsection{Constructing \texorpdfstring{$\force$}{T} from a vector of Gaussian variables}

Let $\force$ be the forcing evolution that we want to evaluate and $\{ \force_i\}_{i \in [1, N]}$ its discrete counterpart, defined at times $\{ t_i\}_{i \in [1, N]}$. We assume the statistics of $\force$ to be stationary, so that all $\force_i$ share the same distribution and the covariance matrix of $\{ \force_i\}_{i \in [1, N]}$ is such that $\Cov(\force_i, \force_j)$ only depends on $\abs{t_i - t_j}$. We introduce the following notations:
\begin{itemize}
\item $\Phi_{\sigma_\force}$: probability density function of $\force$;
\item $\sigma_\force$: standard deviation of $\force$;
\item $A(t)$: autocorrelation function of $\force$;
\item $Y_i = \force_i/\sigma_\force$;
\item $\Sigma$: covariance matrix of $\{Y_i\}_{i \in [1, N]}$.
\end{itemize}
By definition of $\{Y_i\}_{i \in [1, N]}$, we have the following properties:
\begin{itemize}
\item $Y_i$ has distribution $\Phi_1$, denoted $\Phi$ in what follows;
\item $\Sigma_{i, j} = A\parens[\big]{\abs{t_i - t_j}}$;
\item $\Sigma_{ii} = 1$, ie the standard deviation of $Y_i$ is 1 for all i.
\end{itemize}
Finally, we introduce $F_\Phi$, the cumulative distribution of $\Phi$, and $F_\Phi^{-1}$ its inverse.

There is no method to directly compute $\{Y_i\}_{i \in [1, N]}$. Instead, we will construct $\{Y_i\}_{i \in [1, N]}$ using a vector of centered Gaussian variables $\{X_i\}_{i \in [1, N]}$: this is a Gaussian copula method. Since we aim at constructing a stationary process, all the $X_i$ also share the same normal distribution $\mathcal{N}$ with average 0 and standard deviation $\sigma^\prime$. We introduce these additional notations:
\begin{itemize}
\item $\Sigma^\prime$: covariance matrix of $\{X_i\}_{i \in [1, N]}$;
\item $F_\mathcal{N}$: cumulative distribution function of $\mathcal N$;
\item $F^{-1}_\mathcal{N}$: inverse cumulative distribution function of $\mathcal N$.
\end{itemize}
We can easily generate a Gaussian random vector $\{X_i\}_{i \in [1, N]}$ of any covariance matrix $\Sigma^\prime$. The idea is then to transform $\{X_i\}_{i \in [1, N]}$ into $\{Y_i\}_{i \in [1, N]}$. We know that we can construct $Z_i$ of distribution~$\Phi$ through
\[
Z_i = F^{-1}_\Phi \parens[\big]{F_{\mathcal{N}}(X_i)}.
\]
However, the covariant matrix of $\{Z_i\}_{i \in [1, N]}$ is not known a priori. As a consequence, the difficulty is to find the right $\Sigma^\prime$ so that after applying $F^{-1}_\Phi \parens[\big]{F_{\mathcal{N}}(\,\cdot\,)}$ to each $X_i$ we get a vector with covariance matrix $\Sigma$. In what follows we will determine $\Sigma^\prime$ and then deduce $Y_i$ using
\begin{equation}\label{eq:XtoY}
Y_i = F^{-1}_\Phi \parens[\big]{F_{\mathcal{N}}(X_i)}.
\end{equation}
The following expression will also be useful:
\begin{equation}\label{eq:YtoX}
X_i = F^{-1}_{\mathcal{N}} \parens[\big]{F_{\Phi}(Y_i)}.
\end{equation}

\subsubsection{Explicit expression of each function}

The distribution $\Phi$ is the hyperbolic secant distribution with standard deviation 1:
\[
\Phi (x) = \frac{1}{2}\frac{1}{\cosh\parens[\big]{\frac{\pi x}{2}}}.
\]
Its associated cumulative distribution function is
\[
F_\Phi = \frac{2}{\pi} \arctan\parens[\Bigg]{\exp\parens[\Bigg]{\frac{\pi x}{2}}},
\]
whose inverse is
\[
F_\Phi^{-1} = \frac{2}{\pi} \log\parens[\Bigg]{\tan\parens[\Bigg]{\frac{\pi x}{2}}}.
\]
The centered normal distribution $\mathcal{N}$ is
\[
\mathcal{N}(x) = \frac{1}{\sigma^\prime \sqrt{2\pi}} \exp \bracks[\Bigg]{-\frac{1}{2}\parens[\Bigg]{\frac{x}{\sigma^\prime}}^2}.
\]
Its cumulative distribution function is
\[
F_\mathcal{N}(x) = \frac{1}{2}\bracks[\Bigg]{1 + \erf \parens[\Bigg]{-\frac{x}{\sqrt{2}\sigma^\prime}}}
\]
whose inverse writes
\[
F^{-1}_\mathcal{N}(x) = \sigma^\prime \sqrt{2}\erf^{-1}[2x -1].
\]
We also introduce the joint law of two centered Gaussian variables $X_i$, $X_j$ such that $\Var(X_i) = \Var(X_j) = \Sigma^\prime_{ii} = (\sigma^\prime)^2$ and $\Cov(X_i, X_j) =\Cov(X_j, X_i) =\Sigma^\prime_{ij}$:
\begin{equation}\label{eq:jointlaw}
f_{ij}(x, y) = \frac{1}{2\pi \sqrt{\Delta_{ij}}}\exp\parens*{- \frac{\Sigma^\prime_{ii}(x^2 + y^2) - 2 \Sigma^\prime_{ij}xy}{2\sqrt{\Delta_{ij}}}}
\end{equation}
where $\Delta_{ij} = (\Sigma^\prime_{ii})^2 - (\Sigma^\prime_{ij})^2$.

Finally, the autocorrelation function of the forcing is,
\[
A(t) = \exp(-2\pi \cdot 2^{2/3}t)(1 + 2\pi \cdot 2^{2/3} t).
\]

\subsection{Determination of the covariance matrix \texorpdfstring{$\Sigma^\prime$}{Sigma'}}

\subsubsection{Diagonal terms}

We can explicitly compute $\Sigma_{ii}^\prime$ from $\Sigma_{ii}$ using equation~\eqref{eq:YtoX}:
\[
\begin{split}
\Sigma^\prime_{ii}
	& = \Var(X_i)
\\	& = (\sigma^\prime)^2
\\	& = \Var\parens[\big]{F^{-1}_{\mathcal{N}} \parens[\big]{F_{\Phi}(Y_i)}}
\\	& = \int_\mathrm{R}F^{-1}_{\mathcal{N}} \parens[\big]{F_{\Phi}(y)}^2 \Phi(y) \boldif y - \bracks[\Bigg]{\int_\mathrm{R}F^{-1}_{\mathcal{N}} \parens[\big]{F_{\Phi}(y)} \Phi(y) \boldif y}^2
\\	& = I_1(\Sigma_{ii}) - \parens[\big]{I_2(\Sigma_{ii})}^2.
\end{split}
\]
We estimate numerically these two integrals and find $I_1 = 1$ and $I_2 = 0$ (indeed the integrand of $I_2$ is odd). It follows that
\begin{equation}\label{eq:diag}
\Sigma_{ii}^\prime = (\sigma^\prime)^2 = \Sigma_{ii} = 1,
\quad \forall i.
\end{equation}

\subsubsection{Off-diagonal terms}

This time, we want to evaluate $\Sigma_{ij}^\prime$ with $i \neq j$. We do not know the joint law of $(Y_i, Y_j)$, which is too costly to be evaluated, but we have an expression for the joint law of two Gaussian variables (equation~\eqref{eq:jointlaw}). We then use equation~\eqref{eq:XtoY} to express $\Sigma_{ij}$ as a function of $\Sigma_{ij}^\prime$:
\[
\begin{split}
\Sigma_{i, j}
	& = \Cov(Y_i, Y_j)
\\	& = \Cov \bracks[\big]{F^{-1}_\Phi \parens[\big]{F_{\mathcal{N}}(X_i)}, F^{-1}_\Phi \parens[\big]{F_{\mathcal{N}}(X_j)}}
\\	& = \iint_{\mathrm{R}^2} F^{-1}_\Phi \parens[\big]{F_{\mathcal{N}}(x_1)}F^{-1}_\Phi \parens[\big]{F_{\mathcal{N}}(x_2)}f_{ij}(x_1,x_2) \boldif x_1  \boldif x_2 - \bracks[\Bigg]{\int_\mathrm{R}F^{-1}_\Phi \parens[\big]{F_{\mathcal{N}}(x)}\mathcal{N}(x) \boldif x}^2
\\	& = I_a(\Sigma_{ij}^\prime) + I_b^2.
\end{split}
\]
Numerically, we find $I_b = 0$ (the integrand is odd).

We estimate $I_a$ numerically (see Figure~\ref{fig:cov}) and find that $I_a(\Sigma_{ij}^\prime) = \Sigma^\prime_{ij}$ is an excellent approximation. It follows that
\begin{equation}\label{eq:offdiag}
\Sigma_{ij} = \Sigma_{ij}^\prime, \quad \forall i \neq j.
\end{equation}

\begin{figure}[!h]
\includegraphics[width=196pt]{SM1-Covij.pdf}
\caption{Variation of $\Sigma_{ij}$ with $\Sigma_{ij}^\prime$. The blue line is $I_a$. The black dotted line follows $\Sigma_{ij} = \Sigma_{ij}^\prime$.}\label{fig:cov}
\end{figure}

\subsubsection{Conclusion}

Combining equations~\eqref{eq:diag} and~\eqref{eq:offdiag} we have
\begin{equation}
\Sigma = \Sigma^\prime
\end{equation}
making the evaluation of $\Sigma^\prime$ in our case extremely fast and cheap.

\subsection{Verification}

To summarize, we first generate a Gaussian random vector $\{X_i\}_{i \in [1, N]}$ with covariant matrix,
\begin{equation}
\Sigma_{ij}^\prime = A\parens[\big]{\abs{t_i - t_j}}.
\end{equation}
Then we use equation~\eqref{eq:XtoY} to get $\{Y_i\}_{i \in [1, N]}$. Finally, $\force_i$ is obtained using
\begin{equation}
\force_i = \sigma_\force Y_i.
\end{equation}
Figures~\ref{fig:autocorr} and~\ref{fig:distrib} compare the autocorrelation function and the distribution of $\force_i$, respectively, obtained with the above procedure, to their theoretical expressions. The agreement is excellent.

\begin{figure}[!h]
\includegraphics[width=184pt]{SM2-autocorr.pdf}
\caption{Autocorrelation of $\force_i$, ie $\Sigma_{ij}^\prime$ as a function of time. The black dotted line is the targeted function. The blue line corresponds to the generated forcing signal.}\label{fig:autocorr}
\end{figure}

\begin{figure}[!h]
\includegraphics[width=190pt]{SM3-distribution_force.pdf}
\caption{Comparison between the distribution of the generated noise (in blue) and the theoretical target distribution (in black). The agreement is excellent.}\label{fig:distrib}
\end{figure}

\printbibliography

\end{document}