Banner appropriate to article type will appear here in typeset article

# Scale-by-scale kinetic energy budgets in multiphase turbulence

F. Thiesset¹*, J. Vahé¹

¹CNRS, CORIA, UMR 6614, Normandy Univ., UNIROUEN, INSA Rouen, 76000 Rouen, France

(Received xx; revised xx; accepted xx)

The present work aims at exploring the scale-by-scale kinetic energy exchanges in multiphase turbulence. For this purpose, we derive the Kármán-Howarth-Monin equation which accounts for the variations of density and viscosity across the two phases together with the effect of surface tension. We consider both conventional and phase conditional averaging operators. This framework is applied to numerical data from detailed simulations of forced homogeneous and isotropic turbulence covering different values for the liquid volume fraction, the liquid/gas density ratio, the Reynolds, and Weber numbers. We confirm the existence of an additional transfer term due to surface tension. Part of the kinetic energy injected at large scales is transferred into kinetic energy at smaller scales by classical non-linear transport while another part is transferred to surface energy before being released back into kinetic energy, but at smaller scales. The overall kinetic energy transfer rate is larger than in single phase flows. Kinetic energy budgets conditioned in a given phase show that the scale-by-scale transport of turbulent kinetic energy due to pressure is a gain (loss) of kinetic energy for the lighter (heavier) phase. Its contribution can be dominant when the gas volume fraction becomes small or when the density ratio increases. Building on previous work, we hypothesize the existence of a pivotal scale above which kinetic energy is stored into surface deformation and below which the kinetic energy is released by interface restoration. Some phenomenological predictions for this scale are discussed.

**Key words:** Turbulent flow, Multiphase flow, Turbulence theory.

## 1. Introduction

Turbulence is a state of fluids characterized by some seemingly erratic motions, although not Brownian, of many sized, yet hierarchized, eddies. It is ubiquitous in our everyday life and recognized as one of the most important open problems of modern physics. The study of single fluid turbulence has concentrated most of the researchers' efforts over the last century and remains even nowadays a field of intense research.

<footnote>* Email address for correspondence: fabien.thiesset@cnrs.fr</footnote>
The case of turbulent flows involving a mixture of immiscible fluids separated by an interface have received much less attention.

Most of available fundamental knowledge on multiphase fluid turbulence applies to the dispersed regime (Balachandar & Eaton 2010), i.e. to cases where one phase, the carrier phase, is much more abundant than the other, the dispersed phase. When the minority phase represents more than say 5% of total volume, we talk about multiphase fluid turbulence in the dense regime. Investigating this regime of multiphase flows has been made possible only recently with the aid of high fidelity numerical simulations (see e.g. Duret et al. 2012; Dodd & Ferrante 2016; Lu & Tryggvason 2018, 2019; Rosti et al. 2019; Crialesi-Esposito et al. 2022, among others).

Multiphase turbulence distinguishes from its single-phase counterpart by several aspects. First, in most situations, the two phases have different density and different viscosity. Despite fundamental differences in their physical mechanisms, multiphase turbulence shares key features with compressible, reacting, and gas mixture flows, notably spatial and temporal variations in density and viscosity. There is potential for valuable cross-fertilization, particularly through the use of analytical tools that have been independently developed for these distinct classes of flow. Density/viscosity variations are known to influence the whole flow dynamics and in particular the transfer of kinetic energy between the different scales of the turbulent spectrum (see e.g. Taguelmimt et al. 2016; Gauding et al. 2018; Lai et al. 2018; Galtier & Banerjee 2011; Aluie 2011; Eyink & Drivas 2018; Whitman et al. 2019; Sabelnikov et al. 2019a,b; Kolla et al. 2014). Secondly, the other salient feature of multiphase turbulence is the presence of interfacial forces associated to surface tension. The presence of the interface couples the interface geometry to the flow dynamics, and the interface may either store energy or release energy from/to the rest of the domain (e.g. Dodd & Ferrante 2016; Treffzt-Posada & Ferrante 2023). The interface is also known to increase internal intermittency (Crialesi-Esposito et al. 2023a).

As for its single-phase counterpart, multiphase turbulence gives rise to a wide and continuous spectrum of eddies. According to the classical picture of the turbulent cascade of Richardson (1922), velocity fluctuations are injected at large scales, then transferred between scales of decreasing sizes, down to the small scales where energy is dissipated. The greatest challenge in the fluid turbulence community is thus to describe how the velocity fluctuations vary with eddy size. This requires what is referred to as a scale-by-scale theory that explicitly specifies that notion of scale. This is traditionally done using Fourier transforms (Lance & Bataille 1991; Bunner & Tryggvason 2002; Dodd & Ferrante 2016; Risso 2018; Hellinger et al. 2021a,b; Crialesi-Esposito et al. 2022; Innocenti et al. 2021; Zamansky et al. 2024; Ramirez et al. 2024), structure (or correlation) functions (Galtier & Banerjee 2011; Banerjee & Galtier 2013; Lai et al. 2018; Gauding et al. 2018; Ferrand et al. 2020; Hellinger et al. 2021a,b; Kritsuk et al. 2013; Wagner et al. 2012; Bunner & Tryggvason 2002; Lu & Tryggvason 2018, 2019; Trautner et al. 2021), coarse-grained approaches (Aluie 2011, 2013; Eyink & Drivas 2018; Wang et al. 2018; Pandey et al. 2020; Innocenti et al. 2021; Pandey et al. 2023; Narula et al. 2025), wavelet decomposition (Freund & Ferrante 2019). Note that the above cited literature is not exhaustive and relates to either multiphase flows or to flows with density/viscosity variations.

As far as multiphase flows are concerned, the aforementioned analyses of the scale-by-scale kinetic energy budgets have revealed some crucial departures from single-phase flows. One of the most important departures is that surface tension acts as an additional kinetic energy transfer between the different turbulent eddies. This has been shown for both liquid-gas flows (Pandey et al. 2020, 2023; Ramirez
et al. 2024) and liquid-liquid emulsions (Perlekar 2019; Crialesi-Esposito et al. 2022, 2023b; Cannon et al. 2024). The interface was found to absorb kinetic energy from the large scales and release it at smaller scales. For turbulent emulsions, Perlekar et al. (2012); Crialesi-Esposito et al. (2022, 2023b); Cannon et al. (2024) explain this behaviour under the framework proposed by Kolmogorov (1949) and Hinze (1955) (see also Ni 2024). They suggest that the large turbulent eddies are intense enough to overcome the cohesive force due to surface tension, thereby leading to the fragmentation of the fluid parcels. Through this fragmentation process, the interface area increases, and kinetic energy is transferred into surface energy through the process identified by Dodd & Ferrante (2016). In this range of scales, the interface thus absorbs kinetic energy and stores it as surface energy. On the contrary, at smaller scales, the fluid parcels remain stable under the effect of surface tension. These parcels can however coalesce, which leads to a decrease of the interface area. In this range of scales, the surface energy is thus released into kinetic energy. As per Crialesi-Esposito et al. (2022, 2023b); Cannon et al. (2024), the pivotal scale between kinetic energy absorption and kinetic energy release is thus expected to correspond to the pivotal scale between breakup and coalescence which is nothing but the Kolmogorov-Hinze scale. However, there are situations with no break-up nor coalescence, in particular in gravity driven bubbly flows (Pandey et al. 2020, 2023), where the scale-by-scale contribution of the surface tension term exhibits a similar pivotal scale. Hence, the Kolmogorov-Hinze framework and its underlying mechanisms cannot explain the effect of surface tension in general situations. In bubbly flows, Pandey et al. (2020) interpret the role of the surface tension term as a mechanism by which large-scale kinetic energy is absorbed through bubble deformation and stretching. When relaxing, this stored energy is released, but primarily at smaller scales, effectively acting as a scale-conversion process mediated by surface tension. Building on this, we also hypothesize that the breakup/coalescence range could likely be replaced by the more general mechanisms of interface deformation at large scales and interface restoration at smaller scales. This idea originates from Perlekar (2019) who conjectured that the interface may alter the turbulent kinetic energy transfer similarly to some situations of viscoelastic fluid turbulence (see e.g. Valente et al. 2014; Nguyen et al. 2016) where polymers store kinetic energy when being stretched and release kinetic energy when relaxing, these two processes occurring in different ranges of scales. A similar idea was also conjectured by Pandey et al. (2020).

Many studies on multiphase turbulence focus on the case of equal density between the two phases (e.g. Mukherjee et al. 2019; Crialesi-Esposito et al. 2022; Cannon et al. 2024, among others). However, the presence of density difference leads to an additional contribution in the scale-by-scale kinetic energy budget term due to pressure (Pandey et al. 2020, 2023; Narula et al. 2025) even when the flow is homogeneous, which is generally referred to as the baropycnal work (Aluie 2013). Contrasts of density may also change the whole flow dynamics as the phase with higher density carries the flow momentum. The role played by density ratio remains yet only partly understood, and more dedicated studies are required.

Here, the general objective is to explore the scale-by-scale kinetic energy exchanges in multiphase turbulence by use of a detailed analysis of the structure function transport equation (known as the Kármán-Howarth-Monin equation, abbreviated KHM). Only very few studies of this type exist (Pan & Banerjee 2022, 2025; Narula et al. 2025) which motivates the present work. In particular, we aim at using this framework to unravel the effect of density and/or viscosity variations, together with surface tension on the kinetic energy scale distribution. Another objective of the
present study is to infer whether the turbulence properties taken in one given phase differ from classical single-phase constant-density, constant-viscosity turbulence. This requires conditioning the structure functions by the phase, in the lines of Lu & Tryggvason (2018, 2019); Trautner et al. (2021) and derive their associated budget equation in a somewhat similar fashion as Freund & Ferrante (2019) who used the wavelet approach. In other words, we aim at extending the conditionally averaged kinetic energy budgets derived by Dodd & Ferrante (2016); Trefftz-Posada & Ferrante (2023) to the two-point statistics. We will investigate the peculiar case of homogeneous and isotropic turbulence using numerical simulations data obtained by a standard front capturing code named archer.

The paper is organized as follows. First, the general KHM equation applying to multiphase flows is derived in section §2. We will consider both unconditional and conditional averaged equations. All terms of the KHM equation are also analysed in the asymptotic limit of large separations. The numerical setup is detailed in section §3. Results are presented in section §4 and conclusions are drawn in a last section §5. Technical aspects regarding the derivation of the KHM equation and the numerical simulations are gathered in appendices A and B, respectively.

## 2. Kármán-Howarth-Monin equation for multiphase turbulence

### 2.1. Derivation

We seek for a generalized Kármán-Howarth-Monin (KHM) equation (sometimes also referred to as the generalized Kolmogorov equation) that applies to two-phase flows. The latter should account for the variations of density and viscosity together with surface tension. As per Galtier & Banerjee (2011); Lai et al. (2018); Hellinger et al. (2021b) and reference therein, one possible definition of the scale-by-scale kinetic energy in variable-density flows is:

$$|\delta u|_{\rho}^{2} := (\delta(\rho u)) \cdot (\delta u), \qquad (2.1)$$

where $(\delta(\rho u)) := (\rho u)^{+} - (\rho u)^{-}$ and $(\delta u) := u^{+} - u^{-}$ are the increment of $\rho u$ and $u$, respectively, between two points $x^{+}$ and $x^{-}$ arbitrarily separated in space by a distance $r := x^{+} - x^{-}$. The superscript $+(-)$ is used to denote that quantities are taken at point $x^{+}(x^{-})$. It is also convenient to define $(\bar{\delta}_{\bullet}) := (\bullet^{+} + \bullet^{-})/2$, i.e. the arithmetic mean of $\bullet$ between point $x^{+}$ and $x^{-}$. Note that some alternative definitions for the scale-by-scale kinetic energy in variable-density flows can be found in the literature (Ferrand et al. 2020; Brahami 2020). This is discussed in more details in Appendix A.3 and in subsection §2.2.

The transport equation for $|\delta u|_{\rho}^{2}$ is obtained directly from the one-fluid formulation for the two-phase incompressible Navier-Stokes (NS) equation. Each term can be derived from the following general expression:

$$(\delta(\rho u)) \cdot (\delta T_u) + (\delta u) \cdot (\delta T_{\rho u}), \qquad (2.2)$$

where $T_u$ and $T_{\rho u}$ denote the transport equation for $u$ and $\rho u$, respectively. The one-fluid formulation of the two-phase incompressible NS-equation writes:

$$T_{\rho u} := \partial_t \rho u = -\nabla \cdot \rho u u - \nabla P + \nabla \cdot \mathbf{t} + S + F. \qquad (2.3)$$

The incompressibility condition further implies $\nabla \cdot u = 0$. In Eq. (2.3), the mechanical pressure is denoted $P$, $\mathbf{t} := \mu_2 \mathbf{S}$ is the viscous stress tensor with $\mathbf{S} := (\nabla u + \nabla u^T)/2$ the strain rate tensor, $\mu$ the dynamic viscosity. The term noted $F$ corresponds to a
generic forcing, while $S := 2\sigma H\delta_{\Gamma}(x - x_{\Gamma})n_{\Gamma}$ is the surface tension term, where $\sigma$ is the surface tension coefficient, $H$ is the mean curvature of the interface, $n_{\Gamma}$ is the normal to the interface. The surface tension term acts only at the interface, hence the presence of the Dirac $\delta_{\Gamma}(x - x_{\Gamma})$ function allowing the interface position $x_{\Gamma}$ to be tracked. Using the continuity equation, one then obtains the transport equation for $u$:

$$ T_u := \partial_t u = -(u \cdot \nabla)u - v\nabla P + v\nabla \cdot \mathbf{t} + vS + vF, \quad (2.4) $$

where $v := 1/\rho$ is the specific volume.

The quantity $|\delta u|_{\rho}^2$ is a fluctuating field variable that a priori depends on time $t$, the position $X := (x^+ + x^-)/2$ of the midpoint between $x^+$ and $x^-$ and the separation vector $r = x^+ - x^-$ (See Appendix A). The scale-by-scale kinetic energy $|\delta u|_{\rho}^2$ is better assessed in a statistical sense, taking advantage of the statistical symmetry of the flow (if any). Here, we first employ an ensemble average over many statistically independent samples of the flow. These averages will be noted $\langle \bullet \rangle_{\mathbb{E}}$, where $\bullet$ denotes any quantity (in practice, the flow being at statistically stationary state, time averages are computed in place of ensemble average, supposing ergodicity). This ensemble average is then supplemented by a spatial average over a set of points $\mathbb{C}$:

$$ \langle \bullet \rangle_{\mathbb{C}} := \frac{\int_{\mathbb{C}} \langle \bullet \rangle_{\mathbb{E}} dV}{\int_{\mathbb{C}} dV}. \quad (2.5) $$

We will consider four different sets of averaging points $\mathbb{C}$ as detailed below:

(i) The first average covers the ensemble $\mathbb{T}$ of all points of the simulation domain. The symbol $\mathbb{T}$ stands for “total”. This is for instance the classical volume average in Direct Numerical Simulation of the NS equations in a periodic box.

(ii) We also consider the subset $\mathbb{L}$ of the simulation domain which covers the points $X$ such as both $x^+$ and $x^-$ lie in the liquid phase, i.e. $\mathbb{L} := \{X : (x^+ \in \text{liquid}) \wedge (x^- \in \text{liquid})\}$. The symbol $\mathbb{L}$ stands for “liquid”.

(iii) The complementary subset for the gas phase is further considered, i.e. $\mathbb{G} := \{X : (x^+ \in \text{gas}) \wedge (x^- \in \text{gas})\}$. The symbol $\mathbb{G}$ stands for “gas”.

(iv) We finally define conditional averages over the set of points $\mathbb{M}$ such as $x^+$ and $x^-$ lie in different phases, i.e. $\mathbb{M} := \{X : (x^+ \in \text{liquid}) \wedge (x^- \in \text{gas})\}$ (and vice versa). The symbol $\mathbb{M}$ stands for “mixed”.

More details on conditional averages are given in Appendix A, see Eqs. (A 15). Note that for conditional averages the averaging volume depends on $r$. Its boundary, noted $\partial\mathbb{C}$, does not necessarily coincide with the liquid-gas interface.

Once spatially averaged, two-point statistics depend on time $t$ and the separation vector $r$. In case of isotropic flows, two-point statistics are only function of the modulus of the vector $r$, thereby reducing the problem complexity to 2-dimensions ($|r|, t$). When the flow is anisotropic, one could keep the dependence to $r$ (a 3-dimensional vector), thus paying the price of a 4-dimensional problem. The other solution to cope with anisotropy is to apply an angular average over all orientations of the separation vector $r$ (Nie & Tanveer 1999; Hill 2002).

$$ \langle \bullet \rangle_{\mathbb{C}, \Omega} := \frac{1}{4\pi} \iint_{\Omega} \langle \bullet \rangle_{\mathbb{C}} \sin \theta d\theta d\varphi, \quad (2.6) $$

where the set of solid angles is given by $\Omega = \{\varphi, \theta | 0 \le \varphi \le \pi, 0 \le \theta \le 2\pi\}$ with $\varphi = \arctan(r_y/r_x)$ and $\theta = \arccos(r_z/|r|)$ ($r_x, r_y$, and $r_z$ are the components of the
$\mathbf{r}$ vector in $x, y$, and $z$ directions, respectively). Since angularly averaged two-point statistics depend on $|\mathbf{r}|$ only (and not its orientation anymore), the effect of anisotropy is concealed. However, it has the advantage of lightening the analysis since $\langle \bullet \rangle_{\mathcal{C},\Omega}$ depends only on time $t$ and $r = |\mathbf{r}|$, i.e. a 2-dimensional problem, as for isotropic flows. In what follows, an angular average will be applied to the two-point statistics, although for the sake of clarity, the subscript $\Omega$ will be dropped from the notations. Note that the flow under consideration here is isotropic, and hence angular average are only used to increase the statistical convergence (Taylor et al. 2003; Thiesset et al. 2020; Gauding et al. 2022; Federrath et al. 2021).

The transport equation for $\langle |\delta\mathbf{u}|_{\rho}^2 \rangle_{\mathcal{C}}$ is referred to as the KHM equation. Its detailed derivation is given in Appendix A, and we summarize here the main results. The general KHM equation can be written formally as

$$d_t \langle |\delta \mathbf{u}|_{\rho}^2 \rangle_{\mathcal{C}} = \langle \mathcal{T} \rangle_{\mathcal{C}} + \langle \mathcal{P} \rangle_{\mathcal{C}} + \langle \mathcal{V} \rangle_{\mathcal{C}} + \langle \mathcal{S} \rangle_{\mathcal{C}} + \langle \mathcal{F} \rangle_{\mathcal{C}} \quad (2.7)$$

This equation is one of the key result of the present paper. It is a generalization of the KHM equation to account for variable density and the presence of surface tension. In this formulation, it applies irrespectively of the averaging volume $\mathcal{C} \in \{\mathcal{T}, \mathcal{L}, \mathcal{G}, \mathcal{M}\}$ which is another strength of the present work. The newly derived KHM equation contains different terms which are now detailed.

* $d_t \langle |\delta\mathbf{u}|_{\rho}^2 \rangle_{\mathcal{C}}$ denotes the time variations of $\langle |\delta\mathbf{u}|_{\rho}^2 \rangle_{\mathcal{C}}$. In statistically stationary flows, this term is zero.

* $\langle \mathcal{T} \rangle_{\mathcal{C}}$ represents the “transport” of $\langle |\delta\mathbf{u}|_{\rho}^2 \rangle_{\mathcal{C}}$. In incompressible flows, the latter can be decomposed into three contributions:

$$\langle \mathcal{T} \rangle_{\mathcal{C}} = -\langle \nabla_{\mathbf{X}} \cdot (\bar{\delta}\mathbf{u}) |\delta\mathbf{u}|_{\rho}^2 \rangle_{\mathcal{C}} - \langle \nabla_{\mathbf{r}} \cdot (\delta\mathbf{u}) |\delta\mathbf{u}|_{\rho}^2 \rangle_{\mathcal{C}} + \int_{\partial\mathcal{C}} |\delta\mathbf{u}|_{\rho}^2 \mathbf{u}_b \cdot \mathbf{n} \, dS. \quad (2.8)$$

The leftmost term on RHS of Eq. (2.8) is the averaged transfer of $|\delta\mathbf{u}|_{\rho}^2$ in flow position space $\mathbf{X}$. Using the divergence theorem, it can be rewritten as

$$-\langle \nabla_{\mathbf{X}} \cdot (\bar{\delta}\mathbf{u}) |\delta\mathbf{u}|_{\rho}^2 \rangle_{\mathcal{C}} = - \int_{\partial\mathcal{C}} |\delta\mathbf{u}|_{\rho}^2 (\bar{\delta}\mathbf{u}) \cdot \mathbf{n} \, dS. \quad (2.9)$$

This term can thus be interpreted as a flux of the quantity $|\delta\mathbf{u}|_{\rho}^2$ which passes through the averaging volume boundary $\partial\mathcal{C}$ with a velocity $\bar{\delta}\mathbf{u}$, with $\mathbf{n}$ being the outwardly pointing normal vector to $\partial\mathcal{C}$. If the averaging volume is a periodic box ($\mathcal{C} = \mathcal{T}$ is periodic), then this term is zero since all the fluxes coming in and out from the domain cancels out (Hill 2002). In contrast, when $\mathcal{C} \in \{\mathcal{L}, \mathcal{G}, \mathcal{M}\}$, the transfer in $\mathbf{X}$-space cannot be dropped out.

The second term on RHS of Eq. (2.8) represents a flux of $|\delta\mathbf{u}|_{\rho}^2$ in scale-space $\mathbf{r}$, thus corresponding to the transfer of kinetic energy among the different scales of the flow. The rightmost term on RHS of Eq. (2.8) arises when applying the Reynolds transport theorem to the volume-averaged time-derivative of $|\delta\mathbf{u}|_{\rho}^2$, viz.

$$\langle \partial_t |\delta \mathbf{u}|_{\rho}^2 \rangle_{\mathcal{C}} = d_t \langle |\delta \mathbf{u}|_{\rho}^2 \rangle_{\mathcal{C}} - \int_{\partial \mathcal{C}} |\delta \mathbf{u}|_{\rho}^2 \mathbf{u}_b \cdot \mathbf{n} \, dS. \quad (2.10)$$

Here, $\mathbf{u}_b$ is the velocity of the control volume boundary $\partial\mathcal{C}$. When the averaging volume is fixed, for instance when $\mathcal{C} = \mathcal{T}$, one has $\mathbf{u}_b \equiv 0$ and this term vanishes. Otherwise, it represents the averaged flux of $|\delta\mathbf{u}|_{\rho}^2$ which flows through the control volume boundary at a velocity $\mathbf{u}_b$.
• $\langle \mathcal{P} \rangle_{\mathcal{C}}$ represents the effect of pressure on the evolution of the scale-by-scale kinetic energy. In incompressible flows, it can be decomposed into two contributions:

$$ \langle \mathcal{P} \rangle_{\mathcal{C}} = -2 \langle \nabla_{\mathbf{X}} \cdot (\delta \mathbf{u}) (\delta P) \rangle_{\mathcal{C}} - \langle \mathcal{C}(-\nabla P) \rangle_{\mathcal{C}}, \quad (2.11) $$

The leftmost term on RHS of Eq. (2.11) is the scale-by-scale transport of turbulent kinetic energy due to pressure, hereafter simply referred to as pressure transport. It can be rewritten using the divergence theorem

$$ -2 \langle \nabla_{\mathbf{X}} \cdot (\delta \mathbf{u}) (\delta P) \rangle_{\mathcal{C}} = - \int_{\partial \mathcal{C}} (\delta P)(\delta \mathbf{u}) \cdot \mathbf{n} dS \quad (2.12) $$

This term accounts for the pressure transport of scale-by-scale kinetic energy through the averaging volume boundary $\partial \mathcal{C}$. Hence, if $\mathcal{C} = \mathcal{T}$ is periodic, this term is zero and the contribution due to pressure reduces to $\langle \mathcal{P} \rangle_{\mathcal{T}} = -\langle \mathcal{C}(-\nabla P) \rangle_{\mathcal{T}}$ (the meaning of term $\mathcal{C}$ will be described very shortly). It is also sometimes referred to as the scale-by-scale pressure diffusion as it acts in redistributing the scale-by-scale turbulent kinetic energy across position $\mathbf{X}$ and scales $r$.

The second contribution of pressure is due to the term $-\langle \mathcal{C}(-\nabla P) \rangle_{\mathcal{C}}$ which arises due to variations of the density (see Eq. (A5) in Appendix A). It is likely to be analogous to what Aluie (2013) identifies to the baropycnal work. By coarse-graining the compressible NS-equation, Aluie demonstrated that a term similar to $-\langle \mathcal{C}(-\nabla P) \rangle_{\mathcal{C}}$ should operate even in incompressible flows, provided there are density variations. The analyses by Pandey et al. (2020, 2023); Narula et al. (2025) further support the presence of such a contribution in the energy budget of two-phase incompressible flows. The presence of this term is discussed in more details in subsection §2.2.

Note that since the density is constant in either the liquid or gas phase, the term $-\langle \mathcal{C}(-\nabla P) \rangle_{\mathcal{C}}$ should be accounted for only if $\mathcal{C} \in \{\mathcal{T}, \mathcal{M}\}$ and if the density is different in the liquid and gas phase. In summary, when $\mathcal{C} \equiv \mathcal{T}$, only the term $-\langle \mathcal{C}(-\nabla P) \rangle_{\mathcal{C}}$ contributes to the budget since the pressure transport is zero by homogeneity. In contrast, when $\mathcal{C} \in \{\mathcal{L}, \mathcal{G}\}$, we have $-\langle \mathcal{C}(-\nabla P) \rangle_{\mathcal{C}} = 0$ because the density is constant in each phase and only the pressure transport term contributes.

• The terms $\langle \mathcal{V} \rangle_{\mathcal{C}}$, $\langle \mathcal{S} \rangle_{\mathcal{C}}$ and $\langle \mathcal{F} \rangle_{\mathcal{C}}$ represent the scale-by-scale contribution of viscous diffusion, surface tension and forcing, respectively. Their general compact expression is given in Eq. (A9) of Appendix A. The expanded version of the viscous term has been derived by Lai et al. (2018). As with the pressure term, each of these terms can be decomposed into a contribution that looks similar to its constant-density counterpart, to which is added the correction, noted $\mathcal{C}$, which accounts for variations of density within the averaging volume. For the same reason as above, these correction terms vanish when $\mathcal{C} \in \{\mathcal{L}, \mathcal{G}\}$ because density is constant per phase. Note also that the term due to surface tension $\langle \mathcal{S} \rangle_{\mathcal{C}}$ contributes to the budget only if $\mathcal{C} = \mathcal{T}$.

It is worth stressing that the viscous term $\langle \mathcal{V} \rangle_{\mathcal{C}}$ can be decomposed into different contributions: a transport term in scale space $r$, a transport term in position space $\mathbf{X}$ and a contribution due to the kinetic energy dissipation rate (see for instance Galtier & Banerjee 2011; Lai et al. 2018; Hellinger et al. 2021b). Such a decomposition will not be carried out here, although it may provide many important insights into the role of viscosity in the scale-by-scale kinetic energy budget. This aspect could be addressed in future work.

In appendix A, we also address the asymptotic behaviour of all terms in the KHM equation when the separation $r$ tends to infinity. Under the condition of statistical homogeneity, it is proved that all terms of the KHM equation tend towards (four
times) their counterpart in the one-point scale-integrated kinetic energy budgets. For the flow setup studied here, they write:

$$d_t \langle k \rangle_{\mathrm{T}} = \langle F \rangle_{\mathrm{T}} - \langle \epsilon \rangle_{\mathrm{T}} + \langle S \rangle_{\mathrm{T}}, \qquad (2.13a)$$

$$d_t \langle k \rangle_{\mathrm{L}} = \langle F \rangle_{\mathrm{L}} - \langle \epsilon \rangle_{\mathrm{L}} + \langle T_{\nu} \rangle_{\mathrm{L}} + \langle T_p \rangle_{\mathrm{L}}, \qquad (2.13b)$$

$$d_t \langle k \rangle_{\mathrm{G}} = \langle F \rangle_{\mathrm{G}} - \langle \epsilon \rangle_{\mathrm{G}} + \langle T_{\nu} \rangle_{\mathrm{G}} + \langle T_p \rangle_{\mathrm{G}}, \qquad (2.13c)$$

The turbulent kinetic energy is given by $k = \frac{1}{2}\rho|\mathbf{u}|^2$. The terms $F$, $\epsilon$, $S$ represent the contribution due to forcing, kinetic energy dissipation and surface tension, respectively. The terms $T_\nu$ and $T_p$ correspond to the transport of kinetic energy due to viscous diffusion and pressure. In absence of forcing $F \equiv 0$, one recovers from Eqs. (2.13) the equations firstly derived by Dodd & Ferrante (2016). The recent work by Trefftz-Posada & Ferrante (2023) extended these equations to homogeneous shear two-phase flows. The explicit expression for $T_\nu$ and $T_p$ is given by Dodd & Ferrante (2016) and in Appendix A. It is important to notice that for the total fluctuating field, the pressure term is absent from the one-point kinetic energy budget Eq. (2.13)(a) although the baropycnal work is a priori non-zero in the KHM equation. This means that this term is likely to represent a transfer or conversion of energy between scales although its net contribution when integrated over all scales is zero. This was already anticipated by Aluie (2013) for the coarse-grained kinetic energy budget in compressible or variable-density flows.

Compared to the analysis of Lai et al. (2018) or Yao & Papadakis (2023), we consider here the KHM equation for the total velocity field without carrying any decomposition between its (Reynolds or Favre) mean and its (Reynolds or Favre) fluctuation. The reason is that (i) we want these equations to remain as compact as possible and (ii) the flow we will consider hereafter has zero mean. In flows with strong inhomogeneities, it may however be relevant to incorporate such a decomposition in order to reveal some mechanisms such as the transport and production of turbulent fluctuations by the mean flow.

The KHM equation for variable density multiphase flows Eq. (2.7) provides an explicit representation of the effects of unsteadiness, inertia, pressure, surface tension or viscous effect on the evolution of turbulent kinetic energy at different scales of the turbulent spectrum. Similar to the wavelet approach presented by Freund & Ferrante (2019) or the coarse-grained methodology followed by Pandey et al. (2020), this framework allows the evolution of turbulent kinetic energy to be represented in the compound position($\mathbf{X}$)/scale($\mathbf{r}$) space. Hence, it can easily be conditionally averaged in the same way as Freund & Ferrante (2019), i.e. by considering 3 different situations where the two points lie within either the liquid or within the gas phase, or where the two points separation crosses an interface.

Note that here the evolution for the total kinetic energy is recovered by applying the limit to infinite separations while for the wavelet approach or with the spectral analysis, the latter is obtained by integrating over all scales. Hence, the structure functions should rather be interpreted as a cumulative distribution of energy rather than an energy density distribution. Some authors proposed some definitions for the energy density using scale derivatives of the structure functions (Davidson & Pearson 2005; Danaila et al. 2012a; Hamba 2015, 2018; Arun et al. 2021). This could be addressed in a follow-up study.
## 2.2. Discussions

At this stage of the theoretical developments, it seems important to draw attention to a few points requiring careful consideration.

First-of-all, in several studies, the analysis of scale-by-scale kinetic energy budget is performed in spectral space (Crialesi-Esposito et al. 2022, 2023b; Ramirez et al. 2024; Cannon et al. 2024). It is worth stressing that the use of Fourier transforms may raise some issues. Firstly, spectra can be computed only over homogeneous directions with periodic boundary conditions. This is a very stringent restriction that makes spectra inapplicable to all flow situations, in particular in statistically inhomogeneous flows. In addition, as discussed by notably Lucci et al. (2010); Duret et al. (2012); Ramirez et al. (2024), two-phase flows reveal local discontinuities associated with the presence of the interface. This can result in oscillations in the kinetic energy spectra thereby corrupting the physical interpretation of the associated budget. Note though that (Ramirez et al. 2024) succeeded in distinguishing the spectral contribution due to discontinuities from the one due to the bulk phases allowing them to be interpreted separately. Another issue arises due to the spatial non-locality of Fourier transforms that precludes the definition of phase-conditioned spectral quantities. Ramirez et al. (2024) showed however that spectral approaches might remain valuable to explore the dynamics of the two-phase mixture, i.e. the case $C = T$, once the spectral signatures of discontinuities properly interpreted. If one aims at inferring the influence of the interface on the dynamics of each phase separately, analysis in physical space using either coarse-graining approaches, wavelet decompositions or point-splitting methods based on correlation or structure functions should be preferred.

Secondly, the derivation of the present form of the KHM equation requires the transport equation for both velocity **u** and momentum $\rho u$. The latter is obtained rigorously from the two-fluid formulation of the two-phase incompressible Navier-Stokes equation and is valid in the sense of generalized function. On the other hand, the equation for **u** derives from the additional application of the incompressibility condition. As stressed by Ramirez et al. (2024), this equation is mathematically ill-defined since one has to divide the surface tension term, in which appears a Dirac delta function, by the density, which is a Heaviside function. The KHM equation derived here as many other formulations of scale-by-scale energy budgets (for instance Pandey et al. 2020; Narula et al. 2025; Freund & Ferrante 2019) thus inherit this limitation. However, in most numerical solvers, whether they use front-capturing or diffuse-interface techniques, the discontinuities are regularized: in front-capturing methods at the mesh scale, and in diffuse-interface methods across the prescribed interface thickness. Hence, the one-fluid formulation of the NS equation is solved in a regularized sense. The KHM equation derived here as several other scale-by-scale formulations is thus valid in the same regularized sense. This limitation should be kept in mind when interpreting the results obtained from the present theoretical framework.

The third cautionary note concerns the somewhat arbitrary definition of the scale-by-scale kinetic energy $|\delta\mathbf{u}|_{\rho}^2$. As stressed by Aluie (2013); Narula et al. (2025), in the coarse-grained framework, there is no unique way to define the scale-by-scale kinetic energy in flows with significant fluctuations of density. The same applies to the point-splitting methods used here (see Appendix A.3). The choice made here is motivated by the fact that, at least, it recovers the correct one-point scale-integrated kinetic energy when taking the limit to infinite separations. Other definitions filling this condition exists (see Appendix A.3). In absence of consensus on the most appropriate definition,
it is perhaps wise to consider that the present definition is just one possible choice among others. At this stage, some comparative analyses using different definitions of the scale-by-scale kinetic energy should be carried out since this arbitrary choice may have consequences on the physical interpretation of the results.

Indeed, the definition proposed here yields the presence of a non-zero pressure term in the budget for the liquid-gas mixture ($C \equiv T$). In the lines of Aluie (2011, 2013), it is referred to as the baropycnal work and interpreted as being responsible for the transfer of kinetic energy among scales. It cancels out when integrating over all scales, and thus disappears from the one-point kinetic energy budget. This phenomenon must not be mistaken for the conversion of kinetic to internal energy by pressure-dilatation, which is a feature unique to compressible flows. As detailed in Appendix A.3, regardless of the definition, when a sort of density weighted (i.e. Favre) averaging is used to define the scale-by-scale kinetic energy, a consistent term emerges, correlating pressure gradient, velocity and density variations, underscoring the baropycnal work's role. In the coarse-grained framework, Narula et al. (2025) showed that this term may however contribute differently to the budget depending on the precise definition of the scale-by-scale kinetic energy. In the spectral energy budget derived by Ramirez et al. (2024), where the power of each physical process acting on the kinetic energy is expressed as the product of forces by velocity, there is no term similar to the baropycnal work. Therefore, one cannot rule out that the baropycnal term in the KHM equation arises as an artefact of the chosen definition of scale-by-scale kinetic energy, rather than reflecting an unambiguous physical process. This fundamental question calls for additional research.

# 3. Numerical simulations of multiphase turbulence

## 3.1. Numerical methods

The framework described above is used to analyse data from numerical simulation of turbulent multiphase flows issued from the code archer (Ménard et al. 2007; Vaudor et al. 2017). This code solves the one-fluid formulation of the NS equation Eq. (2.3) on a Cartesian staggered grid. The convective term in Eq. (2.3) is written in conservative form and solved using the improved Rudman (1998) technique (see Vaudor et al. 2017, for more details). The latter allows mass and momentum to be transported consistently enabling better accuracy and stability for flows with large liquid/gas density ratios. The viscous term in Eq. (2.3) is computed following the method presented by Sussman et al. (2007). The surface tension term is resolved using the Ghost-Fluid method (Fedkiw et al. 1999).

The interface is transported using a coupled level-set and volume-of-fluid (CLSWOF) solver, in which the level-set function accurately describes the geometric features of the interface (its normal and curvature) and the volume-of-fluid (VOF) function ensures mass conservation. For more information about the archer solver, the reader can refer to e.g. Ménard et al. (2007); Duret et al. (2012); Vaudor et al. (2017).

For time advancement, we use a slightly modified version of the fast three steps Runge-Kutta algorithm (abbreviated fastRK3) proposed by Aithal et al. (2023). More details on this algorithm are given in Appendix B.1.

## 3.2. Numerical configuration

The theoretical framework described in §2 is applied to numerical data of homogeneous isotropic turbulence that is maintained at steady state by adding a forcing term to the
NS equation. This configuration has now become a standard for studying multiphase turbulence (see e.g Duret et al. 2012; McCaslin & Desjardins 2014; Loisy & Naso 2017; Mukherjee et al. 2019; Thiesset et al. 2020; Boukharfane et al. 2021; Riviere et al. 2021; Cannon et al. 2024, among others). Here, we use the large-scale stochastic forcing of Eswaran & Pope (1988) (see also Cannon et al. 2024). The numerical implementation of this forcing procedure is described in Appendix B.2.

The initial condition consists in one spherical droplet placed at the centre of the domain, and 6 half-cylinders of same diameter centred on each face of the periodic box. This configuration is set from the very beginning of the simulation and not released at a given time. The initial velocity field corresponds to an ABC flow with $u' = 1$ that evolves to finally reach a steady state thanks to the forcing term. We have checked that the initial condition has no effect on the final statistical properties. This initial condition was chosen as it rapidly achieves the statistically steady state, hence saving some computational time. For each simulation, after a transient period of $4s \approx 8T_L$ where $T_L$ is the eddy turn-over time, a steady state is reached. The simulations are then run for 16 additional seconds (physical time), corresponding to roughly $32T_L$. Statistics are gathered throughout the statistically steady period. Statistical convergence is estimated using the 95% confidence level (assuming Gaussian distributions), and is found to be within few percent. In most figures below, the statistical uncertainty will be displayed as shaded region surrounding the averaged value.

The relevant non-dimensional parameters for our study are

the liquid volume fraction $\alpha := \langle \phi_L \rangle_T$, (3.1a)

the density ratio $R_{\rho} := \frac{\rho_L}{\rho_G}$, (3.1b)

the viscosity ratio $R_{\mu} := \frac{\mu_L}{\mu_G}$, (3.1c)

the Reynolds number of phase C $Re_C := \frac{\rho_C u'^2 L}{\mu_C}$, (3.1d)

the Weber number of phase C $We_C := \frac{\rho_C u'^2 L}{\sigma}$, (3.1e)

where $\phi_L$ is the liquid phase indicator function defined in (A.2). Throughout the study we set $L = 1$, $\rho_G = 1$, $u'^2 = \langle \rho | \mathbf{u}^2 \rangle / \langle \rho \rangle = 1$. In order to simplify the analysis, we decided to use a constant kinematic viscosity $\nu = \mu/\rho$ between the two phases. Hence, $R_{\mu}^{-1} = R_{\rho}$, so that only one Reynolds number needs to be defined $Re \equiv \nu^{-1}$ (recall that $u' = 1$ and $L = 1$). In addition, the effect of surface tension is studied as a function of an averaged Weber number which is defined by $We := \bar{\rho}u'^2L/\sigma$, where $\bar{\rho} = (\rho_L + \rho_G)/2$. Since $u' = 1$ and $L = 1$, the averaged Weber number reduces to $We \equiv \bar{\rho}/\sigma$.

In summary, our parameter space simplifies to 4 non-dimensional numbers $(\alpha, R_{\rho}, Re, We)$, which can be explored by varying the amount of liquid in the box, the liquid density $\rho_L$, the kinematic viscosity $\nu$ and the surface tension $\sigma$. In the present study, we have performed 14 different numerical simulations covering a range of $(\alpha, R_{\rho}, Re, We)$ as described in Table 1. This database is also supplemented by three additional simulations of single-phase turbulence covering the same range of Reynolds numbers. For the sake of clarity in Table 1, we have isolated five different groups from the 14 different simulations. The two first groups corresponds to variations of the density ratio $R_{\rho}$ for two different volume fractions: $\alpha = 25\%$
corresponding to what is referred to as the 'drop' regime, and $\alpha = 75\%$ corresponding to the 'bubble' regime. Note that the Weber number is fixed which means that the surface tension is adjusted for each $R_\rho$ so that $We = \bar{\rho}/\sigma = 25$. The third group corresponds to variations of the liquid volume fraction $\alpha$ between 12.5% and 87.5%. Here, the density ratio is set to $R_\rho = 25$. The case with $R_\rho = 1$ has already been extensively investigated (e.g. Mukherjee et al. 2019; Crialesi-Esposito et al. 2022; Cannon et al. 2024, among others) and is thus not reproduced here. The fourth and fifth groups correspond to variations of the Reynolds number and the Weber number, respectively. Here again, $R_\rho = 25$ while the liquid volume fraction is $\alpha = 25\%$. The last group of Table 1 corresponds to the single-phase flow simulations.

In Table 1, we also report the values for the total energy injection rate per unit mass, noted $\epsilon_f = \langle F \rangle_\mathrm{T}/\langle \rho \rangle_\mathrm{T}$, which by virtue of Eqs. (2.13) is also the total kinetic energy dissipation rate. The associated Taylor-based Reynolds number $R_\lambda = u'^2 \sqrt{5/(3\nu\epsilon_f)}$ and the Kolmogorov microscale $\eta = \nu^{3/4}/\epsilon_f^{1/4}$ are also given. The latter is given in terms of grid cell size $dx$ and is equal to roughly 2 for all cases considered here. Table 1 further reports the surface area of the liquid-gas interface $A_\Gamma$ and some typical scales noted $r_H, r_S$ and $r_c$ which will be defined and discussed later.

The reader is referred to appendix B.3 for some further discussions about the appropriateness of the numerical resolution. Due to resolution constraints, the Reynolds number achieved in the present numerical database remains moderate, preventing the development of an inertial range and even a restricted scaling range. As a result, the flow statistics reflect a superposition of large, intermediate, and small scales precluding exploring the scaling exponent for kinetic energy at intermediate scales. Nevertheless, the scaling for $\langle |\delta\mathbf{u}|_\rho^2 \rangle_\mathrm{C}$ at small scales remains amenable to exploration. This is discussed in more details in Appendix A.5 and in the results' section below. The description of the post-processing procedures for computing the different terms of the scale-by-scale energy budget is given in appendix B.4.

# 4. Results

We now present the application of the KHM equation to the numerical data. We sequentially focus on the effect of the liquid-gas density ratio (in subsections 4.1 and 4.2), the liquid volume fraction (subsection 4.3), the Reynolds number (subsection 4.4) and Weber number (subsection 4.5).

## 4.1. Effect of the liquid-gas density ratio in the drop regime

The scale distribution of kinetic energy $\langle |\delta\mathbf{u}|_\rho^2 \rangle_\mathrm{C}$ for different density ratio $R_\rho$ in the drop regime ($\alpha = 25\%$) is plotted in Fig. 1. The scale is normalized by the width of the computational domain $L$ while $\langle |\delta\mathbf{u}|_\rho^2 \rangle_\mathrm{C}$ is normalized by twice the scale-integrated turbulent kinetic energy $\langle |\mathbf{u}|_\rho^2 \rangle_\mathrm{C}$. For the total fluctuating field $\mathbf{C} \equiv \mathbf{T}$ (Fig. 1(a)), we observe that the small-scale content of kinetic energy is larger in multiphase turbulence than in single phase turbulence. We further note that increasing the density ratio yields an increase of the scale-by-scale kinetic energy at small scales. The activity of turbulence at small scales is linked to the activity of the dissipative scales. Increasing $R_\rho$, and/or moving from single-phase to multiphase flows, thus leads to an increase of the energy dissipation. This is confirmed by the values of $\epsilon_f$ in the first five lines of Table 1.

Looking at the curves for the liquid and gas conditional statistics (Figs. 1(b) and (c) for $\mathbf{C} \equiv \mathbf{L}$ and $\mathbf{C} \equiv \mathbf{G}$, respectively), we note that the scale distribution of kinetic
<figure><img src="image_1.png" /><figcaption>Figure 1: Effect of the density ratio $R_\rho$ on the scale-by-scale kinetic energy for $\alpha = 25\%$. The colours from blue to red correspond to $R_\rho = 1, 5, 25, 125$ while the black curve is for single-phase turbulence. (a) $C \equiv T$, (b) $C \equiv L$ and (c) $C \equiv G$. The dashed line represents the $r^2$ scaling.</figcaption></figure>

<figure><img src="image_2.png" /><figcaption>Figure 2: 2D slices of the vorticity magnitude (in colour, units of $u'/L$) together with the liquid-gas interface (black curves) for $\alpha = 25\%$ and increasing $R_\rho$ from (a) to (d).</figcaption></figure>
<table><caption>Table 1: List of simulation parameters. The number of simulation points per direction is $N$, hence the resolution $dx = L/N$. The Taylor-scale Reynolds number $R_λ = u'^2/\sqrt{5/(3\nu\epsilon_f)}$, where the kinetic energy injection rate per unit mass (= dissipation) $\epsilon_f = \langle F \rangle_\Gamma / \langle \rho \rangle_\Gamma$. $\epsilon_f$ is given in units of $u'^3/L (= 1)$. The Kolmogorov length scale is defined by $\eta = \nu^{3/4}/\epsilon_f^{1/4}$. The liquid-gas surface area $A_Γ$ is given in units of $L^2 (= 1)$. We also report the Kolmogorov-Hinze scale $r_H$, the scale $r_S$ and $r_c$ which are defined and discussed later. All are given in units of $L$.</caption><thead><tr><th>studied parameter</th><th>α [%]</th><th>R<sub>ρ</sub></th><th>ν [m<sup>2</sup>.s<sup>-1</sup>]</th><th>We</th><th>Re</th><th>N</th><th>ε<sub>f</sub></th><th>R<sub>λ</sub></th><th>η/dx</th><th>A<sub>Γ</sub></th><th>r<sub>H</sub></th><th>r<sub>S</sub></th><th>r<sub>c</sub></th></tr></thead><tbody><tr><td>liquid</td><td>25</td><td>1</td><td>1.4×10<sup>-3</sup></td><td>25</td><td>714</td><td>256</td><td>0.99</td><td>34.6</td><td>1.86</td><td>7.10</td><td>0.100</td><td>0.061</td><td>0.010</td></tr><tr><td>density, drop regime</td><td>25</td><td>5</td><td>1.4×10<sup>-3</sup></td><td>25</td><td>714</td><td>256</td><td>1.09</td><td>33.0</td><td>1.81</td><td>6.21</td><td>0.099</td><td>0.057</td><td>0.010</td></tr><tr><td>liquid</td><td>25</td><td>25</td><td>1.4×10<sup>-3</sup></td><td>25</td><td>714</td><td>256</td><td>1.15</td><td>32.2</td><td>1.79</td><td>5.67</td><td>0.096</td><td>0.051</td><td>0.011</td></tr><tr><td>volume fraction</td><td>25</td><td>125</td><td>1.4×10<sup>-3</sup></td><td>25</td><td>714</td><td>256</td><td>1.23</td><td>31.1</td><td>1.76</td><td>5.76</td><td>0.095</td><td>0.049</td><td>0.012</td></tr><tr><td>liquid</td><td>75</td><td>1</td><td>1.4×10<sup>-3</sup></td><td>25</td><td>714</td><td>256</td><td>0.99</td><td>34.6</td><td>1.86</td><td>7.10</td><td>0.100</td><td>0.062</td><td>0.010</td></tr><tr><td>volume fraction</td><td>75</td><td>5</td><td>1.4×10<sup>-3</sup></td><td>25</td><td>714</td><td>256</td><td>0.99</td><td>34.6</td><td>1.86</td><td>7.10</td><td>0.099</td><td>0.057</td><td>0.011</td></tr><tr><td>liquid</td><td>75</td><td>25</td><td>1.4×10<sup>-3</sup></td><td>25</td><td>714</td><td>256</td><td>0.98</td><td>34.9</td><td>1.86</td><td>6.11</td><td>0.098</td><td>0.056</td><td>0.012</td></tr><tr><td>volume fraction</td><td>75</td><td>25</td><td>1.4×10<sup>-3</sup></td><td>25</td><td>714</td><td>256</td><td>0.92</td><td>35.9</td><td>1.89</td><td>3.39</td><td>0.096</td><td>0.047</td><td>0.011</td></tr><tr><td>liquid</td><td>12.5</td><td>25</td><td>1.4×10<sup>-3</sup></td><td>25</td><td>714</td><td>256</td><td>1.20</td><td>31.5</td><td>1.77</td><td>7.80</td><td>0.096</td><td>0.054</td><td>0.012</td></tr><tr><td>volume fraction</td><td>50</td><td>25</td><td>1.4×10<sup>-3</sup></td><td>25</td><td>714</td><td>256</td><td>0.85</td><td>37.3</td><td>1.93</td><td>3.87</td><td>0.098</td><td>0.054</td><td>0.013</td></tr><tr><td>liquid</td><td>87.5</td><td>25</td><td>1.4×10<sup>-3</sup></td><td>25</td><td>714</td><td>256</td><td>0.85</td><td>37.3</td><td>1.93</td><td>3.87</td><td>0.098</td><td>0.054</td><td>0.013</td></tr><tr><td>volume fraction</td><td>25</td><td>25</td><td>3.8×10<sup>-3</sup></td><td>25</td><td>263</td><td>128</td><td>1.36</td><td>18.0</td><td>1.81</td><td>5.62</td><td>0.106</td><td>0.063</td><td>0.013</td></tr><tr><td>Reynolds number</td><td>25</td><td>25</td><td>5.3×10<sup>-4</sup></td><td>25</td><td>1887</td><td>512</td><td>0.99</td><td>57.3</td><td>1.81</td><td>5.55</td><td>0.089</td><td>0.042</td><td>0.109</td></tr><tr><td>liquid</td><td>25</td><td>25</td><td>1.4×10<sup>-3</sup></td><td>12.5</td><td>714</td><td>256</td><td>1.09</td><td>33.3</td><td>1.82</td><td>3.92</td><td>0.136</td><td>0.061</td><td>0.015</td></tr><tr><td>volume fraction</td><td>25</td><td>25</td><td>1.4×10<sup>-3</sup></td><td>50</td><td>714</td><td>256</td><td>1.12</td><td>32.6</td><td>1.80</td><td>8.12</td><td>0.070</td><td>0.042</td><td>0.088</td></tr><tr><td>liquid</td><td>-</td><td>-</td><td>3.8×10<sup>-3</sup></td><td>-</td><td>263</td><td>128</td><td>0.91</td><td>21.9</td><td>2.00</td><td>-</td><td>-</td><td>-</td><td>-</td></tr><tr><td>volume fraction</td><td>-</td><td>-</td><td>1.4×10<sup>-3</sup></td><td>-</td><td>714</td><td>256</td><td>0.66</td><td>42.0</td><td>2.04</td><td>-</td><td>-</td><td>-</td><td>-</td></tr><tr><td>liquid</td><td>-</td><td>-</td><td>5.3×10<sup>-4</sup></td><td>-</td><td>1887</td><td>512</td><td>0.60</td><td>71.3</td><td>2.02</td><td>-</td><td>-</td><td>-</td><td>-</td></tr></tbody></table>

energy in the liquid is mildly changing contrary to the gas phase where kinetic energy substantially increases with $R_ρ$ when the separation $r$ decreases. Hence, the increase of $\langle |\delta \mathbf{u}|_\rho^2 \rangle_\Gamma$ at small scales is mostly associated to an increase of $\langle |\delta \mathbf{u}|_\rho^2 \rangle_\mathcal{G}$.

Appendix A.5 discusses the small-scale behaviour of the two-point kinetic energy. It is shown that the mixed structure function $\langle |\delta \mathbf{u}|_\rho^2 \rangle_\mathcal{C}$ scales as $r^2$ and is proportional to the entrophy $\langle |\omega|^2 \rangle_\mathcal{C}$ averaged over $\mathcal{C} \in \{\mathbb{T}, \mathbb{L}, \mathcal{G}\}$. This quadratic evolution with respect to $r$ is shown in Fig. 1(a-c) as gray dashed lines. As $r$ decreases, $\langle |\delta \mathbf{u}|_\rho^2 \rangle_\mathcal{C}$ approaches the expected scaling, although some small deviations are observed, particularly in the gas phase (Fig. 1(c)). Previous studies (Estivalezes et al. 2022; Ling et al. 2019) have noted that accurately resolving the entrophy field in multiphase flow simulations is demanding. Therefore, while mass and kinetic energy are well resolved in our data (see Appendix B.3), the observed deviations from the $r^2$-scaling might reflect the need for an even finer mesh to fully capture the entrophy field. Consequently, in the remainder of this paper, the small-scale scaling of $\langle |\delta \mathbf{u}|_\rho^2 \rangle_\mathcal{C}$ will not be further addressed, and the vorticity field will be discussed only qualitatively.
Keeping this consideration in mind, our data indicate that the increase of the small-scale kinetic energy in the gas phase observed in Fig. 1(c) is related to an increase of the entrophy field in the gas phase. From $R_ρ = 1$ to $R_ρ = 125$, the averaged entrophy in the gas phase is doubled, as shown in Fig. 1(c). To confirm this, one can scrutinize the visualizations of the flow portrayed in Fig. 2. We observe an increased amplitude of the vorticity magnitude, hence the small scale activity, when one approaches the interface from the gas phase. This behaviour was already reported by Dodd & Ferrante (2016). The growing production of kinetic energy and vorticity is occurring in the turbulent boundary layers and wakes surrounding the liquid structures. As per Dodd & Ferrante (2016), this results from an increase of the liquid structures Stokes number which increases with $R_ρ$. The observed increased vorticity magnitude is further consistent with the theory developed by Terrington et al. (2022).

The different terms of the scale-by-scale kinetic energy budget are displayed in Fig. 3. Each contribution is presented in units of $\langle F \rangle_C$ which is the contribution of the forcing term in the one-point kinetic energy budget. By doing this, the scale distribution of forcing term to the scale-by-scale budget appears to be independent of the investigated physical parameters and of the averaging volume $C \in \{\mathbb{T}, \mathbb{L}, \mathbb{G}\}$, as shown in Fig. 3(a-c). The forcing term also appears to be the same as in single-phase flows. Consequently, all the other terms of the scale-by-scale budget can be studied in proportion of an invariant forcing contribution. This significantly eases the interpretation of the results.

Let us start by analysing the unconditionally averaged budget, i.e. $C \equiv \mathbb{T}$, portrayed in Fig. 3(a). We observe that the forcing term $\langle \mathcal{F} \rangle_{\mathbb{T}}$ contributes positively to the budget and acts at rather large-scales. At intermediate scales, the kinetic energy is transferred from large to small scales through the non-linear energy transport term $\langle \mathcal{T} \rangle_{\mathbb{T}}$ and the surface tension term $\langle \mathcal{S} \rangle_{\mathbb{T}}$. The kinetic energy budget is finally equilibrated by the viscous term $\langle \mathcal{V} \rangle_{\mathbb{T}}$ which is negative. Therefore, our analysis based on the KHM framework confirms previous conclusions using other approaches such as e.g. spectral energy budgets (see Crialesi-Esposito et al. 2022; Cannon et al. 2024, among others).

We further note that the non-linear transfer term $\langle \mathcal{T} \rangle_{\mathbb{T}}$ is drastically reduced compared to its single-phase counterpart. The kinetic energy exchange in multiphase turbulence thus appears to comply with the following scenario. (i) Part of the kinetic energy injected at large scales is transferred into kinetic energy at smaller scales by the classical non-linear transport process $\langle \mathcal{T} \rangle_{\mathbb{T}}$. This transfer is however less important than in single-phase flows because (ii) part of the kinetic energy injected at large scale is also transferred into surface energy and then released back into kinetic energy, but at smaller scales. Finally, (iii) the increase of the overall kinetic energy transfer is compensated by a shift of the viscous term $\langle \mathcal{V} \rangle_{\mathbb{T}}$ in the direction of smaller scales. The pressure term also contributes as a loss of energy at small scales as soon as $R_ρ \neq 1$. Its contribution is though rather marginal.

Although not shown in Fig. 3(a), summing up the contribution of the surface tension term and the non-linear transport term indicates that the overall kinetic energy transfer (i.e. $\langle \mathcal{T} \rangle_{\mathbb{T}} + \langle \mathcal{S} \rangle_{\mathbb{T}}$) is larger than the one pertaining to single-phase flows and peaks at smaller scales. It suggests that the overall transfer of kinetic energy in multiphase flows is quite likely to compare to the one observed in single-phase flows but at a larger Reynolds number. This explains that numerical simulations of multiphase turbulence require a much finer resolution compared to single-phase turbulence, even when $R_ρ = 1$.
<figure><img src="image_3.png" /><figcaption>Figure 3: Effect of the density ratio $R_\rho$ on the scale-by-scale kinetic energy budgets for $\alpha = 25\%$. The colours from blue to red correspond to $R_\rho = 1, 5, 25, 125$ as in the legend of Fig. 1. The black curves are for the single-phase case at same viscosity. The different lines correspond to
- $\langle T \rangle_c$, - $\langle P \rangle_c$, --- $\langle S \rangle_c$, --- $\langle V \rangle_c$, - $\langle F \rangle_c$, with (a) $C = T$, (c) $C = L$ and (d) $C = G$. The horizontal dashed line indicates $4\langle F \rangle_c$, the limit at large separations of the forcing term. Figure (b) represents the density correction terms $\langle C(a) \rangle_T$. Also represented in (a) are the scales $r_H$ with lines and $r_S$ with filled circles that will be described later on.</figcaption></figure>

When the density ratio $R_\rho$ increases, the contribution of the surface tension term $\langle S \rangle_T$ gets larger while the one due to the non-linear transport term $\langle T \rangle_T$ goes in opposite direction. We note also that $\langle S \rangle_T$ peaks at smaller scales compared to $\langle T \rangle_T$. The scale at which $\langle S \rangle_T$ is maximum is represented by the filled circle symbols in Fig. 3(a). More insights into the behaviour of this scale will be provided later in a dedicated section.

Consistently with our previous conclusions suggesting an increased activity of the smaller scales when $R_\rho$ becomes larger, the contribution of the viscous term starts
being visible at smaller scales when $R_\rho$ increases. Finally, the contribution due to pressure $\langle \mathcal{P} \rangle_\text{T} = -\langle \mathcal{C}(-\nabla \mathcal{P}) \rangle_\text{T}$ is rather small and appears only in the small-scale range where it contributes negatively to the budget. It is zero for $R_\rho = 1$, and increases in amplitude with $R_\rho$.

In Fig. 3(b), we represent the contribution of $-\langle \mathcal{C}(\mathbf{a}) \rangle_\text{T}$ for the viscous, forcing and surface tension terms. The one corresponding to the pressure term is already plotted and discussed in Fig. 3(a). As expected, $\langle \mathcal{C}(\mathbf{a}) \rangle_\text{T}$ is zero for $R_\rho = 1$ and its contribution increases when $R_\rho$ grows. We also have $\langle \mathcal{C}(\mathbf{a}) \rangle_\text{T} \to 0$ when $r \to \infty$ (see appendix A.4). The sign of the density correction for each term is the same as the term itself which means that variations of density tends to amplify the contribution of each process. We note that the density correction term contributes to about 1/3 of the surface tension at $R_\rho = 125$, and about 10% to the viscous term.

The different terms of the conditionally averaged KHM equation are presented in Fig. 3(c) and (d) for the liquid and gas phase, respectively. The first striking observation is for the gas phase (Fig. 3(d)), where we note that the pressure term $\langle \mathcal{P} \rangle_\text{G}$ is positive and thus represents a gain for the scale-by-scale kinetic energy budget. Recall that when conditionally averaged, only the pressure transport contributes to the budget while the term $\langle \mathcal{C} \rangle_\text{L,G}$ is zero. Hence the significant contribution attributed to pressure arises solely from the scale-by-scale pressure transport. Its contribution increases greatly with $R_\rho$ and is even exceeding the one of the forcing term $\langle \mathcal{F} \rangle_\text{T}$ at $R_\rho = 125$. This additional amount of kinetic energy injected into the gas phase yields an amplified contribution of the non-linear transport term $\langle \mathcal{T} \rangle_\text{G}$. When the density ratio $R_\rho = 125$, the latter is observed to be larger than the one pertaining to single-phase flows. For $R_\rho = 1$, the sum of the pressure and non-linear transport, i.e. $\langle \mathcal{T} \rangle_\text{G} + \langle \mathcal{P} \rangle_\text{G}$ is roughly similar to the non-linear transport term of single-phase turbulence. The effect of the viscous stress $\langle \mathcal{V} \rangle_\text{G}$ at each scale also becomes larger with $R_\rho$. Here again, this is explained by an increased activity of the small scales in the gas phase.

In the liquid-phase, our results show that the curves corresponding to $R_\rho = 1$ distinguishes substantially from the others. For instance, it appears that the pressure term $\langle \mathcal{P} \rangle_\text{L}$ (the pressure transport) for $R_\rho = 1$ is very different from the one at $R_\rho = 5$. It seems however to vary only marginally for $R_\rho \ge 25$. The same applies for the viscous term $\langle \mathcal{V} \rangle_\text{L}$ but not for the non-linear transfer term $\langle \mathcal{T} \rangle_\text{L}$ which continues to grow as $R_\rho$ increases. That means that for the liquid phase, the effect of the density contrast between the two phases is not trivial. Some processes seem to saturate while some others continues to evolve. Therefore, the physics observed for $R_\rho = 1$ which is a common simplification in the literature, cannot be readily extrapolated to flows with variations of density, even if the density contrast is moderate.

## 4.2. Effect of the liquid-gas density ratio in the bubble regime

The same analysis is performed for $\alpha = 75\%$ and $R_\rho = 1, 5, 25$. The scale distribution of kinetic energy $\langle |\delta \mathbf{u}|_\rho^2 \rangle_\text{C}/2\langle |\mathbf{u}|_\rho^2 \rangle_\text{C}$ is presented in Fig. 4(a), (b) and (c) for $\mathcal{C} = \text{T}, \text{L}$ and $\mathcal{G}$, respectively. Here, we note that $\langle |\delta \mathbf{u}|_\rho^2 \rangle_\text{T}/2\langle |\mathbf{u}|_\rho^2 \rangle_\text{C}$ is mostly unchanged when the density ratio $R_\rho$ is increased. Hence, the total dissipation rate $\varepsilon_f$ is found to be constant with respect to $R_\rho$ (see Table 1). At small scales, the scale-by-scale kinetic energy is slightly larger than that of single-phase turbulence. This is reflected by a larger dissipation rate $\varepsilon_f$ in multiphase turbulence compared to its single-phase counterpart. In the liquid phase, $\langle |\delta \mathbf{u}|_\rho^2 \rangle_\text{L}$ is almost independent of $R_\rho$ and is very comparable to single-phase turbulence. In the bubble regime, the effect of density
<figure><img src="image_4.png" /><figcaption>Figure 4: Effect of the density ratio $R_\rho$ on the scale-by-scale kinetic energy for $\alpha = 75\%$. The colours from blue to red correspond to $R_\rho = 1, 5, 25$ while the black curve is for single-phase turbulence. (a) $C \equiv T$, (b) $C \equiv L$ and (c) $C \equiv G$.</figcaption></figure>

<figure><img src="image_5.png" /><figcaption>Figure 5: 2D slices of the vorticity magnitude (in colour, units of $u'/L$) together with the liquid-gas interface (black curve) at $\alpha = 75\%$ and for increasing $R_\rho$ from (a) to (c).</figcaption></figure>

ratio is only perceptible in the gas phase. We note again a systematic shift of $\langle |\delta\mathbf{u}|_{\rho}^2 \rangle_G$ towards smaller scales when the density ratio is increased. Here again, this trend is associated with an increased intensity of the small-scales, which is reflected by a larger dissipation rate and a larger viscous transport in the gas phase.

This is indeed confirmed by looking at the visualizations of the vorticity magnitude presented in Fig. 5. Although less apparent than in Fig. 2 for the drop regime, we still observe larger amplitude of vorticity in the gas phase close to the interface. This is in agreement with the analytical considerations of Terrington et al. (2022).

The different contributions to the KHM equation are portrayed in Fig. 6. For the total fluctuating field (Fig. 6(a)), the effect of the density contrast between the two phases is perceptible on the surface tension term $\langle S \rangle_T$ which appears to decrease in amplitude and slides towards smaller scales when $R_\rho$ increases. These variations are quite limited though. We further confirm that the transfer associated to the non-linear transport $\langle T \rangle_T$ is much weaker than in single-phase turbulence. The overall transfer of kinetic energy $\langle T \rangle_T + \langle S \rangle_T$ is substantially larger than single-phase turbulence and is shifted towards smaller scales. This is compensated by a larger contribution of the viscous term in the small scale region. The scale-by-scale contributions of the 'density correction' terms $\langle C(\mathbf{a}) \rangle_T$ (Fig. 6(b)) are similar to what was observed before.
They have same sign as the term to which they apply and increase in amplitude with the density contrast $R_ρ$. Hence, these results do not vary substantially from those related to the drop regime.

Scrutinizing the conditionally averaged budget in the liquid or gas phase reveals that $R_ρ$ has a strong effect on the pressure and viscous term (see Fig. 6(d)). In the gas phase, pressure transport appear to act as an additional gain in the budget which is compensated by the viscous term. These processes dominate at $R_ρ = 25$ as they reach values roughly 5 times larger than the forcing term. In the liquid phase (see Fig. 6(c)), we note that the statistics are mostly unchanged between $R_ρ = 5$ and $R_ρ = 25$. For $R_ρ \ge 5$ and $\alpha = 75\%$, the mass of the liquid represent more than 90% of the total mass. Hence, in the bubble regime, one can readily assume that the dynamics of the liquid phase dominate the whole flow, and thus increasing further $ρ_L$ has only a limited impact. Only the case with $R_ρ = 1$ differ from the two other set of curves. The transfer of kinetic energy in the liquid phase due to the non-linear transport is systematically smaller than in single-phase turbulence. This effect might be the consequence of the smaller $R_λ$ in multiphase turbulence compared to its single-phase counterpart (see Table 1). The Reynolds number $R_λ$ does not vary much when $R_ρ$ varies.

### 4.3. Effect of the liquid volume fraction $\alpha$

The analysis is now repeated to emphasize the effect of the liquid volume fraction $\alpha$. We start by investigating the scale-by-scale kinetic energy which is plotted in Fig. 7(a-c) for $C \equiv T$, $L$ and $G$, respectively. We note that the liquid volume fraction has a non-monotonic effect on the total kinetic energy distribution across scales. Indeed, scrutinizing the small-scale region in Fig. 7(a) reveals that $\langle |\delta\mathbf{u}|^2_\rho \rangle_T$ first increase when $\alpha$ grows from 12.5% to 50% before decreasing again for $\alpha$ from 50% to 87.5%. The scale distribution of kinetic energy in multiphase turbulence tends towards the one of single-phase turbulence when $\alpha \to 1$, i.e. when the liquid drives the whole flow dynamics.

The kinetic energy scale distribution in the liquid (Fig. 7(b)) shows that the case with $\alpha = 12.5\%$ departs significantly from the different curves. At such a small liquid volume fraction, some confinement effects are likely to be at play which precludes turbulent eddies of size larger than the typical size of liquid structures to develop, thereby impeding the turbulent activity in the liquid phase. The scale distribution of kinetic energy in the liquid for the other values of $\alpha$ is roughly invariant. The scale-by-scale turbulent kinetic energy conditioned by the gas phase is presented in Fig. 7(c). The contribution of the small-scales to the kinetic energy is observed to increase with $\alpha$. Here again, the scenario proposed above stating that the liquid phase having much larger inertia yields an increased production of vorticity and thus small-scale kinetic energy in the gas phase is likely to hold true. Hence, the more liquid, the more kinetic energy production at small scales in the gas phase. Despite this increased activity of the small scales in the gas phase, table 1 indicates that the total kinetic dissipation rate $\epsilon_f$ is maximum for $\alpha = 50\%$. This is obviously due to the decreasing proportion in mass and volume of the gas phase when $\alpha$ increases.

The visualizations presented in Fig. 8 do not really show any substantial difference of the vorticity magnitude when $\alpha$ changes. However, when $\alpha \to 1$, the volume of the gas phase which lies close to the interface increases relatively to the total gas volume. Hence, the volume where the vorticity magnitude is larger increases in proportion. This is likely to explain why $\langle |\delta\mathbf{u}|^2_\rho \rangle_G$ increases in the small-scale range when the liquid volume fraction increases.
<figure><img src="image_6.png" /><figcaption>Figure 6: Effect of the density ratio $R_\rho$ on the scale-by-scale kinetic energy budgets for $\alpha = 75\%$. The colours from blue to red correspond to $R_\rho = 1, 5, 25, 125$ as in the legend of Fig. 1. The black curves are for the single-phase case at same viscosity. The different lines correspond to
- $\langle T \rangle_c$, - $\langle P \rangle_c$, --- $\langle S \rangle_c$, --- $\langle V \rangle_c$, - $\langle F \rangle_c$, with (a) $C = T$, (c) $C = L$ and (d) $C = G$. The horizontal dashed line indicates $4\langle F \rangle_c$, the limit at large separations of the forcing term. Figure (b) represents the density correction terms $\langle C(\alpha) \rangle_T$.</figcaption></figure>

The different terms of the KHM equation for different volume fractions are presented in Fig. 9. Figure 9(a) indicates that the unconditionally averaged terms of the KHM equation corresponding to a given $\alpha$ (e.g. 12.5 and 25%) differ greatly from those corresponding to an equivalent $1 - \alpha$ (i.e. 87.5 and 75%). For instance, the surface tension term $\langle S \rangle_T$ for $\alpha = 12.5\%$ appears to be much larger than the one for $\alpha = 87.5\%$, even though the surface area $A_T$ of the interface is similar (see Table 1). The transfer due to the non-linear transport $\langle T \rangle_T$ has also a non-monotonic evolution with respect to the liquid volume fraction. Nevertheless, it seems to tend towards its single phase counterpart when $\alpha \to 1$. The overall transfer of kinetic energy which is due to both
<figure><img src="image_7.png" /><figcaption>Figure 7: Effect of the liquid volume fraction $\alpha$ on the scale-by-scale kinetic energy. The colours from blue to red correspond to $\alpha = 12.5, 25, 50, 75, 87.5\%$.</figcaption></figure>

(a) $C = T$, (b) $C = L$ and (c) $C = G$.

<figure><img src="image_8.png" /><figcaption>Figure 8: 2D slices of the vorticity magnitude (in colour, units of $u'/L$) together with the liquid-gas interface (black curves) for increasing $\alpha$ from (a) to (e).</figcaption></figure>

the non-linear transport and surface tension term is compensated by the viscous term $\langle \mathcal{V} \rangle_T$. The latter is shifted towards larger scales when $\alpha \ge 50\%$. Here again, the term $\langle \mathcal{V} \rangle_T$ tends towards its single-phase counterpart when $\alpha \to 1$.

The contribution of the density correction terms $\langle C(\mathbf{a}) \rangle_T$ is presented in Fig. 9(b). Except for the forcing term, we note that the correction terms are maximum when $\alpha = 50\%$. This suggests that $\langle C(\mathbf{a}) \rangle_T$ should be related to some extent to the surface area of the liquid-gas interface which is minimum for $\alpha = 12.5\%$ and 87.5% and maximum at $\alpha = 50\%$. The contribution of $\langle C(S) \rangle_T$ also peaks at very different scales.
Its contribution is perceptible at small scales when $\alpha$ is small and systematically
moves towards the large-scales when $\alpha$ becomes larger. For the viscous term, the
correction appears to peak around the same scale, independently of $\alpha$.

The terms of the KHM equation conditionally averaged in the liquid and gas phase
are presented in Figs. 9(c) and (d), respectively. For the liquid phase, we note that
when $\alpha \to 1$, the contribution of the pressure transport term $\langle \mathcal{P} \rangle_{\mathrm{L}}$ tends to zero and
the other terms (the viscous $\langle \mathcal{V} \rangle_{\mathrm{L}}$ and non-linear transport term $\langle \mathcal{T} \rangle_{\mathrm{L}}$) tends towards
those obtained in single-phase flows. A careful examination of Fig. 9(d) reveals that
the same holds true for the gas phase when $\alpha \to 0$. For small values of $\alpha$ (i.e. the drop
regime), the pressure term in the liquid phase $\langle \mathcal{P} \rangle_{\mathrm{L}}$ contributes negatively at large
scales which means that it represents a loss for the liquid kinetic energy budget at
those scales. On the other hand, the budget in the gas phase (Fig. 9(d)) indicates that
the pressure term $\langle \mathcal{P} \rangle_{\mathrm{G}}$ contributes positively and very substantially when $\alpha \to 1$, i.e.
in the bubble regime. Note the amplitude of the pressure term $\langle \mathcal{P} \rangle_{\mathrm{G}}$ for $\alpha = 87.5\%$
which is roughly 6 times larger than the forcing term $\langle \mathcal{F} \rangle_{\mathrm{G}}$. This means that when
the liquid-volume fraction is large, most of the kinetic energy production in the gas
phase comes from pressure transport. This process is compensated by an increased
contribution from the viscous term $\langle \mathcal{V} \rangle_{\mathrm{G}}$.

In summary, varying the liquid volume fraction reveals that the physics of the drop
regime ($\alpha < 50\%$) differs significantly from the one of the bubble regime ($\alpha > 50\%$).
This stresses again that the results obtained by setting $R_o = 1$ can very hardly be
representative of what occurs in 'real' multiphase flows characterized by fluctuations,
even moderate, of the fluid density. Our results also reveal that the budget in the
most abundant phase, when the minority phase represents less than say 5%, could
possibly correspond to the one observed in single phase turbulence.

4.4. *Effect of the Reynolds number Re*

The effect of Reynolds number *Re* is studied independently of the other parameters
by varying the kinematic viscosity (see Table 1). The kinetic energy distribution
among the different scales is plotted in Fig. 10. Decreasing the value of the kinematic
viscosity shifts the viscous cut-off towards the small scales. The distribution of kinetic
energy is thus spread over a wider range of scales and the small scales are more
energetic when the Reynolds number increases. Here again, we further note that
for $\mathbb{C} = \mathbb{T}$ (Fig. 10(a)), the small-scale kinetic energy of multiphase turbulence is
larger than its single-phase counterpart at the same Reynolds number. As with Fig.
1(a), the distribution conditioned in the liquid is very comparable to single-phase
turbulence, only the distribution in the gas phase deviates significantly. Hence, the
larger kinetic energy observed for $\mathbb{C} = \mathbb{T}$ at small scales is due to a larger kinetic
energy in the gas phase.

The visualizations of the vorticity field presented in Fig. 11 indicate again that
these small scales are produced in the gas phase in close vicinity of the interface.
When the Reynolds number is increased, there are more vortical structures, they
are of smaller size and their magnitude increases. This is reflected in the increase of
$\langle |\delta\mathbf{u}|_{\rho}^2 \rangle_{\mathrm{T}}$ and $\langle |\delta\mathbf{u}|_{\rho}^2 \rangle_{\mathrm{G}}$ at small-scales. Despite this, the amount of surface area of the
two-fluid interface appears to be independent of the Reynolds number (see Table 1).

The terms of KHM equation are displayed in Fig. 12 for (a,b) the total fluctuating
field, (c) and (d) the liquid and gas phase, respectively. For the total fluctuating field
(Fig. 12(a)), decreasing the viscosity increases the separation between large and small
scales and thus the intermediate scales develop. The viscous term is shifted towards
smaller scales and the terms acting at the intermediate range, i.e. $\langle \mathcal{T} \rangle_{\mathrm{T}}$ and $\langle \mathcal{S} \rangle_{\mathrm{T}}$
<figure><img src="image_9.png" /><figcaption>Figure 9: Effect of the liquid volume fraction $\alpha$ on the scale-by-scale kinetic energy budgets. The colours from blue to red correspond to $\alpha = 12.5, 25, 50, 75, 87.5\%$ as in the legend of Fig. 7, while the black curves correspond to the single-phase case. The different lines correspond to $\langle T \rangle_c, \langle P \rangle_c, \langle S \rangle_c, \langle V \rangle_c, \langle F \rangle_c$, with (a) $C = T$, (c) $C = L$ and (d) $C = G$. The horizontal dashed line indicates $4\langle F \rangle_c$, the limit at large separations of the forcing term. Figure (b) represents correction terms $C(a)$</figcaption></figure>

widen and increase in amplitude. This is qualitatively equivalent to what occurs in single-phase turbulence. The same observation stands true for the additional term due to density variations $\langle C(a) \rangle_T$ presented in Fig. 12(b). The overall transfer of kinetic energy which is composed of both the non-linear transport and the surface tension term, exceeds the one of single-phase turbulence at the same Reynolds number. As observed before, this is compensated by the viscous term being active at smaller scales. The pressure term does not contribute much to the budget of $\langle |\delta u|^2 \rangle_T$.

Looking at the terms of the KHM equation in liquid and gas phase (Fig. 12(c,d)), it is noted that $\langle T \rangle_{L,G}$ increases in amplitude and peaks at smaller scales when the
<figure><img src="image_10.png" /><figcaption>Figure 10: Effect of the Reynolds number $Re \equiv 1/\nu$ on the scale-by-scale kinetic energy. The colours from blue to red correspond to $Re = 263, 714, 1887$. The grey curves correspond to the single-phase (SP) data at same Reynolds number. (a) $C = T$, (b) $C = L$ and (c) $C = G$.</figcaption></figure>

<figure><img src="image_11.png" /><figcaption>Figure 11: 2D slices of the vorticity magnitude (in colour, units of $u'/L$) together with the liquid-gas interface (black curves) for increasing $Re$ from (a) to (c).</figcaption></figure>

Reynolds increases. Its contribution is quite smaller to that of single-phase turbulence. Surprisingly, the pressure term appears to be independent of the Reynolds number. This applies to both the liquid and gas phase. Pressure transport is again found to contribute positively in the gas phase. For the liquid phase, it is a gain at small scales and a loss at large scales.

### 4.5. *Effect of the Weber number We*

We now investigate the effect of Weber number $We$ with $R_\rho = 25$, $\alpha = 25\%$ and $Re = 714$. The scale-by-scale kinetic energy $\langle |\delta\mathbf{u}|^2 \rangle_C$ is displayed in Fig. 13(a), (b) and (c) for $C \equiv T$, $L$, and $G$, respectively.

As an overall observation, the scale-by-scale kinetic energy appears to vary only marginally when the surface tension is varied. One observes though some very slight modifications in the intermediate range of scales for the unconditionally averaged kinetic energy $\langle |\delta\mathbf{u}|^2 \rangle_T$ (Fig. 13(a)). The small-scales remain more or less unchanged consistently with the values of $\epsilon_f$ reported in Table 1. The same small variations in the intermediate range of scales are also visible in the liquid phase in Fig. 13(b). In the gas phase (Fig. 13)(c), the distribution of kinetic energy remain independent of the Weber number in the large and intermediate range of scales. Some marginal
<figure><img src="image_12.png" /><figcaption>Figure 12: Effect of the Reynolds number $Re \equiv 1/\nu$ on the scale-by-scale kinetic energy budgets. The colours from blue to red correspond to $Re = 263, 714, 1887$ as in the legend of Fig. 10. The grey curves correspond to the single-phase data at same Reynolds number. The grey curves correspond to the single-phase data at same Reynolds number. The different lines correspond to $\langle T \rangle_c$, $\langle P \rangle_c$, $\langle S \rangle_c$, $\langle V \rangle_c$, $\langle F \rangle_c$, with (a) $C = T$, (b) $C = L$ and (c) $C = G$. The horizontal dashed line indicates $4\langle F \rangle_c$, the limit at large separations of the forcing term. Figure (b) represents the density correction terms $C(a)$.</figcaption></figure>

differences appear when one travels towards smaller scales where $\langle |\delta\mathbf{u}|^2 \rangle_G$ appears to increase very slightly with $We$. Looking at $\langle |\delta\mathbf{u}|^2 \rangle_T$ in Fig. 13(a), the small-scale kinetic energy of multiphase turbulence is again observed to be larger than the one pertaining to single-phase turbulence. This results from an increased small-scale activity in the gas phase (see Fig. 13(c)).

The direct visualizations of the flow presented in Fig. 14 indicate that the magnitude of the vorticity field does not vary substantially when $We$ is varied. The increase of the Weber yields more pronounced corrugations of the interface and hence a larger
<figure><img src="image_13.png" /><figcaption>Figure 13: Effect of the Weber number $We \equiv \bar{\rho}/\sigma$ on the scale-by-scale kinetic energy. The colours from blue to red correspond to $We = 12.5, 25, 50$, while the black curves correspond to the single-phase case. (a) $C = T$, (b) $C = L$ and (c) $C = G$.</figcaption></figure>

<figure><img src="image_14.png" /><figcaption>Figure 14: 2D slices of the vorticity magnitude (in colour, units of $u'/L$) together with the liquid-gas interface (black curves) for increasing $We$ from (a) to (c).</figcaption></figure>

surface area (see also Table 1). The theory of Terrington et al. (2022) indicates that one particular term of the vorticity production is proportional to the surface tension $\sigma$ and the gradient along the interface of the mean curvature. Hence, when the Weber number increases, the decrease of surface tension is likely compensated by larger gradients of curvature (more corrugations), leading to roughly the same vorticity magnitude, and hence a comparable small-scale kinetic energy distribution in the gas phase.

The effect of surface tension on the scale-by-scale kinetic energy budget is presented in Fig. 15. Scrutinizing the unconditionally averaged terms (Fig. 15(a)), we note that the surface tension term $\langle S \rangle_T$ drops when surface tension decreases. It is thus plausible that the scale-by-scale contribution of the surface tension term will tend to zero at infinite Weber number. The Weber number further appears to have an influence on the pressure term $\langle P \rangle_T$, whose contribution changes sign when the Weber changes. This behaviour is rather surprising. The two other terms of Eq. (2.2) due to the non-linear transport $\langle T \rangle_T$ and the viscous term $\langle V \rangle_T$ are only marginally influenced by variations of the surface tension. However, they are substantially different from the ones obtained in single-phase flows. For instance, the viscous term starts to contribute
at much smaller scales to compensate an increased overall transfer of kinetic energy, which is again much larger in multiphase than in single-phase flows.

The density correction terms are portrayed in Fig. 15(b). As observed previously, we note that $-\langle C(S) \rangle_T$, $-\langle C(V) \rangle_T$ and $-\langle C(F) \rangle_T$ have same sign as $\langle S \rangle_T$, $\langle V \rangle_T$ and $\langle F \rangle_T$, respectively. The contribution of these terms seem to decrease when the Weber number increases.

The terms of Eq. (2.2) conditionally averaged budget in the liquid and gas phase are displayed in Fig. 15(c) and (d), respectively. Noticeable is the evolution of the pressure transport term in the liquid phase $\langle P \rangle_L$ when the surface tension is varied. The latter is positive at small scales (a gain) and negative (a loss) at large scales. The scale at which the change of sign occurs decreases with surface tension. The effect of surface tension is also clearly visible on the pressure term in the gas phase, $\langle P \rangle_G$, which appears to increase when the Weber number increases. Its contribution to the scale-by-scale budget is positive (a gain) and is compensated mainly by an increased contribution of the viscous term $\langle V \rangle_G$.

## 4.6. Phenomenological predictions for the pivotal scale

So far, we have shown that surface tension yields an additional transfer term of kinetic energy between scales. Of particular interest is the scale noted $r_S$ at which the term $\langle S \rangle_T$ is maximum. Indeed, for scales larger than $r_S$, the interface pumps kinetic energy, stores it as surface energy before releasing it back into kinetic energy but at scales smaller than $r_S$.

As already mentioned in the introduction, (Crialesi-Esposito et al. 2023b) interpreted $r_S$ as the pivotal scale between, at large-scales, a range dominated by breakup, and at small scales, a range dominated by coalescence. They thus suggested that $r_S$ should coincide with the Kolmogorov-Hinze scale, noted $r_H$, which is known as the largest fluid parcel which would resist breakup.

Going back to the original idea of Hinze (1955), Mukherjee et al. (2019) extracted the Kolmogorov-Hinze scale from a local, i.e. scale dependent Weber number (see also Perlekar et al. 2012). We proceed similarly here, and the local Weber number is defined from the scale-by-scale kinetic energy, viz.

$$ We(r) := \frac{1}{2} \frac{\bar{\rho} \langle |\delta \mathbf{u}|_{\rho}^2 \rangle_{\mathrm{T}} r}{\langle \rho \rangle_{\mathrm{T}} \sigma} \quad (4.1) $$

Since $\langle |\delta \mathbf{u}|_{\rho}^2 \rangle_T / \langle \rho \rangle_T \to 2u'^2$ when $r \to L$, this definition of the local Weber number is such that $We(r) \to We$ as $r \to L$. The Kolmogorov-Hinze scale is then simply given as $We(r = r_H) = 1$.

Note that the classical definition for the Kolmogorov-Hinze scale can be recovered if the Reynolds number is sufficiently large for an inertial range to exist and for $r_H$ to fall in such an inertial range. In this situation, Kolmogorov scaling might apply for the Favre-averaged structure functions, viz.

$$ \frac{\langle |\delta \mathbf{u}|_{\rho}^2 \rangle_{\mathrm{T}}}{\langle \rho \rangle_{\mathrm{T}}} = C_K (\epsilon_f r)^{2/3}, \quad (4.2) $$

where $C_K$ is known as the Kolmogorov constant. This yields the following relation for Kolmogorov-Hinze scale:

$$ r_H = \left( \frac{2}{C_K} \right)^{3/5} \left( \frac{\sigma}{\bar{\rho}} \right)^{3/5} \varepsilon_f^{-2/5}, \quad (4.3) $$
<figure><img src="image_15.png" /><figcaption>Figure 15: Effect of the Weber number $We = \bar{\rho}/\sigma$ on the scale-by-scale kinetic energy budgets. The colours from blue to red correspond to $We = 12.5, 25, 50$ as in the legend of Fig. 13, while the black curves correspond to the single-phase case. The different lines correspond to $\langle T \rangle_c$, $\langle P \rangle_c$, $\langle S \rangle_c$, $\langle V \rangle_c$, $\langle F \rangle_c$, with (a) $C = T$, (c) $C = L$ and (d) $C = G$. The horizontal dashed line indicates $4\langle F \rangle_c$, the limit at large separations of the forcing term. Figure (b) represents the density correction terms $C(a)$</figcaption></figure>

which is very similar to the classical definition of the Kolmogorov-Hinze scale (see e.g. Hinze 1955; Perlekar et al. 2012; Crialesi-Esposito et al. 2022, 2023b; Cannon et al. 2024). Defining $r_H$ from Eq. (4.2) instead of Eq. (4.3) applies to cases where the Reynolds number is not sufficiently large for an inertial range scaling to be invoked. Note that Hinze (1955) introduced an *ad hoc* critical Weber number, $We_c$, which initially was aimed at accounting for the variations of viscosity. This can readily be done here as well by defining $r_H$ as $We(r=r_H) = We_c$. Note also that Eq. (4.3) differs from the classical definition of the Kolmogorov-Hinze scale through the appearance of the mean density $\bar{\rho}$ instead of the density of the most abundant fluid.
One can thus think of Eq. (4.3) as a likely more general expression for $r_H$ that applies to the dense regime where the notion of carrier versus dispersed phase is irrelevant.

The values for $r_S$ and $r_H$ are reported in Table 1. They are also represented in Figs. 3(a), 6(a), 9(a), 12(a), 15(a) as filled circles and vertical lines, respectively. Comparing the values of $r_S$ and $r_H$ reveals that the Kolmogorov-Hinze scale $r_H$ does not systematically scale with $r_S$. This is particularly visible when the density ratio between the two phases is varied. In this situation, $r_H$ is roughly constant while $r_S$ decreases when $R_\rho$ is varied from 1 to 125. It is also visible when the Weber number is varied. In this situation, both $r_H$ and $r_S$ decrease when We increases. However, $r_H$ does not vary in proportion of $r_S$.

The discrepancy between $r_H$ and $r_S$ was already noted by (Crialesi-Esposito et al. 2023b) at $R_\rho = 1$ but for varying viscosity ratio $R_\mu$. This difference could possibly be taken into account by using an ad hoc critical Weber number function of $R_\rho$. One could otherwise seek for another characteristic scale as done below.

Using some different forms of the scale-by-scale kinetic energy budget equations in gravity driven bubbly flows, Pandey et al. (2020, 2023); Ramirez et al. (2024) have shown that the surface tension term has qualitatively the same scale distribution as the one reported here or in Crialesi-Esposito et al. (2022, 2023b); Cannon et al. (2024). However, the regime of bubbly flows investigated by Pandey et al. (2020, 2023); Ramirez et al. (2024) does not reveal any breakup nor coalescence, and thus does not fall into the Kolmogorov-Hinze framework. This means that breakup/coalescence may not be required to explain the behaviour of the surface tension term.

Here, we extend somewhat the phenomenology proposed by Crialesi-Esposito et al. (2023b) and replace breakup/coalescence events by the more general mechanisms of interface deformation/restoration. We hypothesize that scales larger than $r_S$ contribute to the interface deformation, thereby loosing kinetic energy which is converted into surface energy. The mechanism of interface deformation can, but not necessarily includes breakup. On the other hand, for scales smaller than $r_S$ surface energy is released back into kinetic energy through interface restoration. By restoration, we mean that any corrugations of the interface smaller than $r_S$ are flattened through the action of surface tension. Coalescence can also contribute to the disappearance of scales smaller than $r_S$, thereby releasing kinetic energy at those scales. Breakup can also contribute to restoration. This is true for instance during the Rayleigh-Plateau instability where the deformed ligament have larger surface area than the formed droplets after breakup. For simplicity, the range of scales $r > r_S$ will be now referred to as the deformation range, and the range of scales $r < r_S$, the restoration range.

In recent work (Thiesset et al. 2021; Gauding et al. 2022; Thiesset & Federrath 2023), we introduced a general morphological analysis of corrugated interfaces. It was shown in particular in Gauding et al. (2022); Thiesset & Federrath (2023) that when observed at a sufficiently small scale, an interface appears flat. If observed at larger scales, the interface is seen to be curved though smooth, and can eventually be observed to be fractal at even larger scales. Building on the exact analytical results of Kirste & Porod (1962); Frisch & Stillinger (1963) (see also Gauding et al. 2022; Thiesset & Federrath 2023) a typical scale below which an interface is observed to be flat can be obtained. The latter writes:

$$r_c = \sqrt{\frac{8}{3\langle H^2 \rangle_s - \langle G \rangle_s}}, \quad (4.4)$$
where *H* and *G* represent mean curvature and Gaussian curvature, respectively, while $\langle \bullet \rangle_s$ stands for a surface area weighted average. For iso-surfaces formed by passive scalar Gauding et al. (2022) proved that $r_c$ scales with the Batchelor length-scale.

To say it differently, an observer zooming on the interface down to scale $r_c$ will see a flat surface. The curvature of the interface is visible when it is observed with a field of view larger than $r_c$. In other words, the smallest interface corrugation is related to $r_c$. This scale is thus likely to be a good candidate to represent $r_S$, the pivotal scale between the deformation and the restoration range. The values for $r_c$ computed from Eq. (4.4) are reported in Table 1. Unfortunately, although the above reasoning is appealing on the paper, we find that the variations of $r_c$ do not coincide with the variations $r_S$. The scale $r_S$ is systematically smaller and does not vary in proportion of $r_c$. Nevertheless, it is found that $r_c$ coincides rather satisfactorily with $r_H$ so that one can conclude that the Kolmogorov-Hinze framework provides a rather suitable prediction for the smallest corrugation scales of the interface.

# 5. Conclusions

In the present study, a KHM equation is derived and used to explore the scale-by-scale kinetic energy exchanges in multiphase turbulence. We consider both classical spatial averages and phase conditional averages. This allows the kinetic energy budget to be inferred in each phase separately. This framework is applied to numerical data of forced homogeneous isotropic turbulence simulated using the **archer** code. The database does not restrict to cases with $\rho_L = \rho_G$ and considers different values for the liquid/gas density ratio. The liquid volume fraction, Weber and Reynolds numbers were also varied. The statistics of multiphase turbulence are systematically compared to those of single-phase turbulence at the same Reynolds number.

From this methodology yield several outcomes. The most important are listed below.

Firstly, we confirm that the interface modifies the overall energy transfer between scales. The term associated with the surface tension is found to act similarly to an additional transfer term, although it should be more accurately interpreted as the energy conversion between kinetic and surface energy. The physical picture is then that part of the kinetic energy injected into the system is pumped from the large scales to deform the interface. Surface energy stored by the interface is then released back into kinetic energy but at smaller scales. In the meantime, kinetic energy is also transferred to the smaller scales by the classical non-linear transport, although at smaller rates when compared to its single-phase counterpart. Summing these two contributions, the overall rate of kinetic energy transfer/conversion is larger than in single-phase flows and is compensated by a larger contribution of the viscous term at small scales. In this respect, the visualizations of the vorticity magnitude have revealed that the presence of the interface creates more intense small-scale structures in the close vicinity of the interface in the gas phase.

Secondly, conditioning the KHM equation in either the liquid and gas phase have revealed the role played by pressure transport. The latter was found to act as a gain (loss) of scale-by-scale kinetic energy for the lighter (heavier) fluid. In the gas phase, the contribution due to pressure can also severely exceed the one of the forcing mechanism which is used to maintain the flow at statistically steady state. This additional pressure power is compensated by the viscous term. For the liquid-gas mixture also i.e. $C \equiv T$, a non-zero contribution arises from the pressure term. We refer to this term as the baropycnal work, following Aluie (2011, 2013). It is not
clear however at this stage whether it is merely an artefact of the Favre averaging or
whether it has a genuine physical significance.

Thirdly, we have shown that the physics of bubbles are very different to that of
drops. This conclusion arises as soon as one considers different densities between the
liquid and gas phases. Hence, the case with R<sub>ρ</sub> = 1 is a very special case that can
hardly be extrapolated to situations with some, even moderate, density jumps.

Fourthly, building on previous work, we hypothesize the existence of a pivotal scale, noted r<sub>S</sub>, between surface deformation (kinetic energy pumping) at large scales and surface restoration (kinetic energy release) at small scales. These two processes are believed to be slightly more general than those invoked by Crialesi-Esposito et al. (2023b) as they may apply to situations with no breakup nor coalescence, e.g. gravity driven bubbly flows. We have sought a suitable prediction for this scale. It appears that the Kolmogorov-Hinze scale r<sub>H</sub> may not scale as r<sub>S</sub>. Using a morphological analysis of corrugated interfaces, we extract a new characteristic scale, noted r<sub>c</sub>. The latter is argued to be the smallest scale below which the interface appears flat to an observer. Unfortunately, r<sub>S</sub> and r<sub>c</sub> do not scale with each other. However, the Kolmogorov-Hinze scale r<sub>H</sub> provides satisfactory prediction for r<sub>c</sub>.

The present work may also inspire future dedicated studies. For instance, in the
present database r<sub>H</sub>, r<sub>S</sub> and r<sub>c</sub> do not vary much, except maybe when the Weber
number is varied. A wider range of parameters should thus be considered to provide
further evidence on how these different scales evolve. At this stage, our hypothesis
stating that r<sub>S</sub> separates the deformation range and restoration range still remains
plausible, although no definite proofs are provided yet. Another open question is that
of the asymptotic behaviour of the different terms of the KHM equation in the limit
of large Reynolds and Weber numbers. Dedicated studies are required to fully capture
the behaviour that emerges at higher Reynolds regimes.

The discrepancy between $r_c$ and $r_S$ recalls that the scale-by-scale contribution of the surface tension term is obviously not only related to the geometrical features of the interface. It depends also on the coupling between the dynamical field (the velocity field), the interface geometry and its kinematics. Consequently, more insights could be provided by looking at the alignment between the interface normal, the velocity vector and/or the strain eigenvectors. Further, because the $\langle S \rangle_{\text{T}}$ term is assumed to depend on the creation or destruction of surface area at a given scale, one should probably look at some scale-by-scale (e.g. coarse grained) versions of the transport equation for the surface density.

Our data have also indicated an increased activity of the small scales associated
with an increase of R<sub>ρ</sub> or α. This was reflected in particular in larger vorticity
amplitudes in the gas phase close to the interface. More insight could be provided by
looking at the vorticity production mechanisms in the line of Terrington et al. (2022),
and references therein. In Appendix A.5, we provide first insights into the relation
between the KHM equation and the transport equation for the entrophy. This could
be continued in future analysis in order to shed light onto the vorticity dynamics in
multiphase systems.

Finally, we have heretofore investigated the KHM equation for the total, liquid,
and gas phase. The case referred to as the mixed case ($C = M$) was not considered.
This could be done in follow-up study with the aim of shedding light into the energy
exchanges between the two phases at a given scale r.

**Acknowledgements.** The present work would not have been feasible without the contribution of the members of the **archer** team, viz. J-C. Brändle de Motta, B. Duret, T. Ménard, A. Poux who are gratefully acknowledged. We are very grateful to P. Bousquet-Melou from CRIANN for
helping us optimising our GPU-processing code and R. Zamansky for providing us a first version of the forcing routine. This work has benefited from discussions with C. Federrath, R. Zamansky, A. Burlot, V. Boniou, J-L. Pierson, M. Crialesi-Esposito. Calculations were performed using the computing resources of CRIANN (Normandy, France), under the project 2018002.

**Funding.** This research received funding from the french research national agency (ANR) under the contract ANR-22-CE51-0005.

**Declaration of interests.** The authors report no conflict of interest.

**Data availability statement.** The data and the code can be made available upon reasonable request.

**Author ORCIDs.** F. Thiesset, https://orcid.org/0000-0001-7303-5106

**Author contributions.** Both authors contributed equally to all aspects of the present work.

# Appendix A. Detailed derivation of the scale-by-scale budgets

## A.1. Unaveraged equation

The mathematical machinery for deriving the different terms of the Kármán-Howarth-Monin for incompressible two-phase flows is same as already presented elsewhere by e.g. Galtier & Banerjee (2011); Lai et al. (2018); Hellinger et al. (2021b), initially developed in the context of either compressible flows and/or variable density miscible mixtures. The reader is referred to these papers for more details on the algebra. Here, we followed the exact same method as Hellinger et al. (2021b) but apply it to the one-fluid formulation to the two-phase incompressible NS-equation.

We start by the convective term, noted $\mathcal{T}$, which writes:

$$ \begin{aligned} \mathcal{T} &:= -(\delta \mathbf{u}) \cdot (\delta(\nabla \cdot \rho \mathbf{u} \mathbf{u})) - (\delta(\rho \mathbf{u})) \cdot (\delta((\mathbf{u} \cdot \nabla)\mathbf{u})) \\ &= -(\delta \mathbf{u}) \cdot (\delta(\mathbf{u} \cdot (\nabla \rho \mathbf{u}))) - (\delta(\rho \mathbf{u})) \cdot (\delta((\mathbf{u} \cdot \nabla)\mathbf{u})), \end{aligned} \quad (A 1) $$

where we have used the incompressibility condition $\nabla \cdot \mathbf{u} = 0$. It is then useful to define the derivatives w.r.t the separation $\mathbf{r}$ and the mid-point $\mathbf{X} := (\mathbf{x}^+ + \mathbf{x}^-)/2$ (see e.g. Hill 2002; Danaila et al. 2012b).

$$ \nabla_{\mathbf{r}} = \frac{1}{2} (\nabla^{+} - \nabla^{-}), \quad (A 2a) $$

$$ \nabla_{\mathbf{X}} = (\nabla^{+} + \nabla^{-}). \quad (A 2b) $$

For any quantity $\bullet$, we also have $\nabla^{+} \bullet^{-} = \nabla^{-} \bullet^{+} = 0$. After some manipulations, Eq. (A 1) can be rewritten as (see e.g. Hellinger et al. 2021b)

$$ \mathcal{T} = -\nabla_{\mathbf{X}} \cdot (\bar{\delta}\mathbf{u}) |\delta\mathbf{u}|_{\rho}^2 - \nabla_{\mathbf{r}} \cdot (\delta\mathbf{u}) |\delta\mathbf{u}|_{\rho}^2, \quad (A 3) $$

where $\bar{\delta} \bullet := (\bullet^{+} + \bullet^{-})/2$ is the arithmetic mean of any quantity $\bullet$ between point $\mathbf{x}^{+}$ and $\mathbf{x}^{-}$. Equation (A 3) reveals that the two-point convective term writes as the sum of two divergence terms, with respect to the geometrical space $\mathbf{X}$ and scale space $\mathbf{r}$, respectively.

For the other terms, we use the ruse of Hellinger et al. (2021b) who noticed that for any quantity $\mathbf{a}$, each term of the KHM equation has a generic form $\mathcal{D}(\mathbf{a})$ in the form

$$ \begin{aligned} \mathcal{D}(\mathbf{a}) &:= (\delta\mathbf{u}) \cdot (\delta\mathbf{a}) + (\delta(\rho\mathbf{u})) \cdot (\delta\mathbf{v}\mathbf{a}) \\ &= 2(\delta\mathbf{u}) \cdot (\delta\mathbf{a}) - \mathcal{C}(\mathbf{a}), \end{aligned} \quad (A 4) $$
where

$$ 
\begin{aligned}
\mathcal{C}(\mathbf{a}) &:= (\delta\mathbf{u}) \cdot (\delta\mathbf{a}) - (\delta(\rho\mathbf{u})) \cdot (\delta v\mathbf{a}) \\
&= (\rho^+v^- - 1)\mathbf{u}^+ \cdot \mathbf{a}^- + (\rho^-v^+ - 1)\mathbf{u}^- \cdot \mathbf{a}^+,
\end{aligned}
\quad (A5) $$

is a correction term which accounts explicitly for the variation of density. It is indeed easy to show that in constant density flows, we have $\rho^+v^- = \rho^-v^+ = 1$ and hence $\mathcal{C}(\mathbf{a}) = 0$.

For the term due to the pressure gradient, which is noted $\mathcal{P}$, we set $\mathbf{a} = -\nabla P$ in eq. (A 4), and obtain

$$ 
\begin{aligned}
\mathcal{P} &:= (\delta\mathbf{u}) \cdot (\delta(-\nabla P)) + (\delta(\rho\mathbf{u})) \cdot (\delta(-v\nabla P)) \\
&= 2(\delta\mathbf{u}) \cdot (\delta(-\nabla P)) - \mathcal{C}(-\nabla P).
\end{aligned}
\quad (A6) $$

The first term on RHS of Eq. (A 6) can be written as

$$ -2(\delta\mathbf{u}) \cdot (\delta(\nabla P)) = -2\nabla_{\mathbf{X}} \cdot (\delta\mathbf{u})(\delta P), \quad (A7) $$

where again we have used $\nabla \cdot \mathbf{u} = 0$. The term on RHS of Eq. (A 7) is generally referred to as the pressure transport. Using Eq. (A 7), the pressure term (Eq. A 6) can finally be recast in the form

$$ \mathcal{P} = -2\nabla_{\mathbf{X}} \cdot (\delta\mathbf{u})(\delta P) - \mathcal{C}(-\nabla P). \quad (A8) $$

Proceeding similarly for the surface tension term noted $\mathcal{S}$, the viscous term noted $\mathcal{V}$, and the forcing term $\mathcal{F}$, yields

$$ \mathcal{V} := 2(\delta\mathbf{u}) \cdot (\delta(\nabla \cdot \mathbf{t})) - \mathcal{C}(\nabla \cdot \mathbf{t}) \quad (A9a) $$

$$ \mathcal{S} := 2(\delta\mathbf{u}) \cdot (\delta\mathcal{S}) - \mathcal{C}(\mathcal{S}) \quad (A9b) $$

$$ \mathcal{F} := 2(\delta\mathbf{u}) \cdot (\delta\mathcal{F}) - \mathcal{C}(\mathcal{F}). \quad (A9c) $$

These terms can hardly be recast in a more compact or illustrative notation and are thus kept in this form. Note however that invoking additional hypothesis could lead to further simplifications. For instance, if a linear forcing is used, i.e. $\mathbf{F} := A\rho\mathbf{u}$ (Lundgren 2003; Rosales & Meneveau 2005), with a forcing amplitude $A$ that does not depend on space, one ends up with

$$ \mathcal{F} = 2A|\delta\mathbf{u}|_{\rho}^2 \quad (A10) $$

Lumping together the general formulation of the different terms leads to

$$ 
\begin{aligned}
\mathcal{D}_t &= \mathcal{T} + \mathcal{P} + \mathcal{V} + \mathcal{S} + \mathcal{F} \\
\partial_t |\delta \mathbf{u}|_{\rho}^2 &= -\nabla_{\mathbf{X}} \cdot (\bar{\delta} \mathbf{u}) |\delta \mathbf{u}|_{\rho}^2 - \nabla_{\mathbf{r}} \cdot (\delta \mathbf{u}) |\delta \mathbf{u}|_{\rho}^2 \\
&\quad -2 \nabla_{\mathbf{X}} \cdot (\delta \mathbf{u}) (\delta P) - \mathcal{C}(-\nabla P) \\
&\quad +2(\delta \mathbf{u}) \cdot (\delta(\nabla \cdot \mathbf{t})) - \mathcal{C}(\nabla \cdot \mathbf{t}) \\
&\quad +2(\delta \mathbf{u}) \cdot (\delta \mathcal{S}) - \mathcal{C}(\mathcal{S}) \\
&\quad +2(\delta \mathbf{u}) \cdot (\delta \mathcal{F}) - \mathcal{C}(\mathcal{F})
\end{aligned}
\quad (A11) $$

Equation (A 11) is the general unaveraged KHM equation applying to incompressible two-phase flows. In this formulation, Eq. (A 11) reveals that the effect of density contrasts are accounted for through the different correction terms $\mathcal{C}(\mathbf{a})$. These terms will be zero when the two-phases are assumed to have same density, which is a reasonable approximation in case of emulsions (e.g. Crialesi-Esposito et al. 2022). The final key distinction between Eq. (A 11) and its single-phase variable-density
analogue (see e.g. Galtier & Banerjee 2011; Lai et al. 2018; Hellinger et al. 2021b) lies in an additional term arising from surface tension, inherited from the surface tension contribution in the one-fluid formulation of the two-phase incompressible Navier–Stokes equations.

## A.2. Averaged equations

Eq. (A 11) should be supplemented by some averaging operations that depend on the flow situations. Here, the flow under consideration is stationary and homogeneous, thereby allowing to employ both temporal (over time $t$) and spatial (over geometrical positions in flow $\mathbf{X}$) averages. We will also employ angular averages over all orientations of the separation vector $\mathbf{r}$. These averages will be denoted by $\langle \bullet \rangle_{\mathbf{T}}$. In this situation, Eq. (A 11) simplifies to:

$$ \underbrace{\langle \partial_t |\delta \mathbf{u}|_{\rho}^2 \rangle_{\mathbf{T}}}_{\text{Time deriv } \mathcal{D}_t} = \underbrace{-\langle \nabla_r \cdot (\delta \mathbf{u}) |\delta \mathbf{u}|_{\rho}^2 \rangle_{\mathbf{T}}}_{\text{Transport } \mathcal{T}} \underbrace{-\langle \mathcal{C}(-\nabla P) \rangle_{\mathbf{T}}}_{\text{Pressure } \mathcal{P}} + 2\langle (\delta \mathbf{u}) \cdot (\delta(\nabla \cdot 2\mu\mathbf{S})) \rangle_{\mathbf{T}} - \langle \mathcal{C}(\nabla \cdot 2\mu\mathbf{S}) \rangle_{\mathbf{T}} \\ \underbrace{+2\langle (\delta \mathbf{u}) \cdot (\delta \mathbf{S}) \rangle_{\mathbf{T}} - \langle \mathcal{C}(\mathbf{S}) \rangle_{\mathbf{T}}}_{\text{Viscous } \mathcal{V}} + 2\langle (\delta \mathbf{u}) \cdot (\delta \mathbf{F}) \rangle_{\mathbf{T}} - \langle \mathcal{C}(\mathbf{F}) \rangle_{\mathbf{T}}. \qquad (\text{A 12}) \\ \underbrace{+2\langle (\delta \mathbf{u}) \cdot (\delta \mathbf{S}) \rangle_{\mathbf{T}}}_{\text{Surface tension } \mathcal{S}} + \underbrace{+2\langle (\delta \mathbf{u}) \cdot (\delta \mathbf{F}) \rangle_{\mathbf{T}} - \langle \mathcal{C}(\mathbf{F}) \rangle_{\mathbf{T}}}_{\text{Forcing } \mathcal{F}} $$

For the sake of completeness, we have kept the time derivative term in Eq. (A 12). It will however be zero in the sequel because the flow is at statistically steady state and the averaging volume has fixed, non-moving boundaries. The spatially averaged $\mathbf{X}$-divergence terms is not included in Eq. (A 12) as it vanishes due to periodicity (Hill 2002). The first term on RHS of Eq. (A 12) is a transfer term and writes as the divergence in scale-space of the flux $(\delta \mathbf{u})|\delta \mathbf{u}|_{\rho}^2$. As for its single-phase counterpart, this term represents the cascade process of turbulent kinetic energy across the different scales. The other terms represent the scale-by-scale contribution of pressure, viscous stress, surface tension and forcing.

In two-phase flows, it is also convenient to define conditional averages where statistics are gathered only within a given phase (Dodd & Ferrante 2016; Rosti et al. 2019; Crialesi-Esposito et al. 2022; Trefftz-Posada & Ferrante 2023). When applied to two-point statistics, one has to consider three different situations (Yao et al. 2020; Yao & Papadakis 2023):

* the two points $\mathbf{x}^+$ and $\mathbf{x}^-$ lie in the liquid phase, such conditional averages will be denoted with the subscript $\mathbb{L}$ for “liquid”.

* the two points $\mathbf{x}^+$ and $\mathbf{x}^-$ lie in the gas phase, such conditional averages will be denoted with the subscript $\mathbb{G}$ for “gas”.

* the two points $\mathbf{x}^+$ and $\mathbf{x}^-$ lie in different phases, such conditional averages will be denoted with the subscript $\mathbb{M}$ for “mixed”.

For computing such conditional averages (Yao et al. 2020; Gauding et al. 2021; Yao & Papadakis 2023) one needs to define a phase indicator function which for the liquid phase can be defined by

$$ \phi_L(\mathbf{x}, t) := \begin{cases} 1 & \text{if } \mathbf{x} \text{ lies in the liquid phase} \\ 0 & \text{elsewhere,} \end{cases} \qquad (\text{A 13}) $$
while for the gas phase,

$$ \phi_G(\mathbf{x}, t) := \begin{cases} 1 & \text{if } \mathbf{x} \text{ lies in the gas phase} \\ 0 & \text{elsewhere.} \end{cases} \qquad (A14) $$

Then, for any two-point quantity $\delta\bullet$ one may define the conditional averages as (Yao et al. 2020; Gauding et al. 2021; Yao & Papadakis 2023)

$$ \langle \delta \bullet \rangle_{\mathrm{T}} = \frac{\langle \phi_{L}^{+} \phi_{L}^{-} \delta \bullet \rangle_{\mathrm{T}}}{\gamma_{L}} \qquad (A15a) $$

$$ \langle \delta \bullet \rangle_{\mathrm{G}} = \frac{\langle \phi_{G}^{+} \phi_{G}^{-} \delta \bullet \rangle_{\mathrm{T}}}{\gamma_{G}} \qquad (A15b) $$

$$ \langle \delta \bullet \rangle_{\mathrm{M}} = \frac{\langle \phi_{L}^{+} \phi_{G}^{-} \delta \bullet \rangle_{\mathrm{T}}}{\gamma_{M}}, \qquad (A15c) $$

where,

$$ \gamma_L := \langle \phi_L^+ \phi_L^- \rangle_T \qquad (A16a) $$

$$ \gamma_G := \langle \phi_G^+ \phi_G^- \rangle_T \qquad (A16b) $$

$$ \gamma_M := \langle \phi_L^+ \phi_G^- \rangle_T. \qquad (A16c) $$

We thus have $\gamma_L + \gamma_G + 2\gamma_M = 1$, and for any two-point statistics $\delta\bullet$ (Yao & Papadakis 2023)

$$ \langle \delta \bullet \rangle_{\mathrm{T}} = \gamma_L \langle \delta \bullet \rangle_{\mathrm{L}} + \gamma_G \langle \delta \bullet \rangle_{\mathrm{G}} + 2\gamma_M \langle \delta \bullet \rangle_{\mathrm{M}}. \qquad (A17) $$

Therefore, the conditionally averaged KHM equation is simply obtained by pre-multiplying each term by a correlation of the phase indicator function before averaging over all space, time, and orientations. Any two-point statistics averaged over the whole domain can be retrieved from the sum of the conditional ones weighted by their corresponding phase indicator correlation function $\gamma$ (Yao & Papadakis 2023). These correlations are linked together by (see e.g. Torquato 2002, p 27-28)

$$ \gamma_G = 1 - 2\langle\phi_L\rangle_T + \gamma_L \qquad (A18a) $$

$$ \gamma_L = 1 - 2\langle\phi_G\rangle_T + \gamma_G \qquad (A18b) $$

$$ \gamma_M = \langle \phi_L \rangle_T - \gamma_L = \langle \phi_G \rangle_T - \gamma_G \qquad (A18c) $$

Although rather far from the scope of the present work, it is worth noting the correlation functions of the phase indicator field contains interesting information about the morphology of the liquid and gas phase interface. They can notably be used to infer the interface surface area, its mean and Gaussian curvatures, together with some fractal characteristics, if any (Gauding et al. 2022; Thiesset & Federrath 2023; Thiesset et al. 2020, 2021).

Applying conditional averages has some consequences in the final form of the KHM equation. For instance, the $X$-divergence terms do not vanish a priori, even when averaged over a periodic domain and/or homogeneous flow. The time derivative term should also be retained in the final formulation of the conditionally averaged KHM equation since the conditional averaging volume has non-fixed boundaries. Indeed, let us define the set $C(\mathbf{r}) := \{\mathbf{X} : \phi_C^+ \phi_C^- = 1\}$, where $C \in \{L, G\}$ and $C' \in \{L, G\}$. The conditional volume average of the time derivative term can then be rewritten as

$$ \langle \partial_t | \delta \mathbf{u} | \rho \rangle_{\mathrm{C}} = \frac{1}{\gamma_C} \int_{\mathrm{T}} \phi_C^+ \phi_C^- \partial_t | \delta \mathbf{u} |_{\rho}^2 d\mathbf{X} = \frac{1}{\gamma_C} \int_{\mathrm{C}} \partial_t | \delta \mathbf{u} |_{\rho}^2 d\mathbf{X}. $$
Then, by virtue of the Reynolds transport theorem,

$$ \langle \partial_t | \delta \mathbf{u} |_{\rho}^2 \rangle_{\mathbb{C}} = \frac{1}{\gamma_{\mathbb{C}}} \frac{d}{dt} \int_{\mathbb{C}} | \delta \mathbf{u} |_{\rho}^2 d\mathbf{X} - \frac{1}{\gamma_{\mathbb{C}}} \int_{\partial \mathbb{C}} | \delta \mathbf{u} |_{\rho}^2 \mathbf{u}_b \cdot \mathbf{n} \, dS, \quad (A 19) $$

where $\partial\mathbb{C}$ is the boundary of the control volume $\mathbb{C}$, $\mathbf{n}$ is the outwardly oriented normal vector to $\partial\mathbb{C}$ and $\mathbf{u}_b$ the velocity of the control volume boundary in the laboratory frame of reference. For immiscible fluids, the velocity of the control surface, $\mathbf{u}_b$, is the fluid velocity, $\mathbf{u}$, at the control surface $\partial\mathbb{C}$. While the leftmost term on RHS of Eq. (A 19) is zero when the flow is at statistically steady state, the rightmost term cannot however be dropped as it represents the surface averaged flux of $|\delta\mathbf{u}|_{\rho}^2$ across the control surface. Consequently, even in statistically steady flows, the conditionally averaged time derivative of the quantity $|\delta\mathbf{u}|_{\rho}^2$ term cannot be neglected. Since the rightmost term of Eq. (A 19) writes as a convective flux across the averaging volume boundary, it is worth being added to the convective term. This is what has been done in Eq. (2.7) to form the “transport” term $\langle \mathcal{T} \rangle_{\mathbb{C}}$.

Since the conditional averages are performed in the bulk of each phase which does not comprise the interface, the contribution of the surface tension term in the conditionally averaged two-point kinetic energy budget is zero. The same applied to the one-point kinetic energy budget (Dodd & Ferrante 2016).

In summary, the general conditionally averaged KHM equation reads

$$ \underbrace{\langle \partial_t | \delta \mathbf{u} |_{\rho}^2 \rangle_{\mathbb{C}}}_{\text{Time deriv. } \mathcal{D}_t} = \underbrace{-\langle \nabla_r \cdot (\delta \mathbf{u}) | \delta \mathbf{u} |_{\rho}^2 \rangle_{\mathbb{C}} - \langle \nabla_X \cdot (\bar{\delta} \mathbf{u}) | \delta \mathbf{u} |_{\rho}^2 \rangle_{\mathbb{C}}}_{\text{Transport } \mathcal{T}} \\ -2\langle \nabla_X \cdot (\delta \mathbf{u})(\delta P) \rangle_{\mathbb{C}} - \langle C(-\nabla P) \rangle_{\mathbb{C}} \\ \underbrace{+2\langle (\delta \mathbf{u}) \cdot (\delta (\nabla \cdot 2\mu \mathbf{S})) \rangle_{\mathbb{C}} - \langle C(\nabla \cdot 2\mu \mathbf{S}) \rangle_{\mathbb{C}}}_{\text{Pressure } \mathcal{P}} \\ +2\langle (\delta \mathbf{u}) \cdot (\delta \mathbf{F}) \rangle_{\mathbb{C}} - \langle C(\mathbf{F}) \rangle_{\mathbb{C}}, \qquad (A 20) \\ \underbrace{\text{Viscous } \mathcal{V}}_{\text{Forcing } \mathcal{F}} \\ \underbrace{\text{Forcing } \mathcal{F}}_{\text{Forcing } \mathcal{F}} $$

where the subscript $\mathbb{C} \in \{\mathbb{L}, \mathbb{G}, \mathbb{M}\}$.

Since the density is constant in each phase separately (and equals either $\rho_L$ or $\rho_G$), all variable-density correction terms $C(\mathbf{a})$ should be zero when conditionally averaged within either the liquid or gas phase. Eq. (A 20) can then be simplified in the special case $\mathbb{C} \in \{\mathbb{L}, \mathbb{G}\}$, viz.

$$ \underbrace{\langle \partial_t | \delta \mathbf{u} |_{\rho}^2 \rangle_{\mathbb{C}}}_{\text{Time deriv. } \mathcal{D}_t} = \underbrace{-\langle \nabla_r \cdot (\delta \mathbf{u}) | \delta \mathbf{u} |_{\rho}^2 \rangle_{\mathbb{C}} - \langle \nabla_X \cdot (\bar{\delta} \mathbf{u}) | \delta \mathbf{u} |_{\rho}^2 \rangle_{\mathbb{C}}}_{\text{Transport } \mathcal{T}} \underbrace{-2\langle \nabla_X \cdot (\delta \mathbf{u})(\delta P) \rangle_{\mathbb{C}}}_{\text{Pressure } \mathcal{P}} \\ +2\langle (\delta \mathbf{u}) \cdot (\delta (\nabla \cdot 2\mu \mathbf{S})) \rangle_{\mathbb{C}} \underbrace{+2\langle (\delta \mathbf{u}) \cdot (\delta \mathbf{F}) \rangle_{\mathbb{C}}}_{\text{Viscous } \mathcal{V}} \underbrace{\text{Forcing } \mathcal{F}}_{\text{Forcing } \mathcal{F}}. \quad (A 21) $$

Equation (A 21) applies true only if $\mathbb{C} \in \{\mathbb{L}, \mathbb{G}\}$. It is also valid for $\mathbb{C} = \mathbb{M}$ only if $\rho_L = \rho_G$. In Eq. (A 21), since the density is constant in each phase, it can be dropped from the average, and one recovers the classical single-phase constant-density KHM equation.

### A.3. Different definitions of the two-point kinetic energy

In variable density flows, the definition of the two-point kinetic energy is not unique (Aluie 2013; Narula et al. 2025). In the present work, we have used Eq. (2.1) following
Galtier & Banerjee (2011); Hellinger et al. (2021b); Lai et al. (2018). Different definitions can however be found in the literature. For instance, Hellinger et al. (2021a) define the kinetic energy as $|\delta w|^2$ where $w = \rho^{1/2}u$ as in Kida & Orszag (1990). Hellinger et al. (2021a) found that the KHM equation for this definition of the kinetic energy is virtually the same as Eq. (A 12). There appears a term analogous to the C-term for the pressure which now writes

$$C_{\sqrt{\rho}}(\nabla P) = ((\rho^+ v^-)^{1/2} - 1)u^+ \cdot (\nabla P)^- + ((\rho^- v^+)^{1/2} - 1)u^- \cdot (\nabla P)^+. \quad (A 22)$$

The latter can also be interpreted as the contribution of the baropycnal work. Similar density variation correction terms also arise from the viscous and forcing terms of the Navier-Stokes equation (Hellinger et al. 2021a). More generally, one could use any $0 \le \beta \le 1$ such that $\delta(\rho^\beta u) \cdot \delta(\rho^{1-\beta} u)$ is a possible definition for the scale-by-scale turbulent kinetic energy (Aluie 2013). Based on this definition for the scale-by-scale kinetic energy, the same kind of terms due to density variations are likely to operate in the KHM equation.

Another choice is made by Ferrand et al. (2020) and Brahami (2020) who define the two-point kinetic energy as $\bar{\delta}\rho|\delta u|^2$. In this case, Brahami (2020) derived an additional term due to density variations in the KHM equation (see their Eq. 2.39), which can also be associated to the effect of the baropycnal work. Some additional contributions due to density variations emerge from the viscous (and the forcing) term of the Navier-Stokes equations.

Because in multiphase flows, the density is piecewise constant, one could also have used the single-phase definition for the structure function, i.e. $\langle|\delta u|^2\rangle_C$, i.e. without any sort of density weighted increments. The latter could therefore be interpreted as a Reynolds-averaged version for the velocity fluctuations. This is indeed a possible definition. Although it should better be referred to as the “agitation” rather than the kinetic energy, it may solve the artefact associated with the Favre average (see §2.2). Note that this definition would lead to the exact same KHM equation as Eq. (2.7) for the case $C \in \{\mathbb{L}, \mathbb{G}\}$ since the density is constant per phase and could thus be dropped from the structure function.

In summary, a consensus on the optimal definition of two-point kinetic energy in variable density flows remains elusive. Narula et al. (2025) explored various definitions for scale-by-scale turbulent kinetic energy within the coarse-grained framework. Similar research is required to determine the most suitable definition when a point-splitting method is used.

### A.4. Limit at large separations

It is worth noting that for homogeneous flows, the scale-by-scale kinetic energy (Lai et al. 2018),

$$\langle|\delta u|_{\rho}^{2}\rangle_{\mathbb{T}} = 2\langle\rho|u|^{2}\rangle_{\mathbb{T}} - 4\langle(\bar{\delta}\rho)u^{+} \cdot u^{-}\rangle_{\mathbb{T}} = 4\langle k\rangle_{\mathbb{T}} - 4\langle(\bar{\delta}\rho)u^{+} \cdot u^{-}\rangle_{\mathbb{T}}, \quad (A 23)$$

where the turbulent kinetic energy is defined as $k = \frac{1}{2}\rho|u|^2$. Since, in the limit $|r| \to \infty$, the correlation $\langle(\bar{\delta}\rho)u^+ \cdot u^-\rangle_T \to 0$, the time derivative term in Eq. (A 12) has the following limit

$$\lim_{r \to \infty} \langle \partial_t | \delta \mathbf{u} | \rho \rangle_{\mathbb{T}} = 4 \frac{d}{dt} \langle k \rangle_{\mathbb{T}}. \quad (A 24)$$

Proceeding similarly for the other terms of the KHM equation, it can be shown that, in the limit of large separations, each term of the total KHM equation Eq. (A 12)
tends towards four times its counterpart in the one-point total kinetic energy budget. It is also easy to show that the variable-density corrections $C(a) \to 0$ when $r \to \infty$. As an example, the two-point viscous term can be proven to asymptote the kinetic energy dissipation rate (Hellinger et al. 2021b), viz.

$$ \begin{aligned} & \lim_{r \to \infty} +2\langle(\delta\mathbf{u}) \cdot (\delta(\nabla \cdot 2\mu\mathbf{S}))\rangle_{\mathrm{T}} - \langle C(\nabla \cdot 2\mu\mathbf{S})\rangle_{\mathrm{T}} \\ & = 4\langle \mathbf{u} \cdot (\nabla \cdot 2\mu\mathbf{S}) \rangle_{\mathrm{T}} - 0 \\ & = 4\langle \nabla \cdot 2\mu(\mathbf{u}\mathbf{S}) \rangle_{\mathrm{T}} - 4\langle \nabla \mathbf{u} : 2\mu\mathbf{S} \rangle_{\mathrm{T}}, \quad (\mathrm{A} 25) \end{aligned} $$

By periodicity, the first term on RHS of Eq. (A 25) cancels out, and remains only the kinetic energy dissipation rate which reads $\langle\epsilon\rangle_{\mathrm{T}} := \langle\nabla\mathbf{u} : 2\mu\mathbf{S}\rangle_{\mathrm{T}} = \langle 2\mu\mathbf{S} : \mathbf{S}\rangle_{\mathrm{T}}$ (see e.g. Wilcox 1998, p 251). The other terms of the total kinetic energy transport equation can be retrieved similarly to finally obtain the one-point total kinetic energy budget, which for the present configuration writes:

$$ \frac{d}{dt}\langle k\rangle_{\mathrm{T}} = \langle F\rangle_{\mathrm{T}} - \langle \epsilon\rangle_{\mathrm{T}} + \langle S\rangle_{\mathrm{T}}, \quad (\mathrm{A} 26) $$

where $\langle F\rangle_{\mathrm{T}} := \langle \mathbf{u} \cdot \mathbf{F} \rangle_{\mathrm{T}}$ and $\langle S\rangle_{\mathrm{T}} := \langle \mathbf{u} \cdot \mathbf{S} \rangle_{\mathrm{T}}$ represent the contribution of forcing and surface tension to the kinetic energy budget, respectively. Note that, in statistically steady flows, the time derivative cancels out together with surface tension term. The latter can indeed be written as (Dodd & Ferrante 2016; Trefftz-Posada & Ferrante 2023)

$$ \langle S\rangle_{\mathrm{T}} \sim -\sigma d_t A_{\Gamma}, \quad (\mathrm{A} 27) $$

which is zero since the surface area of the liquid/gas interface $A_{\Gamma}$ is on average constant. The other terms vanish because the flow under consideration is periodic.

The same reasoning can be applied to the different terms of the conditionally averaged KHM equation Eq. (A 21). Here again, in the limit of large separations, each term of Eq. (A 21) tends towards four times their counterpart in the one-point conditionally averaged kinetic energy budget derived by Dodd & Ferrante (2016). For instance, the two-point viscous term conditionally averaged in either the liquid or gas phase has the following limit,

$$ \begin{aligned} & \lim_{r \to \infty} 2\langle(\delta\mathbf{u}) \cdot (\delta(\nabla \cdot 2\mu\mathbf{S}))\rangle_{\mathrm{C}} = 4\langle\mathbf{u} \cdot (\nabla \cdot 2\mu\mathbf{S})\rangle_{\mathrm{C}} \\ & = -8\mu_{\mathrm{C}}\langle\nabla\mathbf{u} : \mathbf{S}\rangle_{\mathrm{C}} + 8\mu_{\mathrm{C}}\langle\nabla \cdot (\mathbf{u}\mathbf{S})\rangle_{\mathrm{C}} \\ & = -8\mu_{\mathrm{C}}\langle\mathbf{S} : \mathbf{S}\rangle_{\mathrm{C}} + 8\mu_{\mathrm{C}}\langle\nabla \cdot (\mathbf{u}\mathbf{S})\rangle_{\mathrm{C}} \\ & := -4\langle\epsilon\rangle_{\mathrm{C}} + 4\langle T_{\nu}\rangle_{\mathrm{C}}, \quad (\mathrm{A} 28) \end{aligned} $$

where the kinetic energy dissipation rate in phase $C \in \{\mathbb{L}, \mathbb{G}\}$ is $\langle\epsilon\rangle_{\mathrm{C}} := 2\mu_{\mathrm{C}}\langle\mathbf{S} : \mathbf{S}\rangle_{\mathrm{C}}$ and the viscous transport is $\langle T_{\nu}\rangle_{\mathrm{C}} := 2\mu_{\mathrm{C}}\langle\nabla \cdot (\mathbf{u}\mathbf{S})\rangle_{\mathrm{C}}$. The conditional one-point average of a quantity $\bullet$ is defined by:

$$ \langle \bullet \rangle_{\mathrm{C}} := \frac{\langle \phi_{\mathrm{C}} \bullet \rangle_{\mathrm{T}}}{\langle \phi_{\mathrm{C}} \rangle_{\mathrm{T}}}. \quad (\mathrm{A} 29) $$

Proceeding with the pressure gradient term in Eq. (A 21) yields

$$ \lim_{r \to \infty} -2\langle\nabla_{\mathbf{X}} \cdot (\delta\mathbf{u})(\delta P)\rangle_{\mathrm{C}} = -4\langle\nabla \cdot \mathbf{u}P\rangle_{\mathrm{C}} := 4\langle T_p\rangle_{\mathrm{C}}, \quad (\mathrm{A} 30) $$

again with $C \in \{\mathbb{L}, \mathbb{G}\}$. These terms are constitutive of the one-point conditional kinetic energy budget (see Eqs. (3.8) to (3.11) in Dodd & Ferrante 2016) and are
recovered here from the KHM equation in the limit of large separations. In presence of a forcing term, the latter writes:

$$ \frac{d}{dt} \langle k \rangle_{\mathcal{C}} = \langle F \rangle_{\mathcal{C}} - \langle \epsilon \rangle_{\mathcal{C}} + \langle T_{\nu} \rangle_{\mathcal{C}} + \langle T_{p} \rangle_{\mathcal{C}}, \quad (A 31) $$

For obtaining the above equation, it was necessary to note that,

$$ \langle \partial_t k \rangle_{\mathcal{C}} + \langle \nabla \cdot \mathbf{u} k \rangle_{\mathcal{C}} = \frac{d}{dt} \langle k \rangle_{\mathcal{C}} - \int_{\Gamma} k \mathbf{u}_b \cdot \mathbf{n} \, dS + \int_{\Gamma} k \mathbf{u} \cdot \mathbf{n} \, dS, \quad (A 32) $$

where we have employed the Reynolds transport theorem for the time derivative term and the Green-Ostrogadski theorem for the transport of kinetic energy. In case of immiscible fluids with no phase change, one has $\mathbf{u}_b = \mathbf{u}$ at the interface between the two phases, and thus the sum of the time derivative term and the transport term reduces to $d\langle k \rangle_{\mathcal{C}}/dt$.

Let us now focus on the large-scales asymptotic behaviour when $\mathcal{C} = \mathbb{M}$. For this, it is first worth noting that for any one-point statistics ●, one has (Dodd & Ferrante 2016)

$$ \langle \bullet \rangle_{\mathcal{T}} = \langle \phi_L \rangle_{\mathcal{T}} \langle \bullet \rangle_{\mathcal{L}} + \langle \phi_G \rangle_{\mathcal{T}} \langle \bullet \rangle_{\mathcal{G}}. \quad (A 33) $$

In addition, when the separation $r$ goes to infinity (Torquato 2002),

$$ \lim_{r \to \infty} \gamma_L = \langle \phi_L \rangle_{\mathcal{T}}^2 \quad (A 34a) $$

$$ \lim_{r \to \infty} \gamma_G = \langle \phi_G \rangle_{\mathcal{T}}^2 \quad (A 34b) $$

$$ \lim_{r \to \infty} \gamma_M = \langle \phi_L \rangle_{\mathcal{T}} \langle \phi_G \rangle_{\mathcal{T}} \quad (A 34c) $$

Then using Eq. (A 20) (with $\mathcal{C} = \mathbb{M}$), Eqs. (A 16), (A 17), (A 18), (A 34), and after some straightforward manipulations, it can be obtained that for e.g. the forcing term:

$$ \lim_{r \to \infty} +2\langle(\delta\mathbf{u}) \cdot (\delta\mathbf{F})\rangle_{\mathbb{M}} - \langle\mathcal{C}(\mathbf{F})\rangle_{\mathbb{M}} = 4\frac{\langle F \rangle_{\mathcal{L}} + \langle F \rangle_{\mathcal{G}}}{2} \quad (A 35) $$

which is thus simply (four times) the average of $\langle F \rangle$ in the liquid and gas phase. The same applies to the other terms.

### A.5. Limit at small separations

We are interested in the limit of the two-point kinetic energy at small separations. For this purpose, we use the following decomposition:

$$ \langle (\delta \rho \mathbf{u}) \cdot (\delta \mathbf{u}) \rangle_{\mathcal{C}} = \underbrace{\langle \delta \rho \mathbf{u} (\mathbf{x} + \mathbf{r}) \cdot \delta \mathbf{u} \rangle_{\mathcal{C}}}_{T_1} + \underbrace{\langle \rho(\mathbf{x}) | \delta \mathbf{u} |^2 \rangle_{\mathcal{C}}}_{T_2}. \quad (A 36) $$

Let us start with term $T_2$. In the limit of small separations, we can use Taylor expansions to write

$$ |\delta \mathbf{u}|^2 \sim r^2 (\nabla \mathbf{u} : \nabla \mathbf{u}) + \mathcal{O}(r^3). \quad (A 37) $$

Decomposing the velocity gradient tensor into symmetric and antisymmetric components, it is then possible to obtain

$$ T_2 = \frac{r^2}{3} \left\langle \rho \left( \mathbf{S} : \mathbf{S} + \frac{1}{2} |\omega|^2 \right) \right\rangle_{\mathcal{C}} + \mathcal{O}(r^3). \quad (A 38) $$

where $\omega$ is the vorticity vector. The $1/3$ coefficient is obtained under the assumption of statistical isotropy.
For term $T_1$, we start by recalling that the increment $\delta\rho$ is active only when the two points cross the interface, i.e. only when $\mathbb{C} = \mathbb{T}$. The density increment $\delta\rho$ is given by:

$$ \delta\rho = (\rho_L - \rho_G)\delta_\Gamma(\mathbf{x} - \mathbf{x}_s)|\mathbf{r} \cdot \mathbf{n}(\mathbf{x}_s)|. \quad (A 39) $$

Here, $\delta_\Gamma(\mathbf{x} - \mathbf{x}_s)|\mathbf{r} \cdot \mathbf{n}(\mathbf{x}_s)|$ is the probability that vector $\mathbf{r}$ crosses the interface. This would lead to

$$ \langle \delta\rho \rangle = \frac{\Delta\rho}{V} \int_{\Gamma} |\mathbf{r} \cdot \mathbf{n}| dA = \Delta\rho \frac{A_{\Gamma} r}{V^2} + \mathcal{O}(r^3), \quad (A 40) $$

where $\Delta\rho = \rho_L - \rho_G$ and $V$ is the averaging volume for $\mathbb{C} = \mathbb{T}$. The factor 1/2 comes from the isotropic assumption for the orientation of the surface. This result is a generalization in 3D-space of the Buffon’s needle problem (Torquato 2002; Thiesset et al. 2021; Thiesset & Federrath 2023; Gauding et al. 2022). Higher order terms would include the effect of mean and Gaussian curvature.

Decomposing further $T_1$ as:

$$ T_1 = \langle \delta\rho (\mathbf{u}(\mathbf{x}) \cdot \delta\mathbf{u} + |\delta\mathbf{u}|^2) \rangle_{\mathbb{C}}, \quad (A 41) $$

we realize that $\delta\rho|\delta\mathbf{u}|^2$ is third-order in $r$ since $\delta\rho \sim r$ while $|\delta\mathbf{u}|^2 \sim r^2$. Therefore, only $\mathbf{u}(\mathbf{x}) \cdot \delta\mathbf{u}$ contributes to leading order. Taylor series expansion for $\mathbf{u} \cdot \delta\mathbf{u}$ yields:

$$ T_1 = \Delta\rho \langle \delta(\mathbf{x} - \mathbf{x}_s) |\mathbf{r} \cdot \mathbf{n}(\mathbf{x}_s)| \mathbf{u} \cdot (\mathbf{r} \cdot \nabla) \mathbf{u} \rangle_{\mathbb{T}} + \mathcal{O}(r^3). \quad (A 42) $$

The term $\mathbf{u} \cdot (\mathbf{r} \cdot \nabla)\mathbf{u}$ can further be rewritten as $\frac{1}{2}\nabla_r|\mathbf{u}|^2$ and therefore, $T_1$ can be expressed as:

$$ T_1 = \frac{\Delta\rho r^2}{V^2} \int_{\Gamma} |\hat{\mathbf{r}} \cdot \mathbf{n}| \nabla_{\hat{\mathbf{r}}} |\mathbf{u}|^2 dA + \mathcal{O}(r^3), \quad (A 43) $$

where $\hat{\mathbf{r}} = \mathbf{r}/r$. Further, it is reasonable (although not proven rigorously) to consider that $\nabla_{\hat{\mathbf{r}}} |\mathbf{u}|^2$ is independent of $\mathbf{n}$. With this assumption, one gets:

$$ T_1 = \frac{\Delta\rho r^2}{V^2} \int_{\Gamma} \nabla_{\hat{\mathbf{r}}} |\mathbf{u}|^2 dA + \mathcal{O}(r^3). \quad (A 44) $$

Because $\int_{\Gamma} |\hat{\mathbf{r}} \cdot \mathbf{n}| dA = 1/2$. By use of isotropy, the gradient $\nabla_{\hat{\mathbf{r}}} |\mathbf{u}|^2$ does not depend on the orientation of $\hat{\mathbf{r}}$. Therefore, there is no preferred direction of $|\mathbf{u}|^2$ with respect to $\hat{\mathbf{r}}$. In addition, in statistically homogeneous flows with no mean flow, the mean velocity and the mean gradients should be zero. Furthermore, for immiscible fluids, the velocity $\mathbf{u}$ is continuous across the interface so that there is no net interfacial exchanges between phases. These arguments constitute a sufficient body of evidence to state that $T_1$ vanishes. Pending a more rigorous proof, we will continue with $T_1 = 0$ for $\mathbb{C} \equiv \mathbb{T}$ as well.

In summary, the small scale limit of the mixed structure function writes:

$$ \langle (\delta\rho\mathbf{u}) \cdot (\delta\mathbf{u}) \rangle_{\mathbb{C}} = \frac{r^2}{3} \left( \left\langle \frac{\rho\epsilon}{2\mu} \right\rangle_{\mathbb{C}} + \left\langle \frac{\rho|\omega|^2}{2} \right\rangle_{\mathbb{C}} \right) + \mathcal{O}(r^3), \quad (A 45) $$

for any $\mathbb{C} \in \{\mathbb{L}, \mathbb{G}, \mathbb{T}\}$. Recall that $\epsilon/2\mu = \mathbf{S} : \mathbf{S}$. Note also that for $\mathbb{C} \in \{\mathbb{L}, \mathbb{G}\}$, the fluid physical properties are constant per phase. Under the assumption of statistical homogeneity, one then has:

$$ \left\langle \frac{\rho\epsilon}{\mu} \right\rangle_{\mathbb{C}} = \langle \rho | \omega |^2 \rangle_{\mathbb{C}}. \quad (A 46) $$
Because this relation holds true per phase, it should also apply to the total field by simple addition of the statistics in the liquid and gas phase. The small-scale limit of the two-point kinetic energy can thus be rewritten as:

$$ \langle (\delta \rho \mathbf{u}) \cdot (\delta \mathbf{u}) \rangle_{\mathcal{C}} = \frac{r^2}{3} \langle \rho |\omega|^2 \rangle_{\mathcal{C}} + \mathcal{O}(r^3) \quad (A 47) $$

$$ = \frac{r^2}{3} \left\langle \frac{\rho \epsilon}{\mu} \right\rangle_{\mathcal{C}} + \mathcal{O}(r^3), \quad (A 48) $$

for any $\mathcal{C} \in \{\mathcal{L}, \mathcal{G}, \mathcal{T}\}$. Eq. (A 48) reveals that, in the limit of small separations, the unconditional and conditional two-point kinetic energy is quadratic in $r$ and is proportional to the enstrophy (or kinetic energy dissipation rate).

In single phase flows in absence of solid boundaries, the evolution equation for the enstrophy field is obtained from the Navier-Stokes equation. When averaged over a domain $\mathcal{C}$, one obtains

$$ d_t \frac{\langle |\omega|^2 \rangle_{\mathcal{C}}}{2} = \langle \omega \omega : \nabla \mathbf{u} \rangle_{\mathcal{C}} + \nu \left\langle \nabla^2 \frac{|\omega|^2}{2} - \nabla \omega : \nabla \omega \right\rangle_{\mathcal{C}} + \langle \omega \cdot (\nabla \times \mathbf{F}) \rangle_{\mathcal{C}} \quad (A 49) $$

Equation (A 49) reveals that the time variation of the averaged enstrophy depends on a production term due to vortex stretching (the first term on RHS of Eq A 49), on the enstrophy viscous diffusion and viscous dissipation (sometimes referred to as the palinstrophy) and on the effect of artificial forcing as represented by the rightmost term on RHS of Eq. (A 49).

Given Eq. (A 48), the limit at small scales of the KHM equation per phase should tend towards the transport equation for the enstrophy, Eq. (A 49). In single-phase flows, this was done for instance by Antonia et al. (2000). Proceeding similarly here, the transport term in scale-space in the KHM equation (Eq. A 21) tends towards the production by vortex stretching in Eq. (A 49), the viscous term in Eq. (A 21) tends towards the viscous terms in Eq. (A 49), while the forcing term in (A 49) is related to its counterpart in Eq. (A 21).

For the case $\mathcal{C} = \mathcal{T}$, the dynamics of vorticity presents a considerably higher level of complexity. Terrington et al. (2022) showed that conservation laws for the volume integrated vorticity reveal some additional source terms which come from surface tension, viscous stresses and pressure gradients. This was achieved through sophisticated algebraic manipulations, which we have not yet been able to leverage to derive an evolution equation for the enstrophy in multiphase systems. Deriving an equation for the enstrophy starting from the KHM equation in the limit of small separations is not an easy task either. Indeed, as for $\delta(\rho\mathbf{u})$, some terms are discontinuous at the interface (e.g. the pressure and viscous terms), which require the same treatment as term $T_1$ above. Handling the derivation for all these terms is left for future analysis.

# Appendix B. Numerical tools and post-processing procedures

## B.1. Time advancement algorithm

We here describe the fastRK3 algorithm used for the time advancement of the numerical solutions. The first step of the algorithm is to transport the coupled level-set and VOF fields at time $t + \Delta t$ and compute the related density $\rho_{n+1}$ and viscosity $\mu_{n+1}$ at time $t + \Delta t$. For ensuring consistency between the transport of VOF and momentum, the convective term in Eq. (2.3) is also computed.
Let us now note $\mathbf{R}$ the right-hand side of Eq. (2.4). Formally, the fastRK3 algorithm is similar to the classical RK3 method. It consists in advancing from the velocity at time $t$ noted $\mathbf{u}_n$ to the velocity at time $t + \Delta t$ noted $\mathbf{u}_{n+1}$, by calculating three intermediate velocity fields, viz

$$ \mathbf{u}_1^* = \mathbf{u}_n + \frac{\Delta t}{3} \mathbf{R}(\mathbf{u}_n) \qquad (B1a) $$

$$ \mathbf{u}_1 = \mathbf{u}_1^* - \frac{\Delta t}{3} \frac{\nabla P_1}{\rho_{n+1}} \qquad (B1b) $$

$$ \mathbf{u}_2^* = \mathbf{u}_n + \Delta t [-\mathbf{R}(\mathbf{u}_n) + 2\mathbf{R}(\mathbf{u}_1)] \qquad (B1c) $$

$$ \mathbf{u}_2 = \mathbf{u}_2^* - \Delta t \frac{\nabla P_2}{\rho_{n+1}} \qquad (B1d) $$

$$ \mathbf{u}_3^* = \mathbf{u}_n + \Delta t \left[ \frac{3}{4} \mathbf{R}(\mathbf{u}_1) + \frac{1}{4} \mathbf{R}(\mathbf{u}_2) \right] \qquad (B1e) $$

$$ \mathbf{u}_{n+1} = \mathbf{u}_3 = \mathbf{u}_3^* - \Delta t \frac{\nabla P_3}{\rho_{n+1}} \qquad (B1f) $$

In the classical RK3 method, the intermediate pressure fields $P_i \in \{P_1, P_2, P_3\}$ are obtained by solving the Poisson equation:

$$ \frac{\nabla^2 P_i}{\rho_{n+1}} = \nabla \cdot \mathbf{u}_i^* \qquad (B2) $$

at all 3 steps of the algorithm. Instead, the fastRK3 method expresses the intermediate pressures $P_1$ and $P_2$ as extrapolation of the pressure at time $t + \Delta t/3$ and $t + \Delta t$, respectively, using $P_n$ and $P_{n-1}$, viz.

$$ P_1 = \frac{5}{3} P_n - \frac{2}{3} P_{n-1} \qquad (B3a) $$

$$ P_2 = 2P_n - P_{n-1}. \qquad (B3b) $$

This is referred to as the mid-point extrapolation in Aithal et al. (2023). By doing this, only one Poisson equation is solved in order to compute the pressure at time $t + \Delta t$, noted $P_3 = P_{n+1}$, while 3 are required for the classical RK3 algorithm (one per sub-step of the Runge-Kutta scheme). Here, $P_3$ is obtained by using a multi-grid preconditioned Conjugate Gradient algorithm (MGCG) (Zhang 1996). We do not use the decomposition proposed by Dodd & Ferrante (2016); Trefftz-Posada & Ferrante (2023) in order to solve a Poisson equation with constant coefficients. Note that the fastRK3 is not self-starting. Hence, in order to get $P_n$ and $P_{n-1}$, the first 3 steps of time advancement are done using the classical RK3 scheme before switching to fastRK3 for the rest of the simulation.

## B.2. Forcing procedure

We here provide the details of the forcing procedure used to maintain turbulence at steady state. The forcing acts only at large scales, in a spectral band $[k_{\min} = 2\pi/L, k_{\max} = 2\pi/(L/3)]$ where $L$ is the width of the simulation domain. We use the same forcing spectrum as Federrath et al. (2010), noted $W(k)$, which corresponds to a paraboloid peaking around $k_f = (k_{\max} + k_{\min})/2$, viz.

$$ W(k) = \begin{cases} \frac{k_f}{k} [1 - K(k - k_f)^2] & \text{where } k_{\min} < k < k_{\max} \\ 0 & \text{elsewhere,} \end{cases} \qquad (B4) $$
where $K = 4/(k_{max} - k_{min})^2$. In order to maintain a statistically steady state with a prescribed kinetic energy, the amplitude of the forcing is adjusted using a PID controller, where the amplification factor of the controller, noted A, is given by

$$A = G \left( \Upsilon(t) + \frac{1}{\tau_i} \int_{0}^{\tau} \Upsilon(t) dt + \tau_d \frac{d}{dt} \Upsilon(t) \right), \quad (B5)$$

where $\Upsilon := \tilde{k}(t) - \tilde{k}_t$ is the error between the kinetic energy at time $t$, i.e. $\tilde{k}(t) = \langle k \rangle_T / \langle \rho \rangle_T$ and the target kinetic energy $\tilde{k}_t$. The controller parameters were set to $G = 200$, $\tau_i = 0.25s$ and $\tau_d = 0.01s$. The correlation time of the stochastic forcing is equal to $L_f/u'$, where $L_f = 2\pi/k_f = L/2$ and $u' = (2\tilde{k}_t)^{1/2}$.

This forcing was chosen because it has some interesting features which are detailed below. First, because it acts at large scales, (i) the intermediate and small scales are believed to behave more naturally, and (ii) the approach towards the asymptotic behaviour at infinite Reynolds number is faster than decaying turbulence or by using a linear forcing (Antonia & Burattini 2006). Further, because this forcing is based on a stochastic process, its contribution to the kinetic energy budget is expected (iii) to be independent of the chosen flow physical parameters and (iv) to be the same in the liquid and gas phase. This was confirmed later by analysing the scale-by-scale budgets. Thirdly, by the use of a PID controller, we minimize the time fluctuations of energy injection, thereby (v) mitigating what is referred to as non-equilibrium effects (Goto & Vassilicos 2015; Fang & Bos 2023).

## B.3. Numerical resolution

To appraise the appropriateness of our numerical solutions, we have monitored two metrics. The first is the difference between the total time-averaged injection and dissipation of kinetic energy. Given Eq. (2.13), injection and dissipation should be on average equal, except if numerical errors (numerical dissipation or dispersion) are at play. In case of insufficient resolutions we have noted that both the transport term and the surface tension term were negative, thereby contributing to the kinetic energy budget as non-physical numerical dissipation. Their contribution were though converging to zero when increasing the resolution. The second metric for inferring the appropriate resolution is the conservation of mass. In case of insufficient resolution, there could be some inconsistency between the geometrical properties obtained from the VOF field and the level-set field. In our solver, if such situations occur, the VOF is readjusted to match the geometrical properties of the level-set, at the price of small mass variations.

In order to keep these metrics of the order of a percent, the viscosity $\nu$ and surface tension $\sigma$ were chosen so that the ratio of the Kolmogorov length-scale $\eta$ to the mesh size $dx$ is about 2 (see Table 1) and the Weber number is not more than 50. While half this resolution is likely sufficient for low density ratios (say $R_\rho \le 5$), we have found that such a fine resolution was necessary for larger density ratio. This results in significantly smaller Taylor scale Reynolds number $R_\lambda$ compared to what can be found in the literature using the same number of simulation points (Crialesi-Esposito et al. 2022; Cannon et al. 2024).

## B.4. Post-processing

Statistics are computed using the post-processing library PyArcher which comes as a submodule of the archer code. Instead of computing the explicit form of the different terms of the KHM equation which are given in Appendix (Eqs. (A 12), (A 20) and Eqs.
(A 21)), we use Eq. (2.2) where all terms of the NS equation are extracted directly from the code. This was found to be the most computationally efficient and accurate. The two-point statistics are computed using homemade routines which work on GPU using OpenACC directives.

In previous studies from the literature (Dodd & Ferrante 2016; Crialesi-Esposito et al. 2022; Treftz-Posada & Ferrante 2023), it is not mentioned how conditional averages are computed. Here, the phase indicator is defined from the distance function used in the level-set method which we note $\Psi$. However, since the code is based on a staggered Cartesian grid, we have decided to interpolate the level-set field (a cell-centred quantity) at the same position as the velocity components, i.e. at the faces of each computational cell. Three different face-centred level-set fields are thus obtained, one for each velocity component. These are noted $\Psi_s^i$ for each direction $i$. The phase indicator function for the liquid phase can subsequently be defined as $\phi_L = 1$ where $\Psi_s^i > 0.5dx$ and for the gas phase, $\phi_G = 1$ where $\Psi_s^i < -0.5dx$. This shift of $0.5dx$ either in the liquid or gas phase was imposed in order to exclude the cells where the interface (and hence the surface tension term) is present.

The budget of the KHM equation is closed up to about $10^{-3}\epsilon_f$. As a projection method is used to enforce the incompressibility condition, the resolution of the NS equation is bounded by the residual of the Poisson equation, which is of the order of $L_2(\nabla \cdot u) \sim 10^{-9}$. This constitutes a small yet finite source of error. The second source of error arises from the transport term in the NS and therefore the KHM equations: in our code, we solve the transport term for $\rho u$ using the Rudman scheme, which needs to be converted to get the transport term for $u$. This operation is likely to come with some error which remains complicated to evaluate. The last, and probably the most dominant, loss of precision can be rooted to the angular averaging operation which requires to interpolate the Cartesian $(r_x, r_y, r_z)$ coordinates to the spherical coordinates $(r, \theta, \varphi)$ as described in (2.6). More details on this operation can be found in Thiesset & Poux (2020).

## REFERENCES

* AITHAL, AB, TIPIRNENI, M & FERRANTE, A 2023 Temporal accuracy of FastRK3. *J. Comput. Phys.* **475**, 111853.

* ALUIE, H 2011 Compressible turbulence: the cascade and its locality. *Phys. Rev. Lett.* **106** (17), 174502.

* ALUIE, H 2013 Scale decomposition in compressible turbulence. *Phys. D* **247** (1), 54–65.

* ANTONIA, RA & BURATTINI, P 2006 Approach to the 4/5 law in homogeneous isotropic turbulence. *J. Fluid Mech.* **550**, 175–184.

* ANTONIA, RA, ZHOU, T, DANAILA, L & ANSELMET, F 2000 Streamwise inhomogeneity of decaying grid turbulence. *Phys. Fluids* **12** (11), 3086–3089.

* ARUN, S, SAMEEN, A, SRINIVASAN, B & GIRIMAJI, S S 2021 Scale-space energy density function transport equation for compressible inhomogeneous turbulent flows. *J. Fluid Mech.* **920**, A31.

* BALACHANDAR, S & EATON, JOHN K 2010 Turbulent dispersed multiphase flow. *Annu. Rev. Fluid Mech.* **42**, 111–133.

* BANERJEE, S & GALTIER, S 2013 Exact relation with two-point correlation functions and phenomenological approach for compressible magnetohydrodynamic turbulence. *Phys. Rev. E* **87** (1), 013019.

* BOUKHARFANE, R, ER-RAIY, A, PARSANI, M & CHAKRABORTY, N 2021 Structure and dynamics of small-scale turbulence in vaporizing two-phase flows. *Scientific reports* **11** (1), 15242.

* BRAHAMI, Y 2020 Active scalar mixing in turbulent jet flows. *Thesis, Normandie Université*.
BUNNER, B & TRYGGVASON, G 2002 Dynamics of homogeneous bubbly flows. Part 2. velocity fluctuations. *J. Fluid Mech.* **466**, 53–84.

CANNON, I., SOLIGO, G. & ROSTI, M. E. 2024 Morphology of clean and surfactant-laden droplets in homogeneous isotropic turbulence. *J. Fluid Mech.* **987**, A31.

CRIALESI-ESPOSITO, M, BOFFETTA, G, BRANDT, L, CHIBBARO, S & MUSACCHIO, S 2023a Intermittency in turbulent emulsions. *J. Fluid Mech.* **972**, A37.

CRIALESI-ESPOSITO, M., CHIBBARO, S. & BRANDT, L. 2023b The interaction of droplet dynamics and turbulence cascade. *Communications Physics* **6** (1), 5.

CRIALESI-ESPOSITO, M, ROSTI, ME, CHIBBARO, S & BRANDT, L 2022 Modulation of homogeneous and isotropic turbulence in emulsions. *J. Fluid Mech.* **940**, A19.

DANAILA, L, ANTONIA, RA & BURATTINI, P 2012a Comparison between kinetic energy and passive scalar energy transfer in locally homogeneous isotropic turbulence. *Phys. D* **241** (3), 224–231.

DANAILA, L, KRAWCZYNSKI, JF, THIESSET, F & RENOU, B 2012b Yaglom-like equation in axisymmetric anisotropic context. *Phys. D* **241**, 216–223.

DAVIDSON, PA & PEARSON, BR 2005 Identifying turbulent energy distributions in real, rather than Fourier, space. *Phys. Rev. Lett.* **95** (21), 214501.

DODD, MS & FERRANTE, A 2016 On the interaction of Taylor length scale size droplets and isotropic turbulence. *J. Fluid Mech.* **806**, 356–412.

DURET, B, LURET, G, REVEILLON, J, MÉNARD, T, BERLEMONT, A & DEMOULIN, FX 2012 DNS analysis of turbulent mixing in two-phase flows. *Int. J. Multiphase Flow* **40**, 93–105.

ESTIVALEZES, J-L, ANISZEWSKI, W, AUGUSTE, F, LING, Y, OSMAR, L, CALTAGIRONE, J-P, CHIRCO, L, PEDRONO, A, POPINET, S, BERLEMONT, A & OTHERS 2022 A phase inversion benchmark for multiscale multiphase flows. *J. Comput. Phys.* **450**, 11080.

ESWARAN, V. & POPE, S. B. 1988 An examination of forcing in direct numerical simulations of turbulence. *Comput and Fluids* **16** (3), 257–278.

EYINK, GL & DRIVAS, TD 2018 Cascades and dissipative anomalies in compressible fluid turbulence. *Phys. Rev. X* **8** (1), 011022.

FANG, L & BOS, W 2023 An EDQNM study of the dissipation rate in isotropic non-equilibrium turbulence. *J. Turbul.* **24** (6-7), 217–234.

FEDERRATH, C, KLESSEN, R S, IAPICHINO, L & BEATTIE, J R 2021 The sonic scale of interstellar turbulence. *Nature Astronomy* **5** (4), 365–371.

FEDERRATH, C, ROMAN-DUVAL, J, KLESSEN, RS, SCHMIDT, W & MAC LOW, M-M 2010 Comparing the statistics of interstellar turbulence in simulations and observations- solenoidal versus compressive turbulence forcing. *Astronomy & Astrophysics* **512**, A81.

FEDKIW, R P, ASLAM, T, MERRIMAN, B & OSHER, S 1999 A non-oscillatory Eulerian approach to interfaces in multimaterial flows (the ghost fluid method). *J Comput Phys* **152** (2), 457–492.

FERRAND, R, GALTIER, S, SAHRAOUI, F & FEDERRATH, C 2020 Compressible turbulence in the interstellar medium: New insights from a high-resolution supersonic turbulence simulation. *Astrophys. J.* **904** (2), 160.

FREUND, A. & FERRANTE, A. 2019 Wavelet-spectral analysis of droplet-laden isotropic turbulence. *J. Fluid Mech.* **875**, 914–928.

FRISCH, H. L. & STILLINGER, F. H. 1963 Contribution to the statistical geometric basis of radiation scattering. *J. Chem. Phys.* **38** (9), 2200–2207.

GALTIER, S & BANERJEE, S 2011 Exact relation for correlation functions in compressible isothermal turbulence. *Phys. Rev. Lett.* **107** (13), 134501.

GAUDING, M, BODE, M, BRAHAMI, Y, VAREA, É & DANAILA, L 2021 Self-similarity of turbulent jet flows with internal and external intermittency. *J. Fluid Mech.* **919**, A41.

GAUDING, M, DANAILA, L & VAREA, E 2018 One-point and two-point statistics of homogeneous isotropic decaying turbulence with variable viscosity. *Int. J. Heat Fluid Flow* **72**, 143–150.

GAUDING, M, THIESSET, F, VAREA, E & DANAILA, L 2022 Structure of iso-scalar sets. *J. Fluid Mech.* **942**, A14.

GOTO, S. & VASSILICOS, J-C. 2015 Energy dissipation and flux laws for unsteady turbulence. *Phys. Lett. A* **379** (16-17), 1144–1148.

HAMBA, F 2015 Turbulent energy density and its transport equation in scale space. *Phys. Fluids* **27** (8).
HAMBA, F 2018 Turbulent energy density in scale space for inhomogeneous turbulence. *J. Fluid Mech.* **842**, 532–553.

HELLINGER, P, PAPINI, E, VERDINI, A, LANDI, S, FRANCI, L, MATTEINI, L & MONTAGUD-CAMPS, V 2021a Spectral transfer and Kármán-Howarth-Monin equations for compressible hall magnetohydrodynamics. *Astrophys. J.* **917** (2), 101.

HELLINGER, P, VERDINI, A, LANDI, S, PAPINI, E, FRANCI, L & MATTEINI, L 2021b Scale dependence and cross-scale transfer of kinetic energy in compressible hydrodynamic turbulence at moderate Reynolds numbers. *Phys. Rev. Fluids* **6** (4), 044607.

HILL, RJ 2002 Exact second-order structure-function relationships. *J. Fluid Mech.* **468**, 317–326.

HINZE, J. O. 1955 Fundamentals of the hydrodynamic mechanism of splitting in dispersion processes. *AIChE Journal* **1** (3), 289–295.

INNOCENTI, A, JACCOD, A, POPINET, S & CHIBBARO, S 2021 Direct numerical simulation of bubble-induced turbulence. *J. Fluid Mech.* **918**, A23.

KIDA, S & ORSZAG, SA 1990 Energy and spectral dynamics in forced compressible turbulence. *J. Sci. Comput.* **5**, 85–125.

KIRSTE, R. & POROD, G. 1962 Röntgenkleinwinkelstreuung an kolloiden systemen asymptotisches verhalten der streukurven. *Kolloid-Zeitschrift und Zeitschrift für Polymere* **184** (1), 1–7.

KOLLA, H., HAWKES, E. R., KERSTEIN, A. R., SWAMINATHAN, N. & CHEN, J. H. 2014 On velocity and reactive scalar spectra in turbulent premixed flames. *J Fluid Mech* **754**, 456–487.

KOLMOGOROV, A. N. 1949 On the breakage of drops in a turbulent flow. *Dokl. Akad. Navk. SSSR* **66**, 825–828.

KRITSUK, AG, WAGNER, R & NORMAN, ML 2013 Energy cascade and scaling in supersonic isothermal turbulence. *J. Fluid Mech.* **729**, R1.

LAI, CCK, CHARONKO, JJ & PRESTRIDGE, K 2018 A Kármán-Howarth-Monin equation for variable-density turbulence. *J. Fluid Mech.* **843**, 382–418.

LANCE, M & BATAILLE, J 1991 Turbulence in the liquid phase of a uniform bubbly air-water flow. *J. Fluid Mech.* **222**, 95–118.

LING, Y, FUSTER, D, TRYGGVASON, G & ZALESKI, S 2019 A two-phase mixing layer between parallel gas and liquid streams: multiphase turbulence statistics and influence of interfacial instability. *J. Fluid Mech.* **859**, 268–307.

LOISY, A & NASO, A 2017 Interaction between a large buoyant bubble and turbulence. *Phys. Rev. Fluids* **2** (1), 014606.

LU, J & TRYGGVASON, G 2018 Direct numerical simulations of multifluid flows in a vertical channel undergoing topology changes. *Phys. Rev. Fluids* **3** (8), 084401.

LU, J & TRYGGVASON, G 2019 Multifluid flows in a vertical channel undergoing topology changes: Effect of void fraction. *Phys. Rev. Fluids* **4** (8), 084301.

LUCCI, F, FERRANTE, A & ELGHOBASHI, S 2010 Modulation of isotropic turbulence by particles of Taylor length-scale size. *J. Fluid Mech.* **650**, 5–55.

LUNDGREN, T S 2003 Linearly forced isotropic turbulence. *Center for Turbulence Research Annual Research Briefs 2003*.

MCCASLIN, JO & DESJARDINS, O 2014 Theoretical and computational modeling of turbulence/interface interactions. In *Proceedings of the Summer Program*, p. 79.

MÉNARD, T, TANGUY, S & BERLEMONT, A 2007 Coupling level set/VOF/ghost fluid methods: Validation and application to 3D simulation of the primary break-up of a liquid jet. *Int. J. Multiphase Flow* **33** (5), 510–524.

MUKHERJEE, S., SAFDARI, A., SHARDT, O., KENJEREŠ, S. & VAN DEN AKKER, H.E.A. 2019 Droplet-turbulence interactions and quasi-equilibrium dynamics in turbulent emulsions. *J. Fluid Mech.* **878**, 221–276.

NARULA, H, PANDEY, V, MITRA, D & PERLEKAR, P 2025 Scale-by-scale energy transfers in bubbly flows. *arXiv preprint arXiv:2506.07693*.

NGUYEN, M Q, DELACHE, A, SIMOËNS, S, BOS, W & EL HAJEM, M 2016 Small scale dynamics of isotropic viscoelastic turbulence. *Physical Review Fluids* **1** (8), 083301.

NI, R 2024 Deformation and breakup of bubbles and drops in turbulence. *Annu. Rev. Fluid Mech.* **56** (1), 319–347.

NIE, Q & TANVEER, S 1999 A note on third-order structure functions in turbulence. *Proc. Roy. Soc. 455* (1985), 1615–1635.
PAN, N & BANERJEE, S 2022 Exact relations for energy transfer in simple and active binary fluid turbulence. *Phys. Rev. E* **106** (2), 025104.

PAN, N & BANERJEE, S 2025 Universal energy cascade in homogeneous binary fluid turbulence: A direct comparison of different exact relations. *Phys. Rev. Fluids* **10** (6), 064615.

PANDEY, V, MITRA, D & PERLEKAR, P 2023 Kolmogorov turbulence coexists with pseudo-turbulence in buoyancy-driven bubbly flows. *Phys. Rev. Lett.* **131** (11), 114002.

PANDEY, V, RAMADUGU, R & PERLEKAR, P 2020 Liquid velocity fluctuations and energy spectra in three-dimensional buoyancy-driven bubbly flows. *J. Fluid Mech.* **884**, R6.

PERLEKAR, P 2019 Kinetic energy spectra and flux in turbulent phase-separating symmetric binary-fluid mixtures. *Journal of Fluid Mechanics* **873**, 459–474.

PERLEKAR, P, BIFERALE, L, SBRAGAGLIA, M, SRIVASTAVA, S & TOSCHI, F 2012 Droplet size distribution in homogeneous isotropic turbulence. *Phys. Fluids* **24** (6).

RAMIREZ, G, BURLOT, A, ZAMANSKY, R, BOIS, G & RISSO, F 2024 Spectral analysis of dispersed multiphase flows in the presence of fluid interfaces. *International Journal of Multiphase Flow* **177**, 104860.

RICHARDSON, LF 1922 *Weather prediction by numerical process*. Cambridge University Press.

RISSO, F 2018 Agitation, mixing, and transfers induced by bubbles. *Annu. Rev. Fluid Mech.* **50**, 25–48.

RIVIERE, A, MOSTERT, W, PERRARD, S & DEIKE, L 2021 Sub-Hinze scale bubble production in turbulent bubble break-up. *J. Fluid Mech.* **917**, A40.

ROSALES, C & MENEVEAU, C 2005 Linear forcing in numerical simulations of isotropic turbulence: Physical space implementations and convergence properties. *Phys Fluids* **17** (9), 095106.

ROSTI, M E, GE, Z, JAIN, SS, DODD, MS & BRANDT, L 2019 Droplets in homogeneous shear turbulence. *J. Fluid Mech.* **876**, 962–984.

RUDMAN, M 1998 A volume-tracking method for incompressible multifluid flows with large density variations. *Int J Numer Methods Fluids* **28** (2), 357–378.

SABELNIKOV, V. A., LIPATNIKOV, A. N., NISHIKI, S. & HASEGAWA, T. 2019a Application of conditioned structure functions to exploring influence of premixed combustion on two-point turbulence statistics. *Proc. Combust. Inst.* **37** (2), 2433–2441.

SABELNIKOV, V. A., LIPATNIKOV, A. N., LAEVICH, NISHIKI, S. & HASEGAWA, T. 2019b Investigation of the influence of combustion-induced thermal expansion on two-point turbulence statistics using conditioned structure functions. *J. Fluid Mech.* **867**, 45–76.

SUSSMAN, M, SMITH, KM, HUSSAINI, MY, OHTA, M & ZHI-WEI, R 2007 A sharp interface method for incompressible two-phase flows. *J Comput Phys* **221** (2), 469–505.

TAGUELMIMT, N, DANAILA, L & HADJADJ, A 2016 Effects of viscosity variations in temporal mixing layer. *Flow Turbul. Combust.* **96**, 163–181.

TAYLOR, M A, KURIEN, S & EYINK, G L 2003 Recovering isotropic statistics in turbulence simulations: The Kolmogorov 4/5th law. *Phys. Rev. E* **68** (2), 026310.

TERRINGTON, S.J., HOURIGAN, K & THOMPSON, MC 2022 Vorticity generation and conservation on generalised interfaces in three-dimensional flows. *J. Fluid Mech.* **936**, A44.

THIESSET, F., DURET, B., MÉNARD, T., DUMOUCHEL, C., REVEILLON, J. & DEMOULIN, F-X. 2020 Liquid transport in scale space. *J. Fluid Mech.* **886**, A4.

THIESSET, F & FEDERRATH, C 2023 Structure of iso-density sets in supersonic isothermal turbulence. *Astronomy and Astrophysics-A&A* **676** (A12).

THIESSET, F, MÉNARD, T & DUMOUCHEL, C 2021 Space-scale-time dynamics of liquid-gas shear flow. *J. Fluid Mech.* **912**, A39.

THIESSET, F. & POUX, A. 2020 Numerical assessment of the two-point statistical equations in liquid/gas flows. *Tech. Rep.* CNRS, Normandy Univ., UNIROUEN, INSA Rouen, CORIA.

TORQUATO, S. 2002 *Random Heterogeneous Materials. Microstructure and Macroscopic Properties*. Springer-Verlag New York.

TRAUTNER, E, KLEIN, M, BRÄUER, F & HASSLBERGER, J 2021 Conditional and unconditional second-order structure functions in bubbly channel flows of power-law fluids. *Phys. Fluids* **33** (5).

TREFFTZ-POSADA, P & FERRANTE, A 2023 On the interaction of Taylor length-scale size droplets and homogeneous shear turbulence. *J. Fluid Mech.* **972**, A9.

VALENTE, PC, DA SILVA, CB & PINHO, FT 2014 The effect of viscoelasticity on the turbulent kinetic energy cascade. *Journal of fluid mechanics* **760**, 39–62.
VAUDOR, G, MÉNARD, T, ANISZEWSKI, W, DORING, M & BERLEMONT, A 2017 A consistent mass and momentum flux computation method for two phase flows. Application to atomization process. *Comput and Fluids* **152**, 204–216.

WAGNER, R, FALKOVICH, G, KRITSUK, A G & NORMAN, ML 2012 Flux correlations in supersonic isothermal turbulence. *J. Fluid Mech.* **713**, 482–490.

WANG, J, WAN, M, CHEN, S & CHEN, S 2018 Kinetic energy transfer in compressible isotropic turbulence. *J. Fluid Mech.* **841**, 581–613.

WHITMAN, S. H. R., TOWERY, C. A. Z., POLUDNENKO, A. Y. & HAMLINGTON, P. E. 2019 Scaling and collapse of conditional velocity structure functions in turbulent premixed flames. *Proc. Combust. Inst.* **37** (2), 2527–2535.

WILCOX, DC 1998 *Turbulence modeling for CFD*, vol. 2. DCW industries La Canada, CA.

YAO, H, ALVES-PORTELA, F & PAPADAKIS, G 2020 Evolution of conditionally averaged second-order structure functions in a transitional boundary layer. *Phys. Rev. Fluids* **5** (9), 093902.

YAO, H & PAPADAKIS, G 2023 On the role of the laminar/turbulent interface in energy transfer between scales in bypass transition. *J. Fluid Mech.* **960**, A24.

ZAMANSKY, R, DE BONNEVILLE, FLR & RISSO, F 2024 Turbulence induced by a swarm of rising bubbles from coarse-grained simulations. *J. Fluid Mech.* **984**, A68.

ZHANG, J 1996 Acceleration of five-point red-black Gauss-Seidel in multigrid for Poisson equation. *Appl. Math. Comput.* **80** (1), 73–93.