# Bubble shape oscillations in a turbulent environment

Aliénor Rivière1  Kamel Abahri1<EMAIL_ADDRESS>Stéphane Perrard1
1PMMH, CNRS, ESPCI Paris, Université PSL, Sorbonne Université, Université de
Paris, 75005, Paris, France

###### Abstract

We investigate bubble deformations in an homogeneous and isotropic turbulent
flow by means of direct numerical simulations of a single bubble in
turbulence. We examine interface deformations by decomposing the local radius
into the spherical harmonics base. We show that the linear dynamics of each
mode, (for low Weber number), can be modeled by a forced stochastic linear
oscillator. We measure the coefficients of the model directly from the modes’
statistics. We find that the natural frequency corresponds to the Rayleigh
frequency, derived in a quiescent flow. However, dissipation increases by a
factor 15 compared to the quiescent case, at $\real_{\lambda}=55$. This
enhanced dissipation originates from a thick boundary layer surrounding the
bubble. We demonstrate that the effective forcing, originating from the
integration of pressure over the bubble surface, is independent on bubble
deformability. Therefore, the interface deformations are only one-way coupled
to the flow. Eventually, we investigate the pressure modes’ statistics in the
absence of bubbles and compare them to the effective forcing statistics. We
show that both fields share the same pdf, characterized by exponential tails,
and a characteristic timescale corresponding to the eddy turnover time at the
mode scale.

##  1 Introduction

###  1.1 Broad context

Bubbly turbulent flows are widely used in industrial processes to enhance mass
transfers and chemical reactions. Examples involve bubble column reactors
(Risso, 2018) and emulsifiers (Håkansson, 2019, 2021) for instance. In
geophysical contexts, bubbles are known to control aerosol productions at the
ocean-atmosphere interface, while playing a major role in gas transfers
(Deike, 2022). In both industrial and natural situations, the knowledge of the
bubble size distribution and its temporal evolution is necessary to predict
mass transfers across bubble interfaces. As a consequence, the study of bubble
breakup in turbulence has received considerable attention since the pioneer
works of Kolmogorov (1949) and Hinze (1955). They predicted that, for bubbles
of size lying within the inertial range of the turbulent cascade, bubble
dynamics and breakup are primarily controlled by the balance between inertial
and capillary forces. This ratio defines the Weber number $\textrm{We}(d)=\rho
U^{2}d/\gamma$, where $\rho$ is the liquid density, $U$ a characteristic
velocity, $d$, the bubble volume equivalent diameter and $\gamma$ the surface
tension between gas and liquid. In turbulence, assuming bubbles break due to
velocity fluctuations at their scale, Kolmogorov and Hinze postulated that the
characteristic velocity $U$ is the average velocity increment at the bubble
scale $\langle\delta u^{2}(d)\rangle^{1/2}$. When kinetic and capillary forces
balance we have $\textrm{We}(d_{h})\approx 1\approx\textrm{We}_{c}$, which
defines a critical size, the Kolmogorov-Hinze scale $d_{h}$ separating
statistically stable bubbles ($d<d_{h}$) from unstable bubbles ($d>d_{h}$).

However, the main physical mechanism leading to breakup remained to be
understood. Sevik & Park (1973) proposed a resonant mechanism, in which bubble
breaks due to series of excitation at its natural frequency, while other
authors argue that large fluctuations are necessary for a bubble to break (Lee
et al., 1987; Luo & Svendsen, 1996; Wang et al., 2003). To address this
question, several authors describe bubble deformation dynamics, either with
the help of a linear damped harmonic oscillator (Risso & Fabre, 1998; Ravelet
et al., 2011; Masuk et al., 2021b) on the bubble Rayleigh modes (Rayleigh,
1879), or via a tensorial equation for the main bubble axis of deformations
(Masuk et al., 2021a). The latter assumes that bubble shape is mostly
ellipsoidal while the former allows any bubble shape and describes each mode
dynamics.

More generally, bubble deformations in turbulence are one over many examples
of the interaction between a turbulent flow and a deformable object. From
plants oscillations in the wind (De Langre, 2008), to disks (Verhille, 2022)
and fibers deformations in water (Rosti et al., 2018; Brouzet et al., 2021),
many studies have aimed at finding a reduced dynamics for the amplitude of the
relevant spatial modes of deformations, in the form of a damped harmonic
oscillator, randomly forced by turbulence. A usual approach is to model the
coefficients of an ordinary differential equation, as well as the statistics
of a random forcing term that model the turbulent forcing. Here, we propose to
directly measure these coefficients and the forcing statistics from the
deformation statistics of bubbles in turbulence, by performing numerical
experiments.

We first review the bubble oscillation dynamics in quiescent flows and their
phenomenological extensions to turbulent flows.

###  1.2 Bubble dynamics in quiescent flows

Rayleigh (1879) investigated the oscillation dynamics of inviscid drops in
vacuum and bubbles in a quiescent inviscid flow. In the linear limit of
deformation, the local radius of a bubble (or a drop) can be decomposed into
axi-symmetric modes using the basis of Legendre functions, which are indexed
by the integer $\ell\in[2,\infty]$. Rayleigh showed that the amplitude of each
mode $\ell$ follows an harmonic oscillator equation, with a characteristic
natural frequency writing,

$\omega_{\ell}^{2}=8(\ell-1)(\ell+1)(\ell+2)\frac{\gamma}{\rho d^{3}}$ (1)

for bubbles, with $f_{\ell}=\omega_{\ell}/(2\pi)$ the characteristic
frequency. Latter on, Lamb (1932) extended this work to gas bubbles
oscillating in a liquid of low kinematic viscosity, $\nu$. He found that
bubbles’ modes oscillate at the Rayleigh frequency with a damping rate
$\lambda_{\ell}$, which reads

$\lambda_{\ell}=8(\ell+2)(2\ell+1)\frac{\nu}{d^{2}},$ (2)

for bubbles of negligible inertia and viscosity. In three dimensions, bubble
shape can be decomposed into the real spherical harmonics base,
$Y_{\ell}^{m}(\theta,\phi)$, indexed by $\ell\in[2,\infty]$ and $m$ an integer
ranging from $-\ell$ to $\ell$, where $\theta$ and $\phi$ are the co-latitude
and longitudinal angles in spherical coordinates. The axi-symmetric modes of
Rayleigh and Lamb correspond to $m=0$. We denote the dimensionless amplitude
of the modes in the spherical harmonics base by $x_{\ell,m}$. The dynamics
found by Lamb (1932) and Rayleigh (1879) applies to each spherical harmonics
mode amplitude $x_{\ell,m}$, so that they follow a damped harmonic oscillator
equation with natural frequency $\omega_{\ell}$ and damping rate
$\lambda_{\ell}$ independent of $m$,

$\ddot{x}_{\ell,m}+\lambda_{\ell}\dot{x}_{\ell,m}+\omega_{\ell}^{2}x=0.$ (3)

When time is made dimensionless using the natural frequency $\omega_{\ell}$,
this equation reads,

$x_{\ell,m}^{\prime\prime}+p(\ell)\textrm{Oh}\,{x}_{\ell,m}^{\prime}+x=0$ (4)

where $p(\ell)=2\sqrt{2}(\ell+2)(2\ell+1)/[(\ell-1)(\ell+1)(\ell+2)]^{1/2}$
and, ′ stands for derivative with respect to the dimensionless time
$\omega_{\ell}t$. The Ohnesorge number $\textrm{Oh}=\mu/\sqrt{\rho\gamma d}$
with $\mu=\nu\rho$ compares viscous to capillary effects, and controls the
quality factor
$Q_{\ell}=\omega_{\ell}/\lambda_{\ell}\sim\textrm{Oh}\,\ell^{-1/2}$ of the
Lamb oscillations.

To estimate the damping rate of small oscillations, Lamb (1932) computed the
velocity gradients of the irrotational inviscid velocity field. Doing so, he
underestimated the dissipation rate, as shown latter by Miller & Scriven
(1968), as most of the dissipation takes place within the bubble boundary
layer, even when viscosity is low. Another approach is given by the normal-
mode analysis (Chandrasekhar, 1959; Reid, 1960; Chandrasekhar, 2013), for the
spherical harmonics modes. This theory predicts an evolution of the bubble
natural frequency and damping rate with the Ohnesorge number. No explicit
formulation can be however derived: one needs to solve a characteristic
equation for each value of Oh. This approach correctly takes into account
viscous effects but only holds at long times, presumably when oscillations
have already been completely damped, and do not describe transient dynamics.
Miller & Scriven (1968) demonstrated that, in the limit of vanishing
viscosity, the normal-mode solution converges to the irrotational one in the
bubble case. For drops, the same demonstration has been done by Chandrasekhar
(1959) and Reid (1960).

Latter on, Prosperetti (1977, 1980) unified the two approaches by studying the
initial-value problem of a drop or a bubble oscillating is an initially
quiescent flow. He demonstrated that, regardless of the value of Oh, the
damped harmonic oscillator dynamics of Lamb (1932) holds at short times
compared to the viscous timescale, $t\ll R_{0}^{2}/\nu$, where $R_{0}=d/2$ is
the bubble equivalent radius. On the other hand, the normal mode description
of Chandrasekhar (1959) holds at long times, $t\gg R_{0}^{2}/\nu$. At
intermediate timescales, Prosperetti (1977, 1980) demonstrated that the
dynamics is more complex due to the existence of a memory term in the equation
of motion of the modes, which couples the dynamics with the past evolution.

###  1.3 Bubbles deformations in turbulence

For a bubble immersed in a turbulent flow, additional dimensionless parameters
may control the deformation. Let us consider a bubble of negligible inertia
and viscosity, equivalent diameter $d$, immersed in a fluid of density $\rho$,
dynamic viscosity $\mu$, with surface tension $\gamma$. When the surrounding
flow field is an homogeneous and isotropic turbulent flow, characterized by an
energy dissipation rate $\epsilon$, and an integral length scale $L_{int}$,
the Buckingham’s $\Pi$ theorem predicts that the dynamics is controlled by
three dimensionless numbers. Choosing a set of dimensionless numbers which
decouples viscous effects from capillary effects, we obtain that a generic
measure of shape deformation $\delta$ can be written as,

$\frac{\delta}{d}=F\left(\textrm{We}(d),\real(d),\frac{d}{L_{int}}\right).$
(5)

where $F$ is a dimensionless function. The Weber number
$\textrm{We}(d)=2\rho\epsilon^{2/3}d^{5/3}/\gamma$ compares inertial and
capillary forces at the bubble scale. The Reynolds number
$\real(d)=\sqrt{2}\rho\epsilon^{1/3}d^{4/3}/\mu$ balances inertial and viscous
forces at the bubble scale. Finally, the ratio $d/L_{int}$ is the scale
separation between the bubble scale and the integral length scale. Note that
using $\epsilon$ and $d$ we can define a characteristic velocity
$U=\sqrt{2}(\epsilon d)^{1/3}=\langle\delta u^{2}(d)\rangle^{1/2}$, the
velocity increment at the bubble scale in homogeneous and isotropic turbulence
(Pope, 2000). When the bubble size lies within the inertial range of the
turbulent cascade, the surrounding flow is scale invariant and we expect the
dynamics to be independent of $d/L_{int}$. The bubble dynamics will then be
primarily controlled by the Weber number. In the presence of gravity $g$, one
must also include the Bond number $\textrm{Bo}=\rho gd^{2}/\gamma$, comparing
gravity to capillary effects. For simplicity, we will not consider gravity in
this study. This assumption is valid for bubble diameter smaller than the
capillary length $\sqrt{\gamma/(\rho g)}\sim 2\textrm{mm}$. In practice,
looking at temporal evolution of bubble deformation, our model may describe
shape oscillations slightly above the capillary length.

In this work we focus on bubbles which do not break, corresponding to a bubble
size $d$ within the inertial range of the turbulent cascade and $d<d_{h}$. For
a typical turbulent flow with $\epsilon=1\textrm{ m}^{2}\textrm{s}^{-3}$, and
$\textrm{We}_{c}\approx 3$,
$d_{h}=(\textrm{We}_{c}\gamma/(2\rho\epsilon^{2/3}))^{3/5}\approx 8\textrm{
mm}$ and $\real(d_{h})\approx 2300$. Note that
$\real(d_{h})\sim\rho^{1/5}\gamma^{4/3}/(\epsilon^{1/5}\mu)$ decreases as
$\epsilon$ increases for a given pair of liquid-gas. It is worth mentioning
that, as a consequence, an increase of the global Reynolds number of the flow
induces more viscous effects at the Hinze scale.

In order to predict bubble breakup, Risso & Fabre (1998) introduced a forced
linear damped oscillator equation to describe the dynamics of sub-Hinze
bubbles. Observing that the average deformation increases linearly with We, up
to the threshold for bubble breakup, they postulated that a linear dynamics
would be valid to describe bubble deformations up to this point. They assumed
that the deformed radius $R(t)$ evolves following

$\ddot{R}+\lambda\dot{R}+\omega^{2}R=F_{ex}(t)$ (6)

where $\lambda$ is a damping rate, $\omega$ a natural frequency and
$F_{ex}(t)$ an instantaneous forcing from turbulence. Bubble deformations and
breakup are mainly controlled by the second spherical harmonics modes
$\ell=2$, which correspond to oblate-prolate oscillation (Risso & Fabre, 1998;
Ravelet et al., 2011; Perrard et al., 2021). As a consequence, as a first
guess, they used the Rayleigh natural frequency of mode 2,
$\omega=\omega_{2}$, equation (1), and the Lamb damping rate
$\lambda=\lambda_{2}$, equation (2), even though these values only hold in a
quiescent irrotational flow. Then, following the original idea from Kolmogorov
(1949) and Hinze (1955), they assumed that the turbulent forcing from
turbulence scales as the square of the instantaneous velocity increment at the
bubble scale $\delta u(d,t)^{2}$, leading to a forcing $F_{ex}(t)=Kd\delta
u(d,t)^{2}$ from dimensional analysis, where $K$ is a numerical constant of
order 1. Doing so, they assumed that the presence of the bubble does not
strongly affect the flow properties, so that the flow statistics correspond to
the single fluid case. Expressing length in units of $d$, and time in units of
$1/\omega_{2}$, equation (6) is now written as

$r^{\prime\prime}+20\sqrt{2/3}\,\textrm{Oh}\,r^{\prime}+r=\tilde{K}\textrm{We}(t)$
(7)

where $\tilde{K}$ is also a constant of order 1 and
$\textrm{We}(t)=2\rho\delta u(d,t)^{2}d/\gamma$ is the instantaneous bubble
Weber number. This model is essentially the same as equation (4), with an
additional random forcing term. This equation has been widely used for bubbles
(Ravelet et al., 2011; Lalanne et al., 2019; Masuk et al., 2021a, b) and drops
(Galinat et al., 2007; Maniero et al., 2012; Håkansson, 2021; Roa et al.,
2023) oscillations in turbulence with the adequate expressions of the damping
rate and natural frequency.

However, there is no guarantee that the bubble natural frequency and damping
rate remain unchanged compared to the quiescent case. They may a priori depend
on both $\real$ and We. Indeed, surrounding flows are known to modify the
natural frequency and the damping rate. For instance, for bubbles in a
uniaxial inviscid straining flow, Kang & Leal (1988) showed that the flow
couples modes $\ell=2$ and 4, inducing a reduction linear in We, of the mode
$\ell=2$ natural frequency at linear order. Rivière et al. (2023) investigated
numerically the deformation dynamics of bubbles in a uniaxial straining flow
at large but finite Reynolds number. Together with the linear We-dependency,
they reported an additional Re-dependency of the natural frequency of mode
$\ell=2$. In addition, in inertial flows, bubble deformations are primarily
driven by Eulerian pressure increments at the bubble scale (Qureshi et al.,
2007), which do not share the same statistics than velocity increments
squared.

###  1.4 Outline of the present work: Infer bubble deformations dynamics from
data

In this paper, following Risso & Fabre (1998), we assume a linear damped
oscillator equation with a stochastic forcing for the oscillations of each
mode of bubble deformation. However, we do not presume any values for the
coefficients of equation (6) and the form of the forcing. Instead, we directly
measure from the deformations dynamics, the effective natural frequency and
damping rate and compare them to the quiescent values. We then deduce the
statistical properties of the effective forcing. To identify the origin of the
effective forcing, we study the statistics of the pressure field evaluated on
a sphere of bubble radius $R_{0}$. Eventually, we investigate the flow
structure around bubbles and the local dissipation rate to discuss the origin
of bubble dynamics dissipation in turbulent flows.

##  2 Bubble deformations in HIT

###  2.1 Numerical set-up: DNS of a single bubble in HIT

We perform direct numerical simulations of an incompressible gas bubble
immersed in an homogeneous and isotropic turbulent flow of an incompressible
liquid, using the open-source software Basilisk (http://basilisk.fr) (Popinet,
2003, 2009; Abu-Al-Saud et al., 2018). Density and viscosity ratios are set to
$850$ and $25$, respectively, close to air-water ratios. The simulation goes
in two steps. We first create an homogeneous isotropic turbulent flow by
solving the one phase incompressible Navier-Stokes equations with an
additional forcing term proportional to the velocity (Rosales & Meneveau,
2005). After a transient regime, the flow reaches a statistically stationary
homogeneous and isotropic turbulent state. The turbulent fluctuations are
characterized by the Taylor Reynolds number $\real_{\lambda}$, defined at the
correlation length of velocity gradients, namely the Taylor micro-scale
$\lambda=\sqrt{15\nu/\epsilon}\,u_{rms}$ Pope (2000), where $u_{rms}$ is the
root mean square of the velocity. The Taylor Reynolds number of the flow is
$\real_{\lambda}=u_{rms}\lambda/\nu=55$. We then extract snapshots of the flow
and use them as flow initial conditions for bubble injection. Snapshots are
separated by at least $6\,t_{c}$ to make sure initial conditions are
independent. The spherical bubble is injected at the center of the simulation
box by changing locally the density and viscosity. The bubble size is chosen
so that it lies within the inertial range of the turbulent cascade where the
flow is scale invariant. The bubble diameter to box length ratio is $0.13$.
During this second stage, forcing is maintained to sustain turbulence, but
only in the liquid phase to guaranty that bubble deformations only come from
the fluid forcing. In both steps, we use adaptive meshgrid refinement in order
to save computational time while resolving all the physical length scales of
the problem. The minimum grid size corresponds to 34 points per bubble radius.
Details of the numerical set-up as well as a convergence study can be found in
Rivière et al. (2021).

In this study we keep the flow Reynolds number constant and we vary the bubble
Weber number by changing the value of the surface tension coefficient. The
bubble Reynolds number is $\real(d)=124$. We explore eight values of We
ranging from $\textrm{We}_{c}\approx 3$ to $0.1\textrm{We}_{c}$. For each
Weber number, we run between 3 and 5 simulations. Except when the bubble
breaks (at $\textrm{We}=2.9$), we run every simulation for at least
$20\,t_{c}$, so that the total time per ensemble is about $100\,t_{c}$. Table
1 summarizes the exact number of simulations and total computational time per
Weber number we perform.

We | 2.9 | 2 | 1.43 | 1 | 0.71 | 0.46 | 0.36 | 0.27
---|---|---|---|---|---|---|---|---
N | 5 | 5 | 3 | 3 | 5 | 3 | 3 | 5
$T_{tot}/t_{c}(d)$ | 62 | 126 | 94 | 94 | 156 | 94 | 84 | 94
Table 1: Number of simulations and total simulated time per values of the
Weber number.

###  2.2 Modes of deformations

To quantify bubble deformations, we decompose the local bubble radius $R$ into
the real spherical harmonics base $Y_{\ell}^{m}(\theta,\phi)$, where $\ell$
and $m$ are the principal and secondary numbers respectively, and $\theta$ and
$\phi$ the co-latitude and longitude,

$R(\theta,\phi,t)=R_{0}\big{[}1+\sum_{\ell=2}^{\infty}\sum_{m=-\ell}^{\ell}x_{\ell,m}(t)Y_{\ell}^{m}(\theta,\phi)\big{]},$
(8)

and we track the modes’ amplitude $x_{\ell,m}$ over time. Bubble shape is
described in the bubble frame of reference so that all harmonics $\ell=1$,
corresponding to bubble translation, are null by definition. Numerically, the
bubble center is determined at each time step recursively by moving the frame
origin to minimize the amplitudes of all modes $\ell=1$. The recursion stops
when the center displacement between two steps is less than
$2.5.10^{-6}R_{0}$. This condition ensures that $|x_{1,m}|<4.10^{-6}$ for all
values of m. Note that the spherical harmonics decomposition holds as long as
the local radius $R$ is mono-valued. The procedure to compute the spherical
harmonics is described in detail in Perrard et al. (2021).

Figure 1: Typical temporal evolution for the mode $(2,0)$ at two different
Weber numbers. Time is made dimensionless using the Rayleigh frequency
$f_{2}$. Modes exhibit random oscillations, with an amplitude increasing with
We.

Figure 1 shows two typical temporal evolution for the mode $(\ell,m)=(2,0)$,
at two different Weber numbers. Time is made dimensionless using the Rayleigh
frequency $f_{2}$. For both We, we observe random oscillations around zero and
the predominance of the bubble resonant frequency $f_{2}$. The amplitude of
the oscillations increases with We.

Since we do not prescribe any special orientation relative to the bubble
shape, all modes with the same principal number $\ell$ are statistically
equivalent. Indeed, one can verify that a rotation of a mode can be expressed
as a linear combination of all the other modes with the same principal number.
As a consequence, we omit $m$ in what follows. For instance $x_{\ell}(t)$
represents a typical temporal evolution of one of the modes $\ell$. In
addition, assuming that $x_{\ell,m}$ are independent, the ensemble averaging
operation $\langle\cdot\rangle$ is computed over different simulations and
over the $m$ values for a given $\ell$. Roa et al. (2023) used a reference
frame dynamically oriented with the bubble principal axis of deformations. In
practice, their reference frame maximizes the amplitude of mode $(2,0)$, such
that the differential elongation can be studied as the invariance by rotation
is broken.

##  3 Determination of the reduced dynamics

###  3.1 Model: a stochastic linear oscillator

Following Risso & Fabre (1998), we introduce a linear stochastic model to
describe each mode dynamics,

$\ddot{x}_{\ell}+\Lambda_{\ell}(\textrm{We})\dot{x}_{\ell}+\Omega_{\ell}(\textrm{We})^{2}x_{\ell}=\mathcal{T}_{\ell}(\textrm{We},t),$
(9)

where $\Lambda_{\ell}$ and $\Omega_{\ell}$ are the damping rate and natural
frequency respectively and $\mathcal{T}_{\ell}$ is a random variable which
models the turbulent forcing. In this section, we aim at measuring
$\Lambda_{\ell}$, $\Omega_{\ell}$ and the statistical properties of
$\mathcal{T}_{\ell}$ from the deformation dynamics. Both parameters
$\Lambda_{\ell}$ and $\Omega_{\ell}$, as well as the forcing statistics, may
depend on the Weber number We. Conversely to what other authors have done,
time is made dimensionless using the eddy turnover time at the bubble scale
$t_{c}(d)=\epsilon^{-1/3}d^{2/3}$ and, from now on, $\dot{\cdot}$ denotes
derivatives with respect to this dimensionless time. This choice avoid a
priori to have a forcing term depending on bubble properties such as surface
tension. It decorrelates the turbulent forcing (righ hand side), from the
bubble response (the left hand side). In these units, the Rayleigh frequency
and the Lamb damping rate write
$\omega_{\ell}^{2}=16(\ell-1)(\ell+1)(\ell+2)/\textrm{We}$ and
$\lambda_{\ell}=8\sqrt{2}(\ell+2)(2\ell+1)\real(d)^{-1}$ respectively. Note
that, in this study, we have not varied the eddy turnover time. When the
bubble size lies within the inertial range of the turbulent cascade its
dynamics is primarily controlled by inertial effects, and the parameters may
not depend explicitly on the bubble Reynolds number, as long as $\real(d)\gg
1$.

In order to measure the coefficients and the force statistics of equation (9),
we make the following assumptions:

(H1):

Modes dynamics are linear and uncoupled, which is valid for $x_{\ell}\ll 1$,
corresponding to $\textrm{We}\ll 1$.

(H2):

The bubble deformation is one way-coupled to the flow. This hypothesis is
discussed in section 3.5.

(H3):

The forcing $\mathcal{T}_{\ell}$ is statistically stationary.

(H4):

The damping rate and the natural frequency do not depend on time.

From hypothesis (H2) the effective $\mathcal{T}_{\ell}$ in unit of the eddy
turnover time, is independent of We. From hypothesis (H3), the effective
forcing is completely determined by its auto-correlation function (or
equivalently its spectrum), and its probability distribution function (pdf).

Under these hypothesis, in the next sections, we will show that

  1. 1.

The natural frequency is not modified by the presence of the flow:
$\Omega_{\ell}=\omega_{\ell}$.

  2. 2.

There is an effective viscosity, driven by turbulence, so that
$\Lambda_{\ell}=0.6(\ell+2)(2\ell+1)$ for $\real(d)=124$.

Combining (1), (2) and equation (9) we will deduce the statistical properties
of the forcing $\mathcal{T}_{\ell}$.

###  3.2 Frequency response of the oscillator - Amplitude of the Fourier
transform

To rationalize the qualitative observations of figure 1 and identify the
angular frequency $\Omega_{\ell}$, we investigate the frequency response of
the bubble. To do so, we compute the temporal Fourier transform,
$\hat{x}_{\ell}$ of $x_{\ell}$ for all $\ell$ and We,

$\hat{x}_{\ell}(\omega)=\int_{-\infty}^{\infty}x(t)e^{-i\omega t}\differential
t,$ (10)

where $\hat{x}_{\ell}$ is also a random variable. Similarly, we introduce
$\hat{\mathcal{T}}_{\ell}$, the Fourier transform of the effective forcing,

$\hat{\mathcal{T}}_{\ell}(\omega)=\int_{-\infty}^{\infty}\mathcal{T}(t)e^{-i\omega
t}\differential t$ (11)

Figures 2 shows the ensemble average $\langle|\hat{x}_{\ell}|\rangle$ as a
function of the frequency, $f=\omega/(2\pi)$, normalized by the corresponding
Rayleigh frequency, $f_{\ell}$.

For $f<f_{\ell}$, for every We, $\langle|\hat{x}_{\ell}|\rangle$ is
approximately constant. The low frequency dynamics is similar to that of a
white noise.

At $f=f_{\ell}$ (back dotted line), for $\textrm{We}\leq 0.46$ we observe a
peak that resembles the resonant response of an oscillator at its natural
frequency. This peak does not exist for larger values of We. Nevertheless, for
every $\ell$, we observe a transition at this very frequency.

For $f>f_{\ell}$, at all We, we report a sharp power-law decay, following at
least $(f/f_{\ell})^{-4}$.

Finally, for $f>3f_{\ell}$, the spectrum amplitude is above the noise level.
Note that this part also corresponds to the end of the inertial range.

Dimensional measurements of bubble deformation dynamics were performed by
Ravelet et al. (2011) in the context of bubbles rising in turbulence. They
measured the temporal spectrum of the horizontal bubble main axis, a proxy for
the amplitude of the second Rayleigh mode. The overall shape of their power
spectrum was similar : weak variation for $f<f_{2}$, no resonance at $f_{2}$
and an a strong decay for $f>f_{2}$. In the absence of gravity, Risso & Fabre
(1998) also reported a transition at $f_{2}$, with a rapid decay for $f>f_{2}$
of the projected area spectrum.

The cut-off frequency being $f_{\ell}$ for all considered cases, we deduce
that the bubble natural frequency in turbulence, $\Omega_{\ell}$ of equation
(9), is not modified by the presence of the surrounding turbulent flow and
that,

$\Omega_{\ell}=\omega_{\ell}=4\left[\frac{(\ell-1)(\ell+1)(\ell+2)}{\textrm{We}}\right]^{1/2}$
(12)

It is surprising that the bubble natural frequency remains unchanged. Indeed,
Prosperetti (1980) showed, for a bubble in an initially quiescent flow, that
viscous effects induces an additional memory term in bubble dynamics. This
memory term can be modeled by an effective natural frequency and damping term.
The surrounding flow field can also modify the natural frequency. In a
uniaxial straining flow for instance, Kang & Leal (1988) demonstrated that a
coupling between modes $\ell=2$ and $\ell=4$ decreases the mode 2 natural
frequency at linear order, with a corrective term linear in We. We hypothesize
that the stochastic nature of turbulence cancels, in average, these
contributions. In the following, we use the theoretical expression of
$\omega_{\ell}$, for the bubble natural frequency, $\Omega_{\ell}$.

Figure 2: Amplitude of the modes’ Fourier transform for all We as a function
of the frequency normalized by the corresponding Rayleigh frequency. The Weber
number value is color-coded.

###  3.3 Zero frequency limit and We-dependency of the forcing

In this section, we investigate the zero frequency limit, and discuss the
consequence for the We-dependency of the forcing. By computing the Fourier
transform of equation (9), combined with (12), we obtain an expression linking
$\hat{x}_{\ell}$ and $\hat{\mathcal{T}}_{\ell}$,

$|\hat{x}_{\ell}|(\textrm{We},\omega)=\frac{|\hat{\mathcal{T}}_{\ell}|(\textrm{We},\omega)}{\sqrt{(\omega^{2}-\omega_{\ell}(\textrm{We})^{2})^{2}+\Lambda_{\ell}(\textrm{We})^{2}\omega^{2}}}$
(13)

The spectral behavior of each $x_{\ell}$ is a combination of the forcing
spectrum $\hat{\mathcal{T}}_{\ell}$ and the bubble response. In the limit case
$\omega=0$, using the expression of the bubble natural frequency (12), we have

$|\hat{x}_{\ell}|(\textrm{We},0)=\frac{|\hat{\mathcal{T}}_{\ell}|(\textrm{We},0)}{\omega_{\ell}(\textrm{We})^{2}}=\frac{\textrm{We}}{16(\ell-1)(\ell+1)(\ell+2)}|\hat{\mathcal{T}}_{\ell}|(\textrm{We},0).$
(14)

We can use this expression to investigate the We-dependency and
$\ell$-dependency of $\hat{\mathcal{T}}_{\ell}$ at $\omega=0$. We extract
$\langle|\hat{x}_{\ell}|\rangle(\textrm{We},0)$ by averaging
$\langle|\hat{x}_{\ell}|\rangle(\textrm{We},\omega)$ over the range
$5.10^{-3}<f/f_{\ell}<10^{-1}$ where the spectrum is constant.

Figure 3a shows $\langle|\hat{x}_{\ell}|\rangle(\textrm{We},0)$ as a function
of We. Solid lines of slope 1 are superimposed, showing that
$\langle|\hat{x}_{\ell}|\rangle(\textrm{We},0)$ increases linearly with We at
all $\ell$. It follows from equation (14) that
$\langle|\hat{\mathcal{T}}_{\ell}|\rangle(\textrm{We},0)$ is independent of We
for all $\ell$. This result justifies that the effective forcing from
turbulence does not depend on bubble deformability at low frequency. The
modification of the flow induced by bubble oscillations does not impinge back
on bubble dynamics. A similar phenomenon has been observed for drops by Vela-
Martín & Avila (2021). They investigated the interfacial stress generated by
eddies depending on their distance to the interface. They concluded that
eddies further that $0.2d$ from the drop interface (outer eddies) generate
most of the stress. They reported that these contributions are, in addition,
independent of We, as these eddies are two far from the interface to be
affected by drop deformations. We can assume that a similar mechanism may hold
also for bubble dynamics so that $\hat{\mathcal{T}}_{\ell}$ does not depend on
We either. These results justify hypothesis (H2) at low frequency, and we
assume that (H2) holds for all frequencies. From now on, we therefore assume
that $\mathcal{T}_{\ell}$ does not depend on We. This hypothesis will be
further validated and tested in section 3.5. The zero frequency limit also
depends on the mode order $\ell$. Figure 3b shows the compensated spectrum
limit $\langle|\hat{x}_{\ell}|\rangle(\omega=0)/\textrm{We}$ as a function of
$\ell$. We find that the zero frequency limit decreases slightly faster than
$\omega_{\ell}^{-2}\sim[(\ell-1)(\ell+1)(\ell+2)]^{-1}$ of equation (14) (red
line). It suggests that $|\hat{\mathcal{T}}_{\ell}|$ weakly decreases with
$\ell$, with $|\hat{\mathcal{T}}_{\ell}|\sim 1/\sqrt{\ell}$. Higher order
modes are associated to smaller scales that are less energetic. However, the
direct investigation of pressure modes in section 5.2 showed a much faster
decreases of the mode energy with $\ell$. The high order modes $\ell\geq 3$
may also be indirectly forced from non linear coupling with the mode 2
changing the $\ell$-dependency of the forcing.

Figure 3: a) Zero frequency limit of the modes’ Fourier transform as a
function of We for all $\ell$. Theoretical prediction
$\langle|\hat{x}_{\ell}|\rangle(\omega=0)\propto\textrm{We}$ is superimposed
in solid lines. Error bars are estimated using the standard deviation of the
spectrum value for $5.10^{-3}<f/f_{\ell}<10^{-1}$. b) Compensated limit
$\langle|\hat{x}_{\ell}|\rangle(\omega=0)/\textrm{We}$, as a function of
$\ell$. Colors encode the We (see figure 2). Assuming $\mathcal{T}_{\ell}$
independent of $\ell$ gives the scaling plotted in red. Assuming
$|\mathcal{T}_{\ell}|\sim 1/\sqrt{\ell}$ gives the scaling plotted in black.

###  3.4 Determination of the effective damping factor: Additional dissipation
due to turbulence

In this section, we present a method to compute the damping factor
$\Lambda_{\ell}$ of equation (9) from the numerical data.

Let $\hat{x}_{a}$ and $\hat{x}_{b}$ be the Fourier transform $\hat{x}_{\ell}$
of the same mode $\ell$ for two Weber numbers $\textrm{We}_{a}$ and
$\textrm{We}_{b}$. For simplicity here, we denote $\omega_{a}$ and
$\lambda_{a}$, the natural frequency and damping rate associated to
$\textrm{We}_{a}$ at this $\ell$. Under hypothesis (H2), the ratio $R_{ab}$

$R_{ab}(\omega)=\left(\frac{\langle|\hat{x}_{a}|\rangle}{\langle|\hat{x}_{b}|\rangle}\right)^{2}=\frac{(\omega^{2}-\omega_{b}^{2})^{2}+\Lambda_{b}^{2}\omega^{2}}{(\omega^{2}-\omega_{a}^{2})^{2}+\Lambda_{a}^{2}\omega^{2}}$
(15)

is independent on $\hat{\mathcal{T}}_{\ell}$.

Since the two natural frequencies $\omega_{a}$ and $\omega_{b}$ are known
(equation (12)), one can estimate the two damping factors, $\Lambda_{a}$ and
$\Lambda_{b}$, using $R_{ab}(\omega_{a})$ and $R_{ab}(\omega_{b})$, the ratios
at the two natural frequencies

$\displaystyle R_{ab}(\omega_{a})$
$\displaystyle=\frac{(\omega_{a}^{2}-\omega_{b}^{2})^{2}+\Lambda_{b}^{2}\omega_{a}^{2}}{\Lambda_{a}^{2}\omega_{a}^{2}}$
(16) $\displaystyle R_{ab}(\omega_{b})$
$\displaystyle=\frac{\Lambda_{b}^{2}\omega_{b}^{2}}{(\omega_{a}^{2}-\omega_{b}^{2})^{2}+\Lambda_{a}^{2}\omega_{b}^{2}}$
(17)

by solving this two-equations system. Note that an optimization of
$\Lambda_{a}$ and $\Lambda_{b}$ on the whole range of frequencies was less
reliable. The signal over noise ratio is optimal near the resonance, and
decreases both at high and low frequencies. Indeed, high frequencies, which
are the more noisy, then dominate the optimization procedure.

Figure 4a illustrates the computation of $\Lambda_{\ell}$. The ratio
$R_{0.71,0.27}$ for $\ell=2$, $\textrm{We}_{a}=0.71$ and
$\textrm{We}_{b}=0.27$ is represented as a function of the frequency $f$ (grey
curve). The black and red vertical lines denote the position of the two
natural frequencies $\omega_{a}$ and $\omega_{b}$ respectively, at which we
measure $R_{0.71,0.27}$. Inverting system (16)-(17) gives an estimate of
$\Lambda_{0.71}$ and $\Lambda_{0.27}$. Using these computed values of
$\Lambda_{0.71}$ and $\Lambda_{0.27}$ we plot the theoretical expression of
equation (15) at all frequencies (black line). This expression captures the
main features of the ratio $R_{0.71,0.27}(\omega)$: the low frequency limit,
the position and amplitude of the peak.

We then follow this procedure for every pair
$(\textrm{We}_{a},\textrm{We}_{b})$ and obtain 14 estimations of
$\Lambda_{\ell}$ per Weber number per mode $\ell$. We did not find a
significant bias on the estimated value of $\Lambda_{\ell}(\textrm{We})$ as a
function of the Weber ratio $\textrm{We}_{a}/\textrm{We}_{b}$. We then average
over all values of $\textrm{We}_{b}$ values to estimate
$\Lambda_{\ell}(\textrm{We}_{a})$. The values of $\Lambda_{\ell}$ for $\ell=2$
and $\ell=3$ as a function of We, and their standard deviation are reported in
table 2.

We | 2.9 | 2 | 1.43 | 1 | 0.71 | 0.46 | 0.36 | 0.27
---|---|---|---|---|---|---|---|---
$\Lambda_{2}$ | 14.2 | 11.7 | 11.9 | 11.1 | 11.0 | 13.8 | 11.9 | 17.1
$\sigma_{\Lambda}^{2}$ | 1.1 | 1.6 | 1.8 | 2.2 | 2.9 | 3.6 | 3.1 | 5.9
$\Lambda_{3}$ | 20.1 | 16.7 | 17.2 | 16.2 | 17.0 | 30.5 | 25.4 | 29.7
$\sigma_{\Lambda}^{3}$ | 3.4 | 4.0 | 4.6 | 4.4 | 5.3 | 13.2 | 12.6 | 18.4
Table 2: Average damping parameter $\Lambda_{\ell}$ and corresponding standard
deviation, for every We.

For $\ell\geq 4$, equation (13) fails to describe the ratio $R_{ab}$. Figure
4b shows $\Lambda_{\ell}$ as a function of We for $\ell=2$ and $\ell=3$ with
errorbars encoding the standard deviation $\sigma_{\Lambda}^{\ell}$. We find
no clear variation of $\Lambda_{\ell}$ with We, especially for $\ell=2$. When
$\ell$ increases, the dissipation also increases, as smaller scales are more
efficient to dissipate energy. The increase of $\Lambda_{\ell}$ with $\ell$ is
compatible with the $\ell$-dependency in a quiescent environment from Lamb
(1932). From our observations we found the following expression for the
damping factor,

$\Lambda_{\ell}=0.6(\ell+2)(2\ell+1).$ (18)

In quiescent environments, dissipation is also independent on We,
$\lambda_{\ell}=8\sqrt{2}(\ell+2)(2\ell+1)\real(d)^{-1}$, as it originates
from molecular diffusion in the liquid. However, we find $\Lambda_{2}\approx
13\lambda_{2}$. The surrounding flow field induces an additional effective
damping. Experimentally, Ravelet et al. (2011) also observed an additional
dissipation for bubbles rising in turbulence but attributed it to the presence
of the wake. Yet, similar observations come from drop oscillations in space.
In the presence of a turbulent internal flow, drop oscillations are
significantly damped (Bojarevics & Pericleous, 2003; Berry et al., 2005). This
additional dissipation is interpreted in terms of a turbulent eddy viscosity
(Xiao et al., 2021). In addition, Vela-Martín & Avila (2021) showed that there
is a transfer of energy from the drop interface to eddies closer than $0.2\,d$
from the drop interface and inside the drop. They call them inner eddies.
These small eddies efficiently dissipate energy. This transfer of energy
suggests that the enhanced dissipation comes from an increase in the local
effective diffusivity.

It is advantageous to estimate the size of an equivalent mixing length
$L_{t}$. This characteristic length of momentum transport has first been
introduced by Prandtl (Boussinesq, 1877; Prandtl, 1949; Xiao et al., 2021) to
describe, in a turbulent flow, the logarithmic profile of velocity near a
wall. By dimensional considerations, one can estimate the effective turbulent
viscosity $\nu_{t}$, using $L_{t}$ and a typical velocity scale of velocity
fluctuations at that scale, $\langle\delta u(L_{t})^{2}\rangle^{1/2}$,

$\nu_{t}=\langle\delta
u(L_{t})^{2}\rangle^{1/2}L_{t}=\sqrt{2}\epsilon^{1/3}L_{t}^{4/3}.$ (19)

Expressing the effective damping rate in terms of this effective dissipation
gives,

$\Lambda_{\ell}=8(\ell+2)(2\ell+1)\frac{\nu_{t}d^{2/3}}{d^{2}\epsilon^{1/3}}=8\sqrt{2}(\ell+2)(2\ell+1)\bigg{[}\frac{L_{t}}{d}\bigg{]}^{4/3}.$
(20)

Injecting equation (18), gives an estimate for $L_{t}$,

$L_{t}=\frac{d}{10}=\frac{R_{0}}{5}.$ (21)

Being of the same order of magnitude as the bubble radius, we hypothesize that
the mixing length originates from a geometric effect, similar to the
separation between inner and outer eddies from Vela-Martín & Avila (2021). We
further investigate the origin of this dissipation in the last section, by
looking at the local velocity gradients close to the bubble interface.

Figure 4: a) Ratio between the Fourier spectrum at $\textrm{We}_{a}=0.71$ and
$\textrm{We}_{b}=0.27$ for the mode $\ell=2$. The red and black vertical lines
denotes the position of the Rayleigh frequency at these two We where we
evaluate $R_{ab}$. The black line is the prediction from equation (15). b)
Damping factor as a function of We for $\ell=2$ and $\ell=3$. The error bars
represent the standard deviations. The solid black line corresponds to
$\Lambda_{2}=12$.

###  3.5 Effective forcing statistics: Temporal correlations and distribution

Since the left hand side of equation (9) is now completely determined, we can
compute the right hand side, and interpret it as a forcing term from the
turbulent flow.

For a bubble of size lying within the inertial range, the bubble deformations
are known to be primarily driven by Eulerian pressure gradients integrated
over the bubble surface (Qureshi et al., 2007). In the limit of small bubble
deformations, the pressure field at the bubble surface can be approximated by
the pressure field on a sphere of radius $R_{0}$. Decomposed in the spherical
harmonics base, this pressure reads

$p(\theta,\phi)=P_{c}\big{[}\sum_{\ell=0}^{\infty}\sum_{m=-\ell}^{\ell}P_{\ell,m}(t)\,Y_{\ell}^{m}(\theta,\phi)\big{]},$
(22)

where $P_{c}=\rho\delta u(d)^{2}$. There is no direct experimental measurement
of these pressure coefficients. Nevertheless, for the modes $\ell=1$, one can
circumvent this difficulty by looking at the acceleration statistics of
finite-size particles in a turbulent flow. Indeed, in the inertial range,
Lagrangian acceleration of particles reflects the statistics of the force
driving them, namely, the pressure gradient averaged over the particle volume
(Qureshi et al., 2007; Calzavarini et al., 2009). This hydrodynamic force
turns out to be the surface averaged pressure, which drives translational
bubble motions. For the higher order modes ($\ell\geq 2$), which drive bubble
deformations, there is no equivalent measurements in turbulence. Therefore,
from time to time, we will compare our statistics of $\mathcal{T}_{\ell}$
($\ell\geq 2$) with statistical quantities closely related to $P_{1}$, namely
the Lagrangian acceleration statistics and the pressure increments. A direct
measure of the statistics of $P_{\ell}$ ($\ell\geq 2$) in the absence of
bubble is provided in section 5.

Practically, we compute $\mathcal{T}_{\ell}$ from the modes’ Fourier transform
$\hat{x}_{\ell}$ using the following relation

$\mathcal{T}_{\ell}(t)=\frac{1}{2\pi}\int_{-\infty}^{\infty}\hat{x}_{\ell}(\omega)(\Omega_{\ell}^{2}-\omega^{2}+i\Lambda_{\ell}\omega)e^{i\omega
t}\differential\omega,$ (23)

where we use the expressions of $\Lambda_{\ell}$ and $\Omega_{\ell}$ from Eqs.
(12) and (18).

As expected from rotational invariance, we find that the average forcing
$\langle\mathcal{T}_{\ell}\rangle$ vanishes for all We. The standard deviation
of $\mathcal{T}_{\ell}$, $\sigma_{\mathcal{T}}^{\ell}$ is shown in figure 5 as
a function of We for $\ell=2$ and 3 (color-coded).
$\sigma_{\mathcal{T}}^{\ell}$ is found to be almost independent of the Weber
number. We found that the effective forcing from the turbulent flow does not
depend on bubble deformability. Therefore, bubble deformations are only one-
way coupled to the flow.

In physical units, the force $\mathcal{T}_{\ell}$ then scales as
$\alpha(\ell)\epsilon^{2/3}d^{-1/3}$, where $\alpha_{\ell}$ is a function of
the mode order. The standard deviation slightly decreases with $\ell$,
compatible with $\alpha_{\ell}\sim{\ell}^{-1/2}$.

In the context of Lagrangian particle acceleration in turbulence, the standard
deviation of acceleration also decreases with particle size as $d^{-1/3}$.
This scaling can be predicted using a scale invariant pressure fluctuations
argument (Voth et al., 2002; Qureshi et al., 2007; Volk et al., 2011). In
addition, Lagrangian acceleration statistics do not depend explicitly on the
Reynolds number at the particle size $\real(d)$, as long as the particle lies
within the inertial range. Only a marginal effect of the flow Taylor Reynolds
number $\real_{\lambda}$ on the variance of the acceleration (Voth et al.,
2002) was found. As a consequence, we expect the effective forcing to be
independent of $\real_{\lambda}$, $\real(d)$ and the Weber number.

Figure 5: Standard deviation of $\mathcal{T}$ as a function of We for $\ell=2$
and $\ell=3$. No We-dependency is observed, while $\sigma_{\mathcal{T}}$
decreases slightly for larger $\ell$. Figure 6: Probability density functions
of $\mathcal{T}_{2}$ (a) and $\mathcal{T}_{3}$ (b) for all We (color-coded),
normalized by their standard deviations. The shape of the distribution is
independent of We. The dashed line represents the hyperbolic secant
distribution, while the solid line is the Gaussian distribution.

Beyond the first two moments of the effective forcing distribution, it is
interesting to look at the full distribution. Figures 6a and 6b show the
probability distribution of $\mathcal{T}_{2}$ and $\mathcal{T}_{3}$
respectively for all We, normalized by their standard deviation
$\sigma_{\mathcal{T}}^{\ell}$. We find that the shape of the distribution is
also independent of the Weber number. These distributions are characterized by
exponential tails, and are well described by the hyperbolic secant
distribution (black dashed line),

$\text{pdf}(\mathcal{T})=\frac{1}{2\sigma_{\mathcal{T}}^{\ell}}\sech(\frac{\pi}{2}\frac{\mathcal{T}}{\sigma_{\mathcal{T}}^{\ell}})$
(24)

which depends on a single parameter, the standard deviation
$\sigma_{\mathcal{T}}^{\ell}$. The probability that a large forcing occurs is
way larger than that of a Gaussian distribution (solid black line).

It is again tantalizing to compare this distribution to Lagrangian
acceleration statistics for both particles and bubbles (Voth et al., 2002;
Qureshi et al., 2007; Volk et al., 2008; Homann & Bec, 2010; Prakash et al.,
2012; Salibindla et al., 2021). For small, neutral tracers and particles of
Kolmogorov scale size, the acceleration distributions exhibit larger tails,
decreasing slower than exponential. However, for larger particles
($d/\eta>10$), the shape exhibits exponential tail, independent of bubble size
and therefore of $\real(d)$ (Voth et al., 2002; Qureshi et al., 2007; Volk et
al., 2011). The pdf shape of the Lagrangian acceleration is well described by
the following expression, initially proposed for tracer particles (Mordant et
al., 2004; Qureshi et al., 2007)

$\text{pdf}(x)=\frac{e^{3s^{2}/2}}{4\sqrt{3}}\left[1-\erf\big(\frac{\log(|x/\sqrt{3}|)+2s^{2}}{\sqrt{2}x}\big{missing})\right]$
(25)

where $x$ is the standardized variable and $s$ an additional fitting
parameter. In the range of resolved scale, the two expressions, equations (24)
and (25), are compatible with our experimental data.

To characterize the temporal evolution of the effective forcing
$\mathcal{T}_{\ell}$, we study its ensemble averaged Fourier transform
$\langle|\hat{\mathcal{T}}_{\ell}|\rangle$. Injecting equations (12) and (18)
within equation (13) we obtain an expression in Fourier space for
$\langle|\hat{\mathcal{T}}_{\ell}|\rangle$:

$\langle|\hat{\mathcal{T}}_{\ell}|\rangle=\langle|\hat{x}_{\ell}|\rangle.\left[{(\omega^{2}-\Omega_{\ell}(\textrm{We})^{2})^{2}+\Lambda_{\ell}^{2}\omega^{2}}\right]^{1/2}.$
(26)

Figure 7a and 7b show $\langle|\hat{\mathcal{T}}_{2}|\rangle$ and
$\langle|\hat{\mathcal{T}}_{3}|\rangle$ respectively as a function of
$f\ell^{-2/3}$, where $\ell^{2/3}$ is the eddy turnover time at scale $d/\ell$
(in units of $t_{c}(d)$). For all frequencies, we found that the effective
forcing spectrum does not depend on the Weber number. At low frequencies
($f<0.2\,\ell^{2/3}$), the forcing amplitude is constant, corresponding to a
white noise. For $f>\ell^{2/3}$, the decay of
$\langle|\hat{\mathcal{T}}_{\ell}|\rangle$ is compatible with $1/f^{2}$. The
limit between these two regimes is set by the eddy turnover time at scale
$d/\ell$. We found that the spectrum of the effective forcing only depends on
the turbulence parameters, and is therefore independent of the bubble
deformations. As was anticipated in section 3.1, model (9) decouples the
turbulent forcing (the right hand side) from the bubble response (the left
hand side). The observation of a cut off frequency at the characteristic time
scale of turbulent fluctuations at the mode scale $d/\ell$ can be interpreted
as a filtering process originating from the integration over the bubble
surface. This filtering operation is further discussed in section 5.

From the previous observations, we propose the following expression for the
forcing spectrum,

$\langle|\hat{\mathcal{T}}_{\ell}|\rangle(f)=\frac{\tau_{\ell}}{1+[f\ell^{-2/3}]^{2}},$
(27)

where $\tau_{\ell}$ is a numerical constant, accounting for the
$\ell$-dependency of $\mathcal{T}_{\ell}$, that is adjusted on the data. From
equation 14 and figure 3, we estimate $\tau_{\ell}\sim\ell^{-1/2}$. The
expression (27) captures quantitatively the effective forcing spectrum (solid
black line in figures 7a and 7b).

Figure 7: Effective forcing spectrum for $\ell=2$ (a) and $\ell=3$ (b) deduced
from equation (26) as a function of the frequency normalized by the eddy
turnover time at scale $d/\ell$. The Weber number is color-coded with the same
colorbar as in figure 2).

In the context of Lagrangian particle accelerations, Voth et al. (2002)
followed by Volk et al. (2008), computed the temporal autocorrelation of
inertial particle accelerations in turbulence. The temporal acceleration
statistics of a finite size particle is usually attributed to a filtering
effect of the small scale turbulent fluctuations at the particle scale
(Qureshi et al., 2007). As a consequence, the correlation time of acceleration
for neutrally buoyant particle is given by the eddy turnover time $t_{c}(d)$.
This result has been recently extended to buoyant particle that exhibits a
modified correlation time $t\sim t_{c}(d)\beta^{-1/2}$ (Fan et al., 2024),
where $\beta=3\rho/(2\rho+\rho_{p})$ is a function of the fluid density $\rho$
and the particle density $\rho_{p}$. For a bubble, we have $\beta=3$,
corresponding to a correlation time of order $t_{c}$. In our case, the
temporal auto-correlation function
$C_{\mathcal{T}_{\ell}}(t)=\langle\mathcal{T}_{\ell}(0)\mathcal{T}_{\ell}(t)\rangle/(\sigma_{\mathcal{T}}^{\ell})^{2}$
for the modes $\ell>1$ can be deduced from the spectrum
$\hat{\mathcal{T}}_{\ell}$ and is written as

$C_{\mathcal{T}_{\ell}}(t)=\exp(-2\pi\ell^{2/3}t)(1+2\pi\ell^{2/3}t).$ (28)

We found that the correlation time in physical units is given by
$t_{c}(d)\ell^{-2/3}/(2\pi)$, which also scales as $t_{c}(d)$, with an
additional dependency in the mode order $\ell$. The prefactor being smaller
than one, the mode oscillations decorrelate faster that the velocity
fluctuations at the bubble scale.

In summary, we found that all the statistics of $\mathcal{T}_{\ell}$ are
independent of We which confirms the initial intuition of Risso & Fabre (1998)
that bubble dynamics and turbulent forcing are decoupled. We found that the
bubble deformation by the flow field can be described by a one-way coupling
model: the flow field generated by bubble oscillations does not significantly
impinges back on bubble dynamics. In addition, experimental results from the
literature suggest that these statistics are likely to be independent on
$\real(d)$, as long as we consider bubbles larger than the Kolmogorov scale.

From the stationary hypothesis (H3), the forcing is completely characterized
by its distribution and temporal autocorrelation function. The combination of
an explicit form for the pdf (eq. (24)) and for the autocorrelation function
(eq. (28)) then provides a complete model of a synthetic stochastic effective
forcing for bubbles deformations in turbulence. Previous modelling approaches
have used two points velocity measurements to model an effective forcing term
(Risso & Fabre, 1998; Lalanne et al., 2019; Masuk et al., 2021b), following
the original idea from Kolmogorov (1949) and Hinze (1955). Here we found that
the statistics of the effective forcing differ significantly from two points
statistics, in particular due to the volumetric filtering effect at the
particle size.

##  4 Model validation

To describe the bubble deformation, we have inferred step by step an equation
including, damping, natural frequency and a statistical model for the
effective forcing term $\mathcal{T}_{\ell}$. To validate our model, we check
against the DNS data the predictions of our linear model.

###  4.1 Modes’ standard deviation and distributions

We first look at the modes’ standard deviation $\sigma_{x}^{\ell}$ and
statistics. Figure 8a shows the modes’s standard deviation as a function of
the Weber number for $\ell\in[2,5]$. We find that $\sigma_{x}^{\ell}$ can be
approximated by
$\sigma_{x}^{\ell}\approx\textrm{We}/[(\ell-1)(\ell+1)(\ell+2)]$, with a
constant of order one. The prediction from our linear model computed from an
integration over the Fourier space, is superimposed in solid line for $\ell=2$
and 3 using eq. (9), showing a quantitative agreement with the numerical data.

A scaling for $\sigma_{x}^{\ell}$ as a function of We and $\ell$ can be
deduced analytically in model cases. One natural case is to consider
$\mathcal{T}_{\ell}$ as a Gaussian white noise of autocorrelation function
$C(t)=D\delta(t)$, where $\delta$ is the Dirac function, and $D$ is
independent of the Weber number. In this case, from the analysis of stochastic
harmonic oscillators (Gitterman, 2005) the standard deviation reads

$\sigma_{x}^{\ell}\sim\left[\frac{D}{\Lambda_{\ell}\Omega_{\ell}^{2}}\right]^{1/2}.$
(29)

From the coefficients $\Lambda_{\ell}$ and $\Omega_{\ell}$ we extracted, this
model predicts $\sigma_{x}^{\ell}\propto\textrm{We}^{1/2}$, which does not
correspond to the observed scaling. A finite correlation time has be taken
into account. We then consider $\mathcal{T}_{\ell}$ as an exponentially
correlated Gaussian noise of autocorrelation function
$\langle\mathcal{T}_{\ell}(t)\mathcal{T}_{\ell}(t^{\prime})\rangle=(\sigma_{\mathcal{T}}^{\ell})^{2}\exp(-|t-t^{\prime}|/t_{\ell})$,
where $t_{\ell}=\ell^{-2/3}/(2\pi)$ is the correlation time of the effective
forcing deduced from equation 28, and $D$ is independent of We. In this case
the mode’s standard deviation reads (Gitterman, 2005),

$\sigma_{x}^{\ell}=\sigma_{\mathcal{T}}^{\ell}\left[\frac{t_{\ell}(1+\Lambda_{\ell}t_{\ell})}{\Omega_{\ell}^{2}\Lambda_{\ell}(1+\Lambda_{\ell}t_{\ell}+\Omega_{\ell}^{2}t_{\ell}^{2})}\right]^{1/2}.$
(30)

The scaling of $\sigma_{x}^{\ell}$ now becomes a function of the ratios
$\Lambda_{\ell}t_{\ell}$ and $\Omega_{\ell}t_{\ell}$. In practice, we have
$\Omega_{\ell}t_{\ell}\gg 1$ and
$\Omega_{\ell}t_{\ell}\gg\Lambda_{\ell}t_{\ell}$ for sufficiently small Weber.
Considering the limit $\Lambda_{\ell}t_{\ell}\gg 1$, equation (30) simplifies
as

$\sigma_{x}^{\ell}=\frac{\sigma_{\mathcal{T}}^{\ell}}{\Omega_{\ell}^{2}}=\frac{\sigma_{\mathcal{T}}^{\ell}}{(\ell-1)(\ell+1)(\ell+2)}\textrm{We}.$
(31)

We then recover the observed scaling for small Weber number. For larger Weber
number, the ratio $\Omega_{\ell}t_{\ell}$ decreases, and we expect a
transition to a shallower increase of $\sigma_{x}^{\ell}$ with We. This
transition should occur for larger Weber number as $\ell$ increases, an
interpretation compatible with the numerical data shown in figure 8a. The
observed scaling of $\sigma_{x}^{\ell}$ with Weber thus corresponds to a
saturation of the bubble deformations dominated rather by the long correlation
time of the forcing (frozen turbulence hypothesis applied to bubble
deformations (Ruth et al., 2019)) than an accumulation of random forcing
events on a time scale $1/\Lambda_{\ell}$. It is worth noticing that the
estimate of the correlation time of the forcing is therefore essential to
predict the amplitude of bubble deformations.

To further check the dependency in $\ell$, figure 8b shows the compensated
standard deviation $\sigma_{x}^{\ell}/\textrm{We}$. We recover that the
decrease of the modes’ amplitude with $\ell$ can be mainly attributed to the
increase of the natural frequency with $\ell$, with a small correction
originating from the weak dependency of $\mathcal{T}_{\ell}$ with $\ell$.
Eventually, we found a quantitative agreement between the standard deviation
$x_{\ell}$ and the predicted value from the linear model. The model captures
the evolution of $\sigma_{x}^{\ell}$ with both We and $\ell$.

The linear increase of $\sigma_{x}^{\ell}$ with We, up to the critical Weber
number ($\textrm{We}\approx 3$ in our case) has important consequences when
modelling bubble breakup. Risso & Fabre (1998) suggested that the threshold
for breakup is close to the value above which the deformations start to be non
linear. A linear model would then be sufficient to describe bubble
deformations up to the breakup threshold.

Figure 8: a) Modes standard deviation as a function of the Weber number for
all $\ell$. The two straights lines are the predictions from our linear model,
for modes $\ell=2$ and $\ell=3$. b) Modes’ standard deviation compensated by
We as a function of the mode principal number.

We then look at the entire statistics of the $x_{\ell}$. Figures 9 show the
probability density functions of the modes $\ell=2$ for all Weber numbers (9a)
normalized by their standard deviation. We find that the shape of the pdf does
not depend on the Weber number and corresponds to the hyperbolic secant
distribution (black dashed line), equivalent to the pdf of the effective force
$\mathcal{T}_{\ell}$. Both the forcing and the mode amplitude share the same
pdf that deviate from gaussianity (solid black line) with exponential tails.
As the distributions exhibit fat tails, the probability that bubbles
experience large deformations leading to breakup is large compared to a
Gaussian distribution (black dotted line). Moreover, for larger $\ell$, the
deviation from gaussian distribution increases, as shown in figure 9b for
$\textrm{We}=1$.

Figure 9: a) Normalized pdf of $x_{2}$ for all We. b) Normalized pdf at
$\textrm{We}=1$ for different $\ell$. a)&b) The black dashed line is the
hyperbolic secant distribution. The solid black line the Gaussian
distribution.

###  4.2 Deformation spectrum

Figure 10 compares the modes’ Fourier transforms with the model prediction
(13) combined with equations (18), (12) and (27) (dotted lines). For all Weber
numbers, the model accurately predicts the zero-limit frequency as well as the
amplitude of the spectrum at the bubble natural frequency $f_{2}$ and the
position and slope of the decay at larger frequencies. At the lowest Weber
number ($\textrm{We}=0.27$), the model overestimates the spectrum just below
the resonance. We remind here that for frequencies larger than $2.5f_{2}$ the
spectrum is dominated by numerical noise. For all the other We, in the absence
of resonance, the model captures the spectrum close to the bubble natural
frequency.

Figure 10: Comparison of the Fourier spectrum amplitude between the DNS and
the model (dotted line). The model spectrum is obtained by combining equations
(18), (12) and (27) in (13), for the mode $\ell=2$. The model captures the low
frequency limit, the position of the transition as well as the high frequency
decay for all We.

###  4.3 Consequences for bubble breakup

Thanks to the quantitative model we develop, we can revisit the breakup
scenario and the criterion for breakup. Two main breakup scenarios have been
proposed for bubbles in turbulence. Bubbles can break either when they
encounter a pressure fluctuation larger than some threshold values (Lee et
al., 1987; Luo & Svendsen, 1996; Wang et al., 2003; Masuk et al., 2021a) or
after series of small excitation at their natural frequency which induce a
resonance (Sevik & Park, 1973; Risso & Fabre, 1998). The ability for a bubble
to store energy on a mode $\ell$, is quantified by the quality factor
$Q_{\ell}=\Omega_{\ell}/\Lambda_{\ell}$. For large quality factor typically
$Q_{\ell}>10$, modes can store energy over several periods of oscillation,
while for lower quality factor, energy is dissipated in few bubble periods at
most. Our linear model provides a quantitative measure of $Q_{\ell}$.
Combining equations (12) and (18) we have an explicit expression for
$Q_{\ell}$ as a function of We and $\ell$,

$Q_{\ell}=4\sqrt{\frac{(\ell-1)(\ell+1)}{0.6(\ell+2)(2\ell+1)^{2}}}\textrm{We}^{-1/2}.$
(32)

In turbulence, bubbles mainly break after oblate-prolate deformations, meaning
deformations along their second modes $\ell=2$ (Risso & Fabre, 1998; Ravelet
et al., 2011; Perrard et al., 2021; Masuk et al., 2021b). For the typical
critical Weber numbers reported in the literature, $0.1<\textrm{We}_{c}<10$
(Sevik & Park, 1973; Risso & Fabre, 1998; Martínez-Bazán et al., 1999; Rivière
et al., 2021), our estimate of the quality factor for the mode $\ell=2$ ranges
from $0.3$ ($\textrm{We}_{c}=10$) to 3 ($\textrm{We}_{c}=0.1$). These quality
factors are too low to observe significant energy storage over several period
of oscillations. We conclude that large pressure fluctuations set the value of
the critical Weber number rather than resonant mechanism.

Note that a sequence of oscillations at the bubble natural frequency may be
observed for sufficiently large quality factor, typically $Q_{2}>10$,
corresponding to $\textrm{We}<8.10^{-3}$. Even though such a Weber number
corresponds to bubbles size much smaller than the Kolmogorov Hinze scale,
which will never break, it may be observed experimentally.

##  5 Link between model coefficients and surrounding turbulent fields

In this section, we aim at connecting the effective variables we identified,
namely the forcing $\mathcal{T}_{\ell}$ and the damping coefficient
$\Lambda_{\ell}$, to relevant flow statistics in turbulence. For the effective
force, bubble deformations are known to originate from pressure differences
along the interface (Qureshi et al., 2007). The presence of a bubble modifies
the flow statistics in its surrounding, through dynamical boundary conditions
at the interface and incompressibility. Nevertheless, for drops, it has been
shown that the outer eddies (further than $0.2d$ from the interface) generate
most of the normal stress (Vela-Martín & Avila, 2021). These outer eddies may
be less affected by the presence of the interface. We will therefore assume
that the pressure statistics in the absence of bubble are a reasonable proxy
to estimate $\mathcal{T}_{\ell}$. On the other hand, the dissipation is
expected to originate from the boundary layer near the interface (Vela-Martín
& Avila, 2021). To rationalise the origin of the additional damping from the
flow statistics, we will therefore study the local dissipation in the bubble
boundary layer.

###  5.1 Point statistics of the pressure field

As a reference case, let us first consider the Eulerian point statistics of
pressure in homogeneous and isotropic turbulence. To compare with the bubble
dynamics, we will still express length scales in units of $d$, timescales in
units of $t_{c}(d)$ and therefore velocity in term of velocity increments at
the bubble scale $\langle\delta u(d)^{2}\rangle^{1/2}$.

We run single phase direct numerical simulations and record the Eulerian
pressure fluctuations $p(x,t)$ at seven different fixed location well
separated in space. We run three simulations for a total of $245t_{c}(d)$.
Resolution is increased compared to the two-phase problem and would be
equivalent to 68 points per bubble radius and 3.6 points per Kolmogorov
length. In this section, ensemble averages are performed over the three
simulations and the seven locations.

Figure 11a illustrates two temporal evolution of pressure, normalized by the
characteristic pressure difference at the bubble scale, $P_{c}=\rho\delta
u(d)^{2}$. We found a pressure standard deviation $\sigma_{p}=0.67P_{c}$.
Pressure exhibits random oscillations of small amplitude around zero, together
with large negative drops. This asymmetry between positive and negative
pressure fluctuations is better observed on the pressure pdf plotted on figure
11b. We recover that negative values are exponentially distributed, while
positive pressure values follow a Gaussian distribution (dashed black line).
The existence of the large negative peaks leading to an asymmetric pdf is well
known and has been reported both in experiments (Abry et al., 1994; Pumir,
1994; Cadot et al., 1995) and direct numerical simulations of homogeneous
isotropic turbulence (Cao et al., 1999; Vedula & Yeung, 1999). It has been
shown that these large negative peaks correspond to vorticity filaments
(Douady et al., 1991; Fauve et al., 1993; Cadot et al., 1995) passing through
the measurement point. As the bubble moves in the fluid, it may experience
different pressure statistics and the Lagrangian pressure statistics could
also be relevant.

Lagrangian pressure statistics have also been investigated numerically.
Numerical studies involve measuring pressure statistics along the paths of
point particles (Bappy et al., 2019), as well as (sub-Kolmogorov) finite-size
bubbles (Bappy et al., 2020b, a) whose dynamics are modeled using a pure
advection or a Maxey-Riley equation (Maxey & Riley, 1983; Toschi &
Bodenschatz, 2009) respectively. They found that larger particles have a
higher probability to be within low pressure regions. Nevertheless, the
overall shape of the pressure pdf, with an exponential tail for negative
values and a Gaussian distribution of positive values, is conserved.

Figure 11: a) Typical temporal evolution of the pressure at two points in
space. We observe small amplitude oscillations together with rare intense
negative peaks. b) Local pressure distribution normalized by its standard
deviation $\sigma_{p}=0.67P_{c}$. The solid black line follows the hyperbolic
secant distribution centered while the black dashed line follows a Gaussian
distribution with standard deviation $4/5$.

To investigate the frequency statistics of the local pressure, we compute the
temporal Fourier transform of each pressure signal $\hat{p}$,

$\hat{p}(\omega)=\int_{-\infty}^{\infty}p(t)e^{-i\omega t}\differential t.$
(33)

The average amplitude of its Fourier transform $\langle|\hat{p}|\rangle$ is
plotted on figure 12. The corresponding inertial range in frequency space is
delimited by the inverse of the eddy turnover time at the integral scale
$f_{c}(L)=1/t_{c}(L)$ (black dotted line) and the inverse of the eddy turnover
time at the Kolmogorov scale, $f_{c}(\eta)$ (dashed line). For low
frequencies, $f<f_{c}(L)$, $\langle|\hat{p}|\rangle$ slowly decreases with
$f$. Abry et al. (1994) have shown that this evolution at low frequencies
originates from the contribution of vorticity filaments, since their typical
lifetime is the integral timescale (Douady et al., 1991; Pumir, 1994).
Removing their contributions flattens the low frequency spectrum (Abry et al.,
1994).

In the inertial range of the turbulent cascade, $f_{c}(L)<f<f_{c}(\eta)$,
$\langle|\hat{p}|\rangle$ decays down to the noise level near $f_{c}(\eta)$.
In the spatial Fourier space, and a fortiori in the temporal Fourier space,
there is no consensus for the scaling of the pressure power spectrum within
the inertial range (Pullin & Rogallo, 1994). A Kolmogorov-like scaling
predicts $|\hat{p}(k)|^{2}\sim\epsilon^{4/3}k^{-7/3}$ (reported by Ishihara et
al. (2003) for instance) but other authors have also reported a $k^{-5/3}$
scaling (Vedula & Yeung, 1999; Gotoh & Rogallo, 1999). To transform the
spatial power spectrum into a temporal power spectrum, a classical way is to
consider that the small scale structures are advected by the large scales.
This is the sweeping hypothesis (Kraichnan, 1964; Tennekes, 1975), which has
been successfully used to reproduce pressure temporal autocorrelation (Yao et
al., 2008). Combining this argument with the Kolmogorov prediction, we find
that $\langle\hat{p}\rangle$ should scale as
$\hat{p}^{K}\sim\epsilon^{2/3}u_{rms}^{5/6}\omega^{-4/3}$, with a
proportionality constant of order 1. We find a reasonable agreement, as shown
by the compensated spectrum $\langle\hat{p}\rangle/\hat{p}^{K}$ in the inset
of figure 12. As evidenced by Pumir (1994); Pullin & Rogallo (1994); Vedula &
Yeung (1999) the Kolmogorov scaling might only hold in a narrow range of
frequencies, corresponding to scales just below the integral scale, due to the
limited inertial range. The proportionality constant is around 3 in our case
(solid black line) lower than the value of 7 proposed by Pumir (1994). The
third regime $f>f_{c}(\eta)$, corresponds to the end of the inertial range and
is close to the limit of resolution of our DNS, as $f_{c}({\Delta
x})=3f_{c}(\eta)$, where $\Delta x$ is the minimum grid size.

Figure 12: a) Amplitude of the local pressure Fourier transform. The vertical
dotted line corresponds to the frequency $f_{c}(L)$ of eddies of integral
length scale in size, while the dashed line corresponds to the frequency
$f_{c}(\eta)$ of eddies of Kolmogorov length scale in size, . Inset plot:
Compensated Fourier transform $\langle|\hat{p}|\rangle p^{-K}$ with $K=-4/3$.
The solid line corresponds to $\langle|\hat{p}|\rangle=3p^{K}$.

###  5.2 Pressure field on a sphere

To compare the pressure forcing with the effective forcing
$\mathcal{T}_{\ell}$, we interpolate the pressure field $p_{S}(\theta,\phi)$
in the single phase DNS on a sphere of radius $R_{0}$, and compute its
spherical harmonics decomposition

$p_{S}(\theta,\phi)=P_{c}\big{[}\sum_{\ell=0}^{\infty}\sum_{m=-\ell}^{\ell}P_{\ell,m}(t)\,Y_{\ell}^{m}(\theta,\phi)\big{]}.$
(34)

Similarly to the modes of deformation $x_{\ell,m}$, the statistics of
$P_{\ell,m}$ are independent of $m$. Ensemble averages are then computed over
the three simulations and the $m$ values. For pressure, the modes $\ell=0$ and
$\ell=1$ are non zero, however we focus in the following on modes $\ell\geq 2$
which are relevant for bubble deformations. Figure 13a shows that the standard
deviation of each mode $\ell$, $\sigma_{P}^{\ell}$, varies exponentially with
$\ell$. A higher $\ell$ is associated to fluctuations at a smaller scale,
which are known to be less energetic. However we have no explanation for the
exponential scaling. We also observed a decay of $\sigma_{\mathcal{T}}^{\ell}$
with $\ell$ (figure 5). The symmetry between positive and negative values is
restored, as shown on figure 13b. Distributions now show exponential tails for
both negative and positive pressure values. The shape of the distribution is
found independent of $\ell$, corresponding to the same hyperbolic secant
distribution (eq. (24)) than the effective forcing distribution we previously
identified.

Figure 13: a) Pressure standard deviation, $\sigma_{P}^{\ell}$ as a function
of $\ell$, showing an exponential decay with $\ell$ (black dotted line) b)
Distributions of $P_{\ell}$, normalized by $\sigma_{P}^{\ell}$, as a function
of $\ell$. All the pressure modes share the forcing distribution given in
equation (24).

Eventually, we compute the temporal Fourier transform $\hat{P}_{\ell}$ of the
spherical pressure modes $P_{\ell}$. Figure 14a shows the ensemble average of
the norm, $\langle|\hat{P}_{\ell}|\rangle$ as a function of the frequency. For
each $\ell$, we recover the three regimes we observed for the point pressure
spectrum and $\mathcal{T}_{\ell}$. The transition between the two first
regimes depends on the mode $\ell$. Considering that the pressure spectrum
share the same characteristic frequency than the effective forcing spectrum,
we expect the transition to occur at $f=\ell^{2/3}$, the frequency associated
with eddies of size $d/\ell$, in units of $t_{c}(d)$. We show in figure 14b
the spectra $\langle|\hat{P}_{\ell}|\rangle$ normalized by their low frequency
limit, $\hat{P}_{\ell}^{0}$, as a function of the frequency normalized by
$\ell^{2/3}$, the eddy turnover time at scale $d/\ell$. All curves collapse on
a single master curve, showing that pressure and effective forcing share the
same time scales. Below the critical frequency ($f<f_{\ell}$), the spectrum
amplitude converges to a constant value, significantly above the integral
frequency $f_{L}$. Similarly to Abry et al. (1994), the pressure spectrum at
low frequency is now constant. We can assume that the averaging over the
sphere has filtered the contribution from localized structures, and in
particular the vorticity filaments. A flat spectrum in the range
$f_{c}(L)<f<\ell^{2/3}$ also indicates that the contribution of eddies larger
than $d/\ell$, which are roughly homogeneous at the mode scale, has also been
filtered out: a bubble is mainly deformed by eddies at its scale. For
$\ell^{2/3}<f<f_{\eta}$, $\langle|\hat{P}_{\ell}|\rangle$ follows $f^{-3}$.
This decay is steeper than the $\ell$-dependency of
$\langle|\hat{\mathcal{T}}_{\ell}|\rangle$ which follows $f^{-2}$. This might
be attributed to the discrepancy between Eulerian and Lagrangian statistics.
From sweeping effect (Kraichnan, 1964), the temporal decorrelation of Eulerian
quantities are expected to occur faster than their Lagrangian counterpart.

To summarize, we have shown that the effective forcing $\mathcal{T}_{\ell}$
deforming a bubble inherits the pdf of pressure modes integrated over a sphere
of same radius. As a consequence of the filtering effect induced by the
integration over a sphere, the characteristic frequency associated to each
mode $\ell$ is the eddy turnover time at scale $d/\ell$, the frequencies
smaller than $\ell^{2/3}$ are well described by a white noise, and the forcing
amplitude decreases with $\ell$.

Figure 14: a) Amplitude of the pressure Fourier transform
$\langle|\hat{P}_{\ell}|\rangle$ for each mode $\ell$ as a function of the
frequency in unit of the eddy turn over time at the bubble scale. The black
dashed line is represents the eddy turnover time at scale $\eta$. The black
dotted line is the eddy turnover time at the integral length scale. b)
Normalised pressure Fourier transform as a function of frequency in unit of
the eddy turnover time at scale $d/\ell$. The black line follows $f^{-3}$.

###  5.3 Dissipation profiles

Our analysis of bubble deformation show that (i) The effective forcing
originates from pressure fluctuations near the bubble, and it is not affected
by bubble deformability. (ii) The damping of bubble oscillations is
significantly enhanced compare to the quiescent case. We interpreted this
additional damping by the presence of a boundary layer of size $L_{t}=R_{0}/5$
independent of We. In this section we investigate the velocity gradient
profile near the bubble, on a distance comparable to bubble typical
deformation. To do so, we need to compute a local distance $r$ to the
interface, which is not provided by the Basilisk VOF algorithm.

The method principle is the following. For every bulk point, we find the
closest grid point on the interface. We then interpolate the bubble surface
around this point, using a quadratic interpolation on the $20$ closest
neighbouring interfacial points. To find the neighbours efficiently, the
interfacial grid points are stored in a k-d tree structure. The distance $r$
to the interface is then found by minimizing the distance from the bulk point
to the quadratic manifold. We follow this procedure for both outside ($r>0$)
and inside ($r<0$) bulk points.

We diagnose the additional dissipative term of the linear model by
investigating the local dissipation rate profile around the bubble. The energy
dissipation rate per unit of mass in a elementary volume is related to the
local velocity gradients by

$\langle\epsilon\rangle(x)=2\nu\langle(\partial_{i}u_{j}+\partial_{j}u_{i})^{2}\rangle,$
(35)

where we use Einstein notations. For each run, we output snapshots of the full
flow field at times separated by at least one eddy turnover time at the bubble
scale, to ensure statistical independence. We then compute profiles of the
local dissipation near the interface by averaging on shells of constant
distance from the bubble interface, as illustrated on figure 15. Eventually,
for each Weber number, we ensemble average the flow snapshots (see table 3) to
extract a mean profile.

Figure 16a shows the average local dissipation, divided by the kinematic
viscosity, $\langle\epsilon\rangle(r)/\nu$, as a function of the distance $r$
to the bubble interface. Far from the bubble interface, for $r>R_{0}/2$ and
$r<-R_{0}/2$, the local dissipation converges to a constant. In the gas,
velocity gradients are maximum at $r=-R_{0}/15$. The existence of a maximum
inside the bubble near the interface originates from the nearly no slip
boundary condition imposed by the denser fluid on the gas inside the bubble.
Similar boundary layer has indeed been observed near solid particle surface
(no slip boundary condition) (Shen et al., 2022; Chiarini & Rosti, 2024). For
bubbles, we therefore expect that decreasing the gas density increases the
amplitude of the peak. The velocity gradients inside and outside the bubble
share the same order of magnitude: the dissipation hence mainly takes place
outside the bubble, in the liquid, where the dynamical viscosity is much
larger. To understand the origin of the additional dissipation we then focus
on the outside boundary layer.

We | 2 | 1.43 | 1 | 0.71 | 0.46 | 0.36 | 0.27
---|---|---|---|---|---|---|---
N | 48 | 68 | 68 | 27 | 46 | 24 | 52
Table 3: Number of snapshots per Weber number used to compute the flow
profiles. Figure 15: Example of the distance computation on a slice of a
bubble at $\textrm{We}=2$. The red points are on the interface. Isocontours
are separated by $0.0625R_{0}$. We also show the association between one bulk
points and the corresponding interfacial point.

For $r>0$, we observe a thick boundary layer of typical size $R_{0}/5$ (see
figure 16b), compatible with our estimation of $L_{t}$ (see equation (21)).
Figure 16b shows the dissipation rate value at the interface, in the liquid
$\langle\epsilon\rangle|_{0^{+}}$ as a function of the Weber number. At
vanishing Weber number, we find a non zero dissipation originating from a
geometrical boundary layer. The interfacial value varies between 3 times
($\textrm{We}=0.27$) and four times ($\textrm{We}=2$) larger than in the bulk.
In addition, we find an increase compatible with a linear dependency of the
interfacial dissipation with We. If we interpret this additional dissipation
as an energy transfer rate from the surface deformation to the flow, it would
scale as $\Lambda_{\ell}(\dot{x}_{\ell})^{2}$. We have
$(\dot{x}_{\ell})^{2}\sim(\omega_{\ell}\sigma_{x}^{\ell})^{2}\propto\textrm{We}$.
This interpretation is therefore compatible with a damping coefficient
$\Lambda$ independent of We.

In the absence of flow, the thickness of the boundary layer of the oscillating
bubble can be estimated by $\sqrt{2\nu/\omega_{2}}$. For a Weber number
ranging from 2.9 to 0.27, this estimation gives a boundary layer of size
ranging from $0.07R_{0}$ to $0.04R_{0}$, which is much thinner than the
boundary layer thickness we measured. We conclude that the boundary layer
originates from a geometrical turbulent boundary layer, and not from bubble
oscillations. The existence of a thick boundary layer was completely
disregarded in the computation of Lamb (1932) for a potential flow far from
the interface. The thick boundary layer we observed for the dissipation
profile is then consistent with a viscous damping one order of magnitude
larger than in the quiescent case.

Figure 16: a) Local velocity gradient inside and outside the bubble as a
function of the distance to the interface, for all We. b) Limit of the
dissipation rate at the bubble interface in the liquid phase as a function of
the Weber number.

##  6 Conclusion

In summary, we have shown that the deformations of bubbles in turbulence can
be described in terms of a stochastic linear oscillator on the Rayleigh modes
of oscillations up to a Weber of order unity. Conversely to previous works, we
have directly measured using DNS of bubbles in turbulence the coefficients of
this reduced model, namely, the damping rate and the natural frequency,
together with the statistical properties of the effective forcing. We have
shown that the natural frequency associated to each mode of deformation is not
modified compared to the quiescent case. For the effective dissipation, we
found that the damping is one order of magnitude larger than the prediction
from Lamb. Looking at the dissipation profiles near the interface, we
confirmed that the additional dissipation originates from a thick geometrical
boundary layer of size $L_{t}\approx R_{0}/5$ in our case. In physical units,
we expect the damping coefficients $\Lambda_{\ell}$ to scale as
$\nu/d^{2}\real(d)$. Eventually, we found that the effective forcing, which
results from the integration of pressure modes on the bubble surface, does not
depend on the Weber number. This observation confirms that bubble deformations
are one-way coupled to the flow: the back-reaction of bubble deformations on
the surrounding turbulent flow can be neglected. This effective forcing is
characterized by a probability distribution with exponential tails and a
typical correlation time which scales with the eddy turnover time at the
mode’s scale $t_{c}(d/\ell)$. We also looked at the statistics of pressure
fluctuations on a sphere in the absence of bubbles, and we showed that the
effective forcing shares the same pdf as the pressure modes’ pdf as well as
the same characteristic timescale. Due to the enhanced dissipation compared to
the quiescent case, we showed that the resonant oscillation mechanism is not
statistically relevant to explain break-ups. Indeed, at Weber of order unity,
the bubble cannot accumulate deformation energy on several periods of
oscillations as the quality factor $Q=\omega_{2}/\lambda$ of the main bubble
oscillations is too small. Instead, the critical Weber number is set by the
interaction between the statistics of large pressure fluctuations and bubble
deformations.

###### Acknowledgements.

We thank Adrian Van Kan, Luc Deike and François Pétrélis for scientific
discussions. We also thank Christophe Josserand and Laurent Duchemin 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 Ile 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. This work was
funded by a PSL Starting Grant 2022, as well as the ANR Lascaturb (reference
ANR-23-CE30-0043-03). Declaration of interests. The authors declare no
conflict of interest.

## References

  * Abry et al. (1994) Abry, Patrice, Fauve, S, Flandrin, P & Laroche, C 1994 Analysis of pressure fluctuations in swirling turbulent flows. Journal de Physique II 4 (5), 725–733.
  * Abu-Al-Saud et al. (2018) Abu-Al-Saud, Moataz O., Popinet, Stéphane & Tchelepi, Hamdi A. 2018 A conservative and well-balanced surface tension model. Journal of Computational Physics .
  * Bappy et al. (2019) Bappy, Mehedi, Carrica, Pablo M & Buscaglia, Gustavo C 2019 Lagrangian statistics of pressure fluctuation events in homogeneous isotropic turbulence. Physics of Fluids 31 (8).
  * Bappy et al. (2020a) Bappy, Mehedi H, Carrica, Pablo M, Vela-Martín, Alberto, Freire, Livia S & Buscaglia, Gustavo C 2020a Pressure statistics of gas nuclei in homogeneous isotropic turbulence with an application to cavitation inception. Physics of Fluids 32 (9).
  * Bappy et al. (2020b) Bappy, Mehedi Hasan, Vela-Martin, Alberto, Buscaglia, Gustavo Carlos, Carrica, Pablo M & Freire, Lívia Souza 2020b Effect of bubble size on lagrangian pressure statistics in homogeneous isotropic turbulence. In Journal of Physics: Conference Series, , vol. 1522, p. 012002. IOP Publishing.
  * Berry et al. (2005) Berry, SR, Hyers, RW, Racz, LM & Abedian, B 2005 Surface oscillations of an electromagnetically levitated droplet. International journal of thermophysics 26, 1565–1581.
  * Bojarevics & Pericleous (2003) Bojarevics, Valdis & Pericleous, Koulis 2003 Modelling electromagnetically levitated liquid droplet oscillations. ISIJ international 43 (6), 890–898.
  * Boussinesq (1877) Boussinesq, Joseph 1877 Essai sur la théorie des eaux courantes. Impr. nationale.
  * Brouzet et al. (2021) Brouzet, Christophe, Guiné, Raphaël, Dalbe, Marie-Julie, Favier, Benjamin, Vandenberghe, Nicolas, Villermaux, Emmanuel & Verhille, Gautier 2021 Laboratory model for plastic fragmentation in the turbulent ocean. Physical Review Fluids 6 (2), 024601.
  * Cadot et al. (1995) Cadot, O, Douady, S & Couder, Y 1995 Characterization of the low-pressure filaments in a three-dimensional turbulent shear flow. Physics of Fluids 7 (3), 630–646.
  * Calzavarini et al. (2009) Calzavarini, Enrico, Volk, Romain, Bourgoin, Mickaël, Lévêque, Emmanuel, Pinton, J-F & Toschi, Federico 2009 Acceleration statistics of finite-sized particles in turbulent flow: the role of faxén forces. Journal of Fluid Mechanics 630, 179–189.
  * Cao et al. (1999) Cao, Nianzheng, Chen, Shiyi & Doolen, Gary D 1999 Statistics and structures of pressure in isotropic turbulence. Physics of Fluids 11 (8), 2235–2250.
  * Chandrasekhar (1959) Chandrasekhar, Subrahmanyan 1959 The oscillations of a viscous liquid globe. Proceedings of the London Mathematical Society 3 (1), 141–149.
  * Chandrasekhar (2013) Chandrasekhar, Subrahmanyan 2013 Hydrodynamic and hydromagnetic stability. Courier Corporation.
  * Chiarini & Rosti (2024) Chiarini, Alessandro & Rosti, Marco Edoardo 2024 Finite-size inertial spherical particles in turbulence. arXiv preprint arXiv:2404.16475 .
  * De Langre (2008) De Langre, Emmanuel 2008 Effects of wind on plants. Annu. Rev. Fluid Mech. 40, 141–168.
  * Deike (2022) Deike, Luc 2022 Mass transfer at the ocean–atmosphere interface: the role of wave breaking, droplets, and bubbles. Annual Review of Fluid Mechanics 54, 191–224.
  * Douady et al. (1991) Douady, S, Couder, Y & Brachet, ME 1991 Direct observation of the intermittency of intense vorticity filaments in turbulence. Physical review letters 67 (8), 983.
  * Fan et al. (2024) Fan, Yaning, Wang, Cheng, Jiang, Linfeng, Sun, Chao & Calzavarini, Enrico 2024 Accelerations of large inertial particles in turbulence. Europhysics Letters .
  * Fauve et al. (1993) Fauve, Stephan, Laroche, C & Castaing, B 1993 Pressure fluctuations in swirling turbulent flows. Journal de physique II 3 (3), 271–278.
  * Galinat et al. (2007) Galinat, Sophie, Risso, Frédéric, Masbernat, Olivier & Guiraud, Pascal 2007 Dynamics of drop breakup in inhomogeneous turbulence at various volume fractions. Journal of Fluid Mechanics 578, 85–94.
  * Gitterman (2005) Gitterman, Moshe 2005 Noisy Oscillator, The: The First Hundred Years, From Einstein Until Now. World Scientific.
  * Gotoh & Rogallo (1999) Gotoh, Toshiyuki & Rogallo, Robert S 1999 Intermittency and scaling of pressure at small scales in forced isotropic turbulence. Journal of Fluid Mechanics 396, 257–285.
  * Håkansson (2019) Håkansson, Andreas 2019 Emulsion formation by homogenization: Current understanding and future perspectives. Annual review of food science and technology 10, 239–258.
  * Håkansson (2021) Håkansson, Andreas 2021 The role of stochastic time-variations in turbulent stresses when predicting drop breakup—a review of modelling approaches. Processes 9 (11), 1904.
  * Hinze (1955) Hinze, Julius O 1955 Fundamentals of the hydrodynamic mechanism of splitting in dispersion processes. AIChE journal 1 (3), 289–295.
  * Homann & Bec (2010) Homann, Holger & Bec, Jeremie 2010 Finite-size effects in the dynamics of neutrally buoyant particles in turbulent flow. Journal of Fluid Mechanics 651, 81–91.
  * Ishihara et al. (2003) Ishihara, Takashi, Kaneda, Yukio, Yokokawa, Mitsuo, Itakura, Ken’ichi & Uno, Atsuya 2003 Spectra of energy dissipation, enstrophy and pressure by high-resolution direct numerical simulations of turbulence in a periodic box. journal of the physical society of japan 72 (5), 983–986.
  * Kang & Leal (1988) Kang, IS & Leal, LG 1988 Small-amplitude perturbations of shape for a nearly spherical bubble in an inviscid straining flow (steady shapes and oscillatory motion). Journal of Fluid Mechanics 187, 231–266.
  * Kolmogorov (1949) Kolmogorov, Andrey 1949 On the breakage of drops in a turbulent flow. In Dokl. Akad. Navk. SSSR, , vol. 66, pp. 825–828.
  * Kraichnan (1964) Kraichnan, Robert H 1964 Kolmogorov’s hypotheses and eulerian turbulence theory. Phys. Fluids 7 (11), 1723–1734.
  * Lalanne et al. (2019) Lalanne, Benjamin, Masbernat, Olivier & Risso, Frédéric 2019 A model for drop and bubble breakup frequency based on turbulence spectra. AIChE Journal 65 (1), 347–359.
  * Lamb (1932) Lamb, Horace 1932 Hydrodynamics. Aufl., Cambridge: Univ. Press1879–1932 427.
  * Lee et al. (1987) Lee, Chung-Hur, Erickson, LE & Glasgow, LA 1987 Bubble breakup and coalescence in turbulent gas-liquid dispersions. Chemical Engineering Communications 59 (1-6), 65–84.
  * Luo & Svendsen (1996) Luo, Hean & Svendsen, Hallvard F 1996 Theoretical model for drop and bubble breakup in turbulent dispersions. AIChE journal 42 (5), 1225–1233.
  * Maniero et al. (2012) Maniero, Riccardo, Masbernat, Olivier, Climent, Eric & Risso, Frédéric 2012 Modeling and simulation of inertial drop break-up in a turbulent pipe flow downstream of a restriction. International journal of multiphase flow 42, 1–8.
  * Martínez-Bazán et al. (1999) Martínez-Bazán, Carlos, Montanes, JL & Lasheras, Juan C 1999 On the breakup of an air bubble injected into a fully developed turbulent flow. part 1. breakup frequency. Journal of Fluid Mechanics 401, 157–182.
  * Masuk et al. (2021a) Masuk, Ashik Ullah Mohammad, Qi, Yinghe, Salibindla, Ashwanth KR & Ni, Rui 2021a Towards a phenomenological model on the deformation and orientation dynamics of finite-sized bubbles in both quiescent and turbulent media. Journal of Fluid Mechanics 920, A4.
  * Masuk et al. (2021b) Masuk, Ashik Ullah Mohammad, Salibindla, Ashwanth KR & Ni, Rui 2021b Simultaneous measurements of deforming hinze-scale bubbles with surrounding turbulence. Journal of Fluid Mechanics 910.
  * Maxey & Riley (1983) Maxey, Martin R & Riley, James J 1983 Equation of motion for a small rigid sphere in a nonuniform flow. The Physics of Fluids 26 (4), 883–889.
  * Miller & Scriven (1968) Miller, CA & Scriven, LE 1968 The oscillations of a fluid droplet immersed in another fluid. Journal of fluid mechanics 32 (3), 417–435.
  * Mordant et al. (2004) Mordant, Nicolas, Crawford, Alice M & Bodenschatz, Eberhard 2004 Three-dimensional structure of the lagrangian acceleration in turbulent flows. Physical review letters 93 (21), 214501.
  * Perrard et al. (2021) Perrard, Stéphane, Rivière, Aliénor, Mostert, Wouter & Deike, Luc 2021 Bubble deformation by a turbulent flow. Journal of Fluid Mechanics 920, A15.
  * Pope (2000) Pope, S. B. 2000 Turbulent flows. Cambridge university press, Cambridge, UK.
  * Popinet (2003) Popinet, S. 2003 Gerris: a tree-based adaptive solver for the incompressible euler equations in complex geometries. J. Comput. Phys. 190 (2), 572–600.
  * Popinet (2009) Popinet, S. 2009 An accurate adaptive solver for surface-tension-driven interfacial flows. Journal of Computational Physics 228, 5838–5866.
  * Prakash et al. (2012) Prakash, Vivek N, Tagawa, Yoshiyuki, Calzavarini, Enrico, Mercado, Julián Martínez, Toschi, Federico, Lohse, Detlef & Sun, Chao 2012 How gravity and size affect the acceleration statistics of bubbles in turbulence. New journal of physics 14 (10), 105017.
  * Prandtl (1949) Prandtl, Ludwig 1949 Report on investigation of developed turbulence. Zeitschrift fuer Angewandte Matematik und Mechanik 5 (NACA-TM-1231).
  * Prosperetti (1977) Prosperetti, Andrea 1977 Viscous effects on perturbed spherical flows. Quarterly of Applied mathematics 34 (4), 339–352.
  * Prosperetti (1980) Prosperetti, Andrea 1980 Free oscillations of drops and bubbles: the initial-value problem. Journal of Fluid Mechanics 100 (2), 333–347.
  * Pullin & Rogallo (1994) Pullin, DI & Rogallo, RS 1994 Pressure and higher-order spectra for homogeneous isotropic turbulence. Stanford Univ., Studying Turbulence Using Numerical Simulation Databases. 5: Proceedings of the 1994 Summer Program .
  * Pumir (1994) Pumir, Alain 1994 A numerical study of pressure fluctuations in three-dimensional, incompressible, homogeneous, isotropic turbulence. Physics of Fluids 6 (6), 2071–2083.
  * Qureshi et al. (2007) Qureshi, Nauman M, Bourgoin, Mickaël, Baudet, Christophe, Cartellier, Alain & Gagne, Yves 2007 Turbulent transport of material particles: an experimental study of finite size effects. Physical review letters 99 (18), 184502\.
  * Ravelet et al. (2011) Ravelet, Florent, Colin, Catherine & Risso, Frédéric 2011 On the dynamics and breakup of a bubble rising in a turbulent flow. Physics of Fluids 23 (10).
  * Rayleigh (1879) Rayleigh, Lord 1879 On the capillary phenomena of jets. Proc. R. Soc. London 29 (196-199), 71–97.
  * Reid (1960) Reid, William Hill 1960 The oscillations of a viscous liquid drop. Quarterly of Applied Mathematics 18 (1), 86–89.
  * Risso (2018) Risso, Frédéric 2018 Agitation, mixing, and transfers induced by bubbles. Annual Review of Fluid Mechanics 50, 25–48.
  * Risso & Fabre (1998) Risso, Frédéric & Fabre, Jean 1998 Oscillations and breakup of a bubble immersed in a turbulent field. Journal of Fluid Mechanics 372, 323–355.
  * Rivière et al. (2023) Rivière, Aliénor, Duchemin, Laurent, Josserand, Christophe & Perrard, Stéphane 2023 Bubble breakup reduced to a one-dimensional nonlinear oscillator. Physical Review Fluids 8 (9), 094004.
  * Rivière et al. (2021) Rivière, Aliénor, Mostert, Wouter, Perrard, Stéphane & Deike, Luc 2021 Sub-hinze scale bubble production in turbulent bubble break-up. Journal of Fluid Mechanics 917, A40.
  * Roa et al. (2023) Roa, Ignacio, Renoult, Marie-Charlotte, Dumouchel, Christophe & Brändle de Motta, Jorge César 2023 Droplet oscillations in a turbulent flow. Frontiers in Physics 11, 1173521.
  * Rosales & Meneveau (2005) Rosales, Carlos & Meneveau, Charles 2005 Linear forcing in numerical simulations of isotropic turbulence: Physical space implementations and convergence properties. Physics of fluids 17 (9).
  * Rosti et al. (2018) Rosti, Marco Edoardo, Banaei, Arash Alizad, Brandt, Luca & Mazzino, Andrea 2018 Flexible fiber reveals the two-point statistical properties of turbulence. Physical review letters 121 (4), 044501.
  * Ruth et al. (2019) Ruth, Daniel J, Mostert, Wouter, Perrard, Stéphane & Deike, Luc 2019 Bubble pinch-off in turbulence. Proceedings of the National Academy of Sciences 116 (51), 25412–25417.
  * Salibindla et al. (2021) Salibindla, Ashwanth KR, Masuk, Ashik Ullah Mohammad & Ni, Rui 2021 Experimental investigation of the acceleration statistics and added-mass force of deformable bubbles in intense turbulence. Journal of Fluid Mechanics 912, A50.
  * Sevik & Park (1973) Sevik, M. & Park, S. H. 1973 The Splitting of Drops and Bubbles by Turbulent Fluid Flow. Journal of Fluids Engineering 95 (1), 53–60.
  * Shen et al. (2022) Shen, Jie, Peng, Cheng, Wu, Jianzhao, Chong, Kai Leong, Lu, Zhiming & Wang, Lian-Ping 2022 Turbulence modulation by finite-size particles of different diameters and particle–fluid density ratios in homogeneous isotropic turbulence. Journal of Turbulence 23 (8), 433–453.
  * Tennekes (1975) Tennekes, Henk 1975 Eulerian and lagrangian time microscales in isotropic turbulence. Journal of Fluid Mechanics 67 (3), 561–567.
  * Toschi & Bodenschatz (2009) Toschi, Federico & Bodenschatz, Eberhard 2009 Lagrangian properties of particles in turbulence. Annual review of fluid mechanics 41, 375–404.
  * Vedula & Yeung (1999) Vedula, Prakash & Yeung, Pui-Kuen 1999 Similarity scaling of acceleration and pressure statistics in numerical simulations of isotropic turbulence. Physics of Fluids 11 (5), 1208–1220.
  * Vela-Martín & Avila (2021) Vela-Martín, Alberto & Avila, Marc 2021 Deformation of drops by outer eddies in turbulence. Journal of Fluid Mechanics 929, A38.
  * Verhille (2022) Verhille, Gautier 2022 Deformability of discs in turbulence. Journal of Fluid Mechanics 933, A3.
  * Volk et al. (2011) Volk, Romain, Calzavarini, Enrico, Leveque, Emmanuel & Pinton, J-F 2011 Dynamics of inertial particles in a turbulent von kármán flow. Journal of Fluid Mechanics 668, 223–235.
  * Volk et al. (2008) Volk, R, Calzavarini, Enrico, Verhille, G, Lohse, Detlef, Mordant, N, Pinton, J-F & Toschi, Federico 2008 Acceleration of heavy and light particles in turbulence: comparison between experiments and direct numerical simulations. Physica D: Nonlinear Phenomena 237 (14-17), 2084–2089.
  * Voth et al. (2002) Voth, Greg A, La Porta, Arthur, Crawford, Alice M, Alexander, Jim & Bodenschatz, Eberhard 2002 Measurement of particle accelerations in fully developed turbulence. Journal of Fluid Mechanics 469, 121–160.
  * Wang et al. (2003) Wang, Tiefeng, Wang, Jinfu & Jin, Yong 2003 A novel theoretical breakup kernel function for bubbles/droplets in a turbulent flow. Chemical Engineering Science 58 (20), 4629–4637.
  * Xiao et al. (2021) Xiao, Xiao, Brillo, Jürgen, Lee, Jonghyun, Hyers, Robert W & Matson, Douglas M 2021 Impact of convection on the damping of an oscillating droplet during viscosity measurement using the iss-eml facility. npj Microgravity 7 (1), 36.
  * Yao et al. (2008) Yao, Hua-Dong, He, Guo-Wei, Wang, Meng & Zhang, Xing 2008 Time correlations of pressure in isotropic turbulence. Physics of Fluids 20 (2).