Modulation of homogeneous and isotropic turbulence in emulsions

We present a numerical study of emulsions in homogeneous and isotropic turbulence at $Re_\lambda=137$. The problem is addressed via Direct Numerical Simulations (DNS), where the Volume of Fluid (VOF) is used to represent the complex features of the liquid-liquid interface. We consider a mixture of two iso-density fluids, where fluid properties are varied with the goal of understanding their role in turbulence modulation, in particular the volume fraction ($0.03<\alpha<0.5$), viscosity ratio ($0.01<\mu_d/\mu_c<100$) and large scale Weber number ($10.6<We_\mathcal{L}<106.5$). The analysis, performed by studying integral quantities and spectral scale-by-scale analysis, reveals that energy is consistently transported from large to small scales by the interface, and no inverse cascade is observed. Furthermore, the total surface is found to be directly proportional to the amount of energy transported, while viscosity and surface tension alter the dynamic that regulates energy transport. We also observe the $-10/3$ and $-3/2$ scaling on droplet size distributions, suggesting that the dimensional arguments which led to their derivation are verified in HIT conditions.


Introduction
Emulsions are multiphase flows of two immiscible (totally or partially) liquid phases with similar densities. Such flows are extremely common in industrial applications such as pharmaceutical (Nielloud 2000;Spernath & Aserin 2006), food processing (McClements 2015), oil production (Kokal & Others 2005;Mandal et al. 2010;Kilpatrick 2012) and waste treatment. Emulsions are also relevant for environmental flows such as oil spilling in oceans, when the oil droplets distribution becomes fundamental for quantifying environmental damages (Li & Garrett 1998;French-McCay 2004;Gopalan & Katz 2010). Many studies have been performed on the rheological behavior of emulsions in the past (Einstein 1906(Einstein , 1911Pal 2000Pal , 2001Jansen et al. 2001;De Vita et al. 2019), while the current knowledge on their behavior in turbulent flows is limited (Yi et al. 2021).
The two fluids are usually referred to as continuous phase (or carrier phase in case of strong advection) and dispersed phase (or droplet-phase) depending on whether the volume fraction α is respectively greater or lower than 0.5; the system is denoted as binary flow when α = 0.5. As the density ratio is usually considered to be close to 1, gravity effects are negligible with respect to the stirring and advection needed to sustain turbulence in the flow. For this reason, four dimensionless numbers can be used to describe these flows, namely the Reynolds number Re, the Weber number W e, the volume fraction of the dispersed phase and the viscosity contrast. Depending on the specific configuration under investigation, the definition of these numbers can change, yet they completely define the case studied provided the two fluid have the same density. Several aspects of fundamental importance in emulsions, such as turbulence modulation, droplet size distributions and inter-phase energy fluxes, are not fully understood. We therefore aim to partially fill this gap by means of numerical simulations. In the following we provide an overview of the main results available in literature. Results for bubble/droplet laden flows are also discussed when relevant to the present work.

Observations on droplet size distribution
The Droplet Size Distribution (DSD) is a key aspect of emulsions, as its prediction becomes fundamental in most applications. In his early seminal work, Kolmogorov (1949) discussed the criteria under which a droplet undergoes breakup when subject to surrounding turbulence. Kolmogorov first proposed a dimensional argument according to which surface tension forces need to be locally balanced by turbulent energy fluctuations. This idea was later addressed in Hinze (1955) and translated into a critical Weber number W e c of order 1 at which breakup occurs, leading to the definition of the Hinze scale d H as the minimum droplet diameter at which breakup may occur due to pressure fluctuations. A general definition for this scale is: where σ is the surface tension coefficient, ρ c is the carrier phase density and ε is the energy dissipation rate. This estimate proved valid for bubbles (Masuk et al. 2021;Chan et al. 2021) and emulsions (Perlekar et al. 2012;Mukherjee et al. 2019;Rosti et al. 2020;Yi et al. 2021). Different O(1) values have been reported for W e c in numerical (Rivière et al. 2021) and experimental works (Deane & Stokes 2002;Lemenand et al. 2017), from 0.5 up to 5; for dilute emulsions in turbulence W e c ≈ 1.17, according to the values from both numerical (Perlekar et al. 2012) and experimental (Yi et al. 2021) data. For bubbles larger than the Hinze scale, Garrett et al. (2000) found that, in isotropic turbulent conditions, droplets break with a cascade process, and the diameter distribution follows a d −10/3 power-law. This deterministic process can accurately describe bubble size distributions in breaking waves obtained in experiments (Garrett et al. 2000;Deane & Stokes 2002;Qi et al. 2020) and numerical simulations (Deike et al. 2016;Chan et al. 2021). The same power law has also been proposed for emulsions, based on diffuseinterface numerical simulations (Skartlien et al. 2013;Mukherjee et al. 2019;Soligo et al. 2019). For bubbles smaller than the Hinze scale, Deane & Stokes (2002) suggested the existence of a fragmentation process; in this case, a d −3/2 power-law is used to accurately fit experimental data. Agreement with this empirical power-law has been observed in Homogeneous and Isotropic Turbulence (HIT) both for bubbles (Rivière et al. 2021) and emulsions (Mukherjee et al. 2019). The transition between the two power-laws is defined by the Hinze scale. A consequence of this transition is that droplets with d d H generate both local and non-local bubble/droplet production, as they can fragment in both droplets larger or smaller than the Hinze scale (Rivière et al. 2021). Although both power-laws have been derived under the hypothesis of dilute conditions (α 0.05) they have been recently observed in HIT studies of dense emulsions (Mukherjee et al. 2019), rising the question on the effective role of coalescence in the process.
The connection between bubbles and emulsions is non-trivial and deserves special attention. In his work, Hinze (1955) discussed how W e c depends on the fluid properties of the dispersed phase. He assumed that W e c = C[1 − f (N V i )], with f a generic function of the viscosity group N V i = µ d / √ ρ d σd, where µ d is the dispersed phase viscosity. On the other hand, d H was derived under the assumption of a dilute emulsion, hence the density in Equation (1.1) refers to the carrier phase, as the phase where the energy dissipation rate ε could be measured in experiments. This allows the direct application of the Hinze criteria in flows where density/viscosity ratios are significant as in air-water flows. However, significant uncertainties are discussed in literature about the properties of the function f and the role of the dispersed phase properties remains mostly unknown (Masuk et al. 2021). Also unknown is the role of turbulence inhomogeneity and anisotropy, which, according to Hinze (1955), may be a further source of non-linear effects in the determination of W e c . In fact, in flows where the energy dissipation rate shows strong spatial variations, W e c varies for each bubble/droplet and it assumes meaning only on an average sense, making it difficult to disentangle the effects of turbulence anisotropy and property contrast. Despite all these uncertainties, correlations from Hinze (1955); Garrett et al. (2000); Deane & Stokes (2002), derived for isotropic turbulent conditions, applies in most studies with strong property contrasts and large-scale anisotropy. This is likely due to the underlying assumption that the breakup process is purely inertial, as it only depends on ε (Garrett et al. 2000). Thus, bubble breakup studies become relevant also for the present study.
It is finally worth noticing that the flow configuration appears to have a significant impact on DSD and experimental observations in shear flows can depart quite substantially from the discussed power-law behaviors. The recent work of Yi et al. (2021) presents strong experimental evidences of gamma/log-normal DSD in Taylor-Couette flow, confirming the previous findings of Pacek et al. (1998). These configurations are characterized by strong anisotropy, making the comparison with data obtained for emulsions and bubbles in HIT difficult. On the other hand, Soligo et al. (2019) studied breakup and coalescence of emulsions dynamic in a turbulent channel flow. These authors observed the appearance of the −10/3 power-law for the droplet size distribution in presence of surfactants. It is interesting to observe that, in this numerical study, the scaling from Garrett et al. (2000) seems to apply in anisotropic configurations. Fortunately, there has been a significant effort in recreating local HIT conditions in experiments in the latest years (Debue et al. 2018;Dubrulle 2019;Knutsen et al. 2020) and new studies are expected to provide new insights on these aspects.

Studies of two-fluid turbulence
With the advent of more powerful computational resources, a significant number of studies have considered droplets in turbulent flows, yet almost only through diffuseinterface methods which may display significant mass loss. In their study of emulsions in HIT turbulence, Perlekar et al. (2012) show that a statistical stationary state can be reached for the droplet size distributions. In the study, the authors used the Pseudo-Potential Lattice Boltzmann method (Biferale et al. 2011), which compensates mass losses (due to droplets dissolution) by artificially re-inflating existing ones. Simulations of the Cahn-Hilliard-Navier-Stokes formulation are presented in Perlekar et al. (2014) for binary fluids. These authors found that enforcing large-scale HIT arrests coarsening. This result is particularly significant for emulsions (of which binary fluids represent a special case) as it shows that turbulence is the main factor to determine the droplets size. Furthermore, these authors report modified energy spectra for the mixtures, with a crossover in correspondence to the Hinze scale. Komrakova et al. (2015) used a free-energy lattice Boltzmann method to numerically simulate emulsion breakup in HIT, induced by an external large-scale linear forcing. Their findings show that energy spectra present deviations with respect to the singlephase configuration and that the numerical method employed may alter the small-scale dynamics of the flow. Finally, increased coalescence is found for volume fractions α > 0.05 also owing to the nature of the diffused interface method.
Droplet interactions with turbulence have been studied by Dodd & Ferrante (2016) in decaying isotropic turbulence. Amongst several observations, these authors discuss the effects of droplet breakup and coalescence on the turbulent kinetic energy budget. Droplet coalescence lowers the total amount of area, hence decreases the surface energy and consequently increases the kinetic energy locally, while the opposite occurs in the case of breakup. More recently, Mukherjee et al. (2019) have studied emulsions in HIT conditions using a pseudo-potential Lattice Boltzmann method, discussing droplet statistics and their correlation with the surrounding turbulence. They confirm the findings of Perlekar et al. (2014) for energy spectra pivoting at the Hinze scale, demonstrating that energy is subtracted from large scales and injected at small scales, while no direct observation of the underlying mechanism is presented. These authors also show that the droplet generation can be described through the Weber number spectra. In the same work, Mukherjee and co-workers discuss and demonstrate the need of using a forcing scale smaller than the turbulent-box size in order to achieve a polydisperse droplet distribution. It is yet important to note that Mukherjee et al. (2019) used a pseudopotential lattice-Boltzmann method, which leads to a significant loss of the dispersed mass during the simulation, as fairly discussed by the authors.
As concerns binary fluids, Perlekar (2019) shows how the presence of interfaces leads to a different energy transfer mechanism, confirming the conclusions in Dodd & Ferrante (2016). The author uses the scale-by-scale (SBS) energy balance to show that the energy absorption at larger scales is mainly given by the interface source term in the Cahn-Hilliard equation used by the author to describe the multiphase nature of the flow. Furthermore, Perlekar (2019) shows that small-scale statistics are almost unchanged when changing W e while they are affected by the Reynolds number. This study complements the previous findings in binary fluids (Perlekar et al. 2014(Perlekar et al. , 2017, where coarsening was analyzed in 3D and 2D turbulence by means of a spinoidal decomposition. Rosti et al. (2020) study droplets in Homogeneous Shear Turbulence (HST), focusing on the effect of the droplet initial diameter and the shear-rate magnitude; the results show that a statistically stationary regime (i.e. balance of coalescence and breakup events and energy balance convergence) can be reached, while the Taylor-scale Reynolds number Re λ decreases with increasing surface tension.
Despite the growing literature on the subject, many issues remain to be fully understood. In particular, most of the studies have been carried out using diffuse-interface approaches, which cannot exactly represent the surface-terms effects yet key in many occasions. In this sense, our work complements the very recent one by Rivière et al. (2021) focused on the bubble break-up dynamics.

Objectives of the present study
In the present work, we use Direct Numerical Simulations (DNS) to study the effects of viscosity ratio, volume fraction and surface tension on the emulsion turbulent behavior. The chosen setup is tri-periodic HIT, with turbulence sustained throughout the simulation time. The analysis is performed at Re λ ≈ 137, large enough to represent realistic turbulent flows, while volume fraction, viscosity ratio and surface tension are varied to cover most relevant applications (Jansen et al. 2001). We analyze the turbulence through global and phase-averaged energy balance, energy spectra, SBS energy budget, and Probability Density Functions (PDF) for the intermittency analysis. Furthermore, we discuss droplet size distributions for all cases. In summary, we will show that (i) the energy balance is significantly altered by the properties of the dispersed phase; (ii) surface tension forces induce an additional mechanism for energy transfer from larger scale towards the energy dissipation range; (iii) the modified energy transport mechanism alters the energy spectra; (iv) the presence of the interface increases intermittency and alters the small scale statistics; (v) the droplet size distribution displays both the d −3/2 and d −10/3 power-laws, with remarkable accuracy also for d < d H .

Governing equations and numerical method
We consider an incompressible flow obeying the continuity and Navier-Stokes equations: where u i is the velocity in the i-th direction, p is the pressure, ρ and µ the local density and viscosity. The forcing term f σ i = σξδ s n i represents the surface tension force, where σ is the surface tension, ξ the local interface curvature, n i the i-th component of the surface normal vector and δ s the Dirac delta function that ensures the surface force is applied only at the interface (Tryggvason et al. 2011). The last term in equation Equation (2.1b) is the forcing needed to sustain turbulence by injecting energy at the large scales; among the several algorithms available to force sustained homogeneous and isotropic turbulence (e.g. Eswaran & Pope 1988;Rosales & Meneveau 2005;Mallouppas et al. 2013;Bassenne et al. 2016), we use here the Arnold-Beltrami-Childress (ABC) forcing (Mininni et al. 2006), (2.2) with x, y and z ∈ [0, 2π]. As reported by Podvigina & Pouquet (1994), the ABC forcing creates an unstable single-phase flow for 1/ν > 20, with ν the kinematic viscosity and κ 0 the forcing wavelength.
The description of the code and the algorithm used can be found in Rosti et al. (2019Rosti et al. ( , 2020, together with several validations. The method is therefore only shortly described here, see also Costa (2018) for references to the code structure.
The equations are discretized on a staggered uniform Cartesian mesh: the spatial derivatives are computed using second-order centered finite differences and a second-order Adam-Bashford scheme is used for the time integration. The pressure splitting method presented in Dodd & Ferrante (2014) is used to obtain a constant-coefficient Poisson equation, which we then solve with the direct FFT-based pressure solver presented in Costa (2018).
The interface between the two fluids is described with the Volume of Fluid (VOF) method, in particular the Multi-dimensional Tangent Hyperbola INterface Capturing (MTHINC) algorithm developed by Ii et al. (2012). The advection equation for the VOF can be written in divergence form as where H is the color function assuming the value of 0 and 1 in each of the fluids, and φ the cell-averaged value of H. In the MTHINC method, the function H is locally approximated using the hyperbolic tangent, where (x , y , z ) ∈ [0, 1] is the cell-cented local coordinate system, β is a sharpness parameter (equal to 1 in the current work), d a normalization factor and P is the threedimensional surface function, assumed here to be quadratic (Ii et al. 2012). The advantage of the method is that Equation (2.4) allows to solve the fluxes in Equation (2.3) by semianalytical integration. Once the VOF function φ is known, we evaluate the local fluid properties as where the subscripts c and d indicates carrier and dispersed phase. Finally, the Continuum Surface Force (CSF) model is used to compute the surface tension force (Brackbill et al. 1992), with the normal evaluated with Youngs' method and the curvature as in Ii et al. (2012).

Flow configuration
All the simulations are performed using the same ABC forcing, injecting energy at wavenumber κ 0 = 2π/L = 2, with A = B = C = 1, corresponding to Re λ ≈ 137 for the single phase flow (see Section 2.4 for the characteristics of the reference single-phase flow and definition of the meaningful observables). As reported in literature (Komrakova et al. 2015;Mukherjee et al. 2019) forcing the second wavelength is recommended in order to avoid coalescence induced by large turbulent structures in periodic domains.
In addition to the Reynolds number, the emulsion flows are characterised by 4 nondimensional parameters. The volume fraction, α = V d /V, defined as the ratio between the volume occupied by the dispersed phase V d and the total volume V = (2π) 3 , the viscosity ratio µ d /µ c , where the subscripts d and c indicate the dispersed and carrier phase, and the Weber number, W e L = ρ c Lu 2 rms /σ, where u rms is the space-time average of the root-mean-square velocity of the single-phase case (which can be related to the forcing amplitude A = B = C) and L the scale of the ABC forcing. Finally, the density ratio, ρ = ρ c /ρ d , is kept constant equal to 1 in this study.
Here, we will vary the dispersed phase volume fraction, the viscosity ratio and the Weber number; the parameters pertaining the different simulations discussed below are presented in Table 1. Note, finally, that the table also indicates the integration time N T required to reach statistical convergence of the turbulent quantities and droplet-sizedistribution (DSD) in units of large eddies turnover times, T = Lu rms (Mininni et al. 2006). The simulations are considered at convergence when global energy production balances dissipation (see Section 2.3 and Equation (2.10) for their definition) with an error of less than 4%, also implying that the area derivative over time is negligible (see Section 2.3 for further details). Interestingly, N T varies significantly with the physical configuration. In particular, starting with the reference cases BE1 and BE2, we observe that increasing µ d /µ c longer times are needed to reach a statistically stationary state, which we will attribute to a decrease of the breakup rate. A similar behavior is observed when decreasing W e L , when higher surface tension forces decrease the probability of breakup events. Finally, large structures become unavoidable when increasing the volume fraction α (Komrakova et al. 2015;Mukherjee et al. 2019), which implies longer simulation times.
Visualisations of the transient phase to reach the final steady state are reported in Figure 1 for the reference case BE1 with α = 0.03. The simulation starts at t 0 using the fully developed single-phase HIT field from case SP2. The dispersed phase is initialised as an ensemble of spheres, which soon deform in the flow as shown in panel b) pertaining time t 1 = T /4. At statistical convergence, t ≈ 10T , when statistics are collected, we observe a poly-dispersed distribution of asymmetric droplets. Note finally that for α 10% the simulations are initialised using spherical droplets of size d 0 ≈ 0.12L, while a single spherical droplet of initial size d 0 = (6αL 3 /π) 1/3 was used for larger values of α. We have checked that the initial distribution has no effect on the final droplet size distribution, as also reported in Mukherjee et al.  . The droplet are initialized at t0 in a developed turbulent field. As turbulence is maintained, breakup and coalescence start occurring (t1), and statistical convergence in the DSD is achieved after a few turnover times (t2)

Observables, phase-averaged energy balance and scale-by-scale budget
In this section, we discuss the theoretical tools and the physical observable that will be discussed throughout the study. The objective of this study is to understand the turbulence modulations induced by a second phase, focusing on comparisons of the energy spectra and the SBS analysis. In particular, we will consider the Taylor scale Reynolds number, the energy spectra, the phase-averaged energy balance, Probability Density Functions (PDF) of velocity fluctuations and vorticity, and the SBS budgets. The Taylor scale Reynolds number is defined as where λ = (5νu i u i /ε) 1/2 is the Taylor scale, with the energy dissipation rate computed as ε = ν∂ i u j ∂ j u i . For the reference single-phase flow Re λ = 137. Here, we compute ε and all the relevant observables at each computational grid point and then average in space and time. This procedure is required due to material properties discontinuities when µ d /µ c is varied. Note that, from now on, ε will denote the space-time averaged value. For the a multiphase flow, the energy balance is obtained by multiplying Equation (2.1b) by the velocity u i We define the volume average as where the subscript m represents an integral over the dispersed phase d, the carrier phase c or the total volume, if omitted. Applying the operator · to Equation (2.7) leads to where k is the turbulent kinetic energy. Due to the homogeneity of the HIT configuration, the transport term arising from the nonlinear transport in Equation (2.7) vanishes. Further details on its derivation for the case of emulsions can be found in Dodd & Ferrante (2016); Rosti et al. (2020). It can be proven that Ψ σ ∝ ∂ t A (Dodd & Ferrante 2016) (with A the total interface area) and that the contribution of the surface tension to the total energy variation also vanishes, since the time derivative is zero at the statistically stationary state. Hence, we obtain that the production term P is perfectly balanced by the energy dissipation ε.
Next, we consider the phase-averaged energy balance, following the approach described in Dodd & Ferrante (2016); Rosti et al. (2020). Averaging Equation (2.7) on each phase (i.e. enforcing eq. 2.8), we obtain the phase-average energy balance: Here, P m and ε m indicate production rate and viscous dissipation rate per unit volume in each phase. The terms T ν m and T p m are the viscous and pressure transport densities and represent the coupling between the two phases; when the sum of these two, T m = T ν m +T p m is positive, energy is absorbed from phase m, when negative energy is transferred to the other phase. Again, in statistical stationary conditions, ∂ t k m ≈ 0. We now move to spectral space and present the SBS balance. This is derived for the two-fluid flows following the formulation in Olivieri et al. (2020a,b); for more details the reader is refereed to Frisch (1995); Alexakis & Biferale (2018). Taking the Fourier transform of the momentum equations (eq. 2.1b), we obtain Note that as the viscosity µ is a function of space and time, we actually compute the dissipation term in physical space to avoid a convolution in the spectral space. We next multiply Equation (2.11) with the complex conjugate of the velocityũ * i and drop the pressure term by imposing the incompressibility condition κ iũi = 0 as in this work ρ c = ρ d = 1. Multiplying the complex conjugate of Equation (2.11) bỹ u i , summing the equations obtained forũ andũ * and averaging in time, we finally obtain i t is the time-averaged kinetic energy in the spectral domain, whose time derivative is zero at statistical steady state; iũ i t is the time-averaged energy transfer due to the non-linear term; iũi t is the time-averaged work of the surface tension force at the different scales; iũ i t is the time-averaged energy input due to the large-scale forcing. All of the above are three-dimensional fields in spectral space. Note that, at steady state when the total interfacial area is constant, S σ integrates to zero (Dodd & Ferrante 2016) so that this term can be effectively seen as an energy transport due to the surface tension. Finally, we perform a spherical-shell integral in spectral space and express each term in the budget as a function of the magnitude of the wavevector κ. This operation results in, e.g., (2.13) This term represents the shell-to-shell energy transfer function for the non-linear term of the momentum equation and similarly for the other terms above. If we instead perform the integration over a sphere (i.e. for all |κ i | < κ) we obtain the cumulated SBS budget: (2.14) In the expression above, the fluxes Π(κ) = |κi|<κ T (κ i ) and Π σ (κ) = |κi|<κ S σ (κ i ), indicate the energy flux from all the largest scale to κ i , and are typically used to study scalings in the inertial range (where Π(κ) = ε) and the direction of the energy cascade (Alexakis & Biferale 2018). The remaining terms represent the energy injected and the dissipation at all scales below κ i . The cumulative SBS budget can easily be related to the energy balance in the physical domain: Equation (2.14) and Equation (2.9) are equivalent for κ = κ max , hence it can be easily demonstrated that Π(κ max ) = Π σ (κ max ) = 0 and ε = D(κ max ) = P = F(κ max ). In this work, we will mostly show the shell-by-shell energy budget (eq. 2.12, integrated using 2.13, referred to SBS if not differently specified), as more suited for detailed comparisons at each scale, while the cumulative energy budget is used for the single phase flow only.

Analysis of the reference single-phase flow and grid convergence
We will now motivate the choice of the resolution adopted for the emulsion simulations, i.e. N = 512 points in each direction. To this aim, we will first analyze the behavior observed in single phase turbulence. The energy spectra and the cumulative SBS balance pertaining the single-phase flow are shown in Figure 2. A good agreement between the cases SP1 and SP2 is evident in panel (a). The κ −5/3 law for the inertial range extends over almost a decade, showing a fully developed turbulent flow at a moderately high Reynolds number. The 2 cases yield the same result in terms of ε, with a relative error of less than 5%. The cumulative SBS balance shows that, due to the moderate Re λ , the viscous term D is not negligible already at large scales, where it dissipates approximately 3% of the injected energy. This is in agreement with the observation that a fully developed inertial range is only partially observable for Re λ 200 in single-phase turbulence (Ishihara et al. 2009). A substantial dissipative range is present for κ 10 2 , indicating an accurate computation of the smallest scales. In this region we can observe that Π(κ max ) = 0 and |κi|<κ D(κ i ) = ε. As a consequence of the imposed ABC forcing, energy injection is clearly observed for κ = 2. The SBS budget shows almost no difference between the results from the two grid resolutions, and that all relevant processes are already accurately captured at the lower resolution, N = 256, down to the smallest scales.
To investigate convergence of the multiphase flows, we consider the case W e L = 42.6, µ d /µ c = 1 and α = 0.5 for three different resolutions, from N = 256 to N = 1024 corresponding to cases C14, C24, and C34, see Figure 3. This configuration was selected because it is the one with the largest interfacial area and largest fluctuations, dA/dt, and hence where larger errors in the averaged energy budget are expected. Nevertheless the spectra in panel (a) do not seem to be significantly affected by the grid resolution and the dissipative range is observed also at the lowest resolution, N = 256.
A more stringent test is the convergence of the SBS budget, here (and hereafter) shown in its shell-by-shell form, see panel (b) where all terms are normalized by ε and premultiplied by the wavelength κ to improve readability. Comparing the data at different resolution, the energy injection at large scales F and the energy transfer by the nonlinear term T are almost identical in the inertial range. The energy dissipation D and the transfer due to interfacial forces S σ display some differences for κ 10. If we consider the integral contribution from the surface tension, which should theoretically be zero, Π σ (κ max )/ε ≈ 0.08 for the lowest resolution, N = 256; this value decreases for N = 512 and remains almost constant at N = 1024, Π σ (κ max )/ε ≈ 0.04. It is worth underlying that this is the largest error encountered among all cases, since Π σ (κ max )/ε ≈ 0.01 for most of the other cases discussed in this study. Overall, the energy not resolved by the the simulation with a grid of N = 512 and the differences in the SBS transfer functions can be considered as negligible for the scope of the present study, where we wish to primarily examine the energy transfer towards the smallest scales.
The data in figure 3 highlight already important features of emulsions in HIT, which will be consistently observed in all cases studied. First, the energy at the smallest scales increases with respect to the single-phase case (panel a). Secondly, the presence of the interface alters the behavior of the turbulent flow, with the surface tension forces S σ transferring energy from large scales towards the dissipative range. Because of this, the total energy transported by the non-linear term T reduces and the dissipative range extends towards smaller scales, where we observe a balance between the work of surface tension and viscous dissipation. These modified flow features are similar to what found by Olivieri et al. (2020a,b) for fiber suspensions in turbulent flows.

Emulsions at different volume fractions
We first examine the influence of the dispersed-phase volume fraction on the turbulent flow, cases BE1 and Cxx in Table 1, corresponding to increasing values of α from 3 to 50%. A render of the cases discussed here is shown in Figure 4, where the isocontour of VOF fields are shown for volume fractions 0.06, 0.2 and 0.5.
The modulation of the turbulence is first quantified in terms of integral quantities. Figure 5 shows Re λ , computed according to Equation (2.6), versus the volume fraction α. Re λ increases almost linearly with α, by approximately 15% for α = 0.5. A similar trend is found for λ, as shown in the inset of Figure 5. Considering that the average of ε and k is approximately constant in all cases (variations of ±3%), the increase of Re λ and λ is therefore due to the local variations of the ratio k/ε. In particular, the increased values (k/ ) for similar averaged values of the two quantities is attributed to the increased correlation between regions of strong turbulent kinetic energy and low dissipation. A graphical evidence is presented in Figure 6, where we show the instantaneous ratio k/ for the single-phase (case SP2 in panel a) and multiphase flow (case BE2 with α = 0.1 in panel b) in logarithmic scale. The figure shows that when the dispersed phase is present, large regions of fluid with higher k/ are observed far from the droplet interface (denoted with white line). This can be explained as follows: as the total dissipation is constant, the local increase of near the interface, as also observed in Dodd & Ferrante (2016), correspond to a decrease of the dissipation rate in large portions of the fluid, those far from an interface. Considering that the turbulent kinetic energy is less affected by the presence of the interface, the ratio (k/ ) increases in average. To support this statement, Figure 6 for the same planes. In the emulsion (panel d ), higher values of ε are found close to the droplet interface and to the clustering regions, while for the single-phase flow (panel c) no specific pattern is observed. The one-dimensional energy spectra E(κ) multiplied by κ 5/3 , i.e. the so-called compensated spectra, are displayed in Figure 7 for the different α considered. The Taylor scale of the single-phase flow is indicated by the dot-dashed line, while the vertical dotted black line is used for the wavenumber κ H corresponding to the Hinze scale, defined as: Note that, the prefactor 0.725 in Equation (3.1) is set accordingly to the original work of Hinze (1955) for emulsions in HIT conditions, corresponding to W e c = 1.17. The data in Figure 7 reveal that the presence of the dispersed phase reduces the energy with respect to the single-phase case (SP1) for κ < κ H . At the same time, the energy content increases at the smaller scales, κ > κ H , in the dissipative range of the spectra. As noted in previous studies (Mukherjee et al. 2019), the amount of energy subtracted to the large scales is proportional to the volume fraction α. Interestingly, the wavenumber at which the curves cross over from reduced to increased energy content corresponds to the Hinze scale. For brevity, we will denote as pivoting point the wavelength where the spectra of the multiphase cases intersect the one from the single-phase reference case.
Pivoting points were not clearly observed in some previous studies on emulsions (Mukherjee et al. 2019;Perlekar 2019), while they are clearly visible in others (Perlekar et al. 2014;Dodd & Ferrante 2016;Rosti et al. 2020); this is possibly due to the different methods used to simulate the dispersed phase: the ability of the VOF to accurately resolve the interface reduces the energy dissipation by the surface tension term in the dissipative range. Such energy dissipation is indeed clearly observed by Perlekar (2019) who present results obtained by solving the Cahn-Hilliard equation in a diffuse-interface formulation. As mentioned in Section 2.4, these numerical artefacts do not have significant effects on the dynamics at the inertial range, while they affect the dissipative range. This aspect will be discussed also later in this section.   Insight on the energy transfer among the different scales is gained be using the SBS analysis. The full SBS energy budegt, i.e. the contributions from the different terms in Equation (2.12), is displayed in Figure 8(a) for case BE2, chosen as an illustrative example with an intermediate value of α = 0.1. The external forcing is injecting energy at κ = 2, which is absorbed by the non-linear transfer term T , for a large majority, and by the surface tension term S σ , for a small part. The non-linear term transfers energy towards smaller scales, larger values of κ. The surface tension term, S σ , acts as a dissipative process at large scales, where it absorbs approximately the same energy as the dissipative term D; however for 10 < κ < 20, we observe a significant change in the energy transport mechanism: S σ becomes positive, hence contributing to transferring energy towards the small scales, similarly to T , a process active until κ max . It is important to note that the surface tension transport remains active also at small scales in the dissipative range, consequently extending the range of wavelengths where the dissipation term remains active. These observations confirm the previous findings obtained in Perlekar (2019) for binary mixtures.
The details of the effect of the volume fraction α on each term of the SBS balance are displayed in Figure 8(b-d ). We first analyze the non-linear transfer term T in panel (b). As α increases, T absorbs progressively less energy at the injection frequency κ = 2.
Consequently, less energy is transferred towards smaller scales by nonlinear advection. The energy flux Π (not shown) do not display an inverse cascade for any α. Furthermore, we notice that no energy is transferred to the far end of the dissipative range, which is resolved over a large range of scales in all cases.
The contribution from the surface tension S σ , see Figure 8(c), confirms that interfacial stresses absorb part of the energy injected into the domain at κ = 2. The energy absorbed by the surface tension term at large scales is approximately proportional to α. The surface tension term becomes positive at smaller scales, where energy is released. The positive peak is reached at approximately the Hinze scale for all cases. As for the energy absorption, the magnitude of the peak scales proportionally to α. We also notice that, for any α, the surface tension terms act also in the dissipative range, where the non-linear term T is zero.
The behavior observed so far for T and S σ provide a clear explanation for the previous observations on the energy spectra. At small wavenumber, the energy cascade produced by the non-linear energy transfer is partially inhibited by the presence of the interfacial forces. For high wavenumbers, T progressively reaches zero, but the energy previously subtracted by the interfacial stresses at large scale is redistributed at small scales, which can be seen in Figure 7 as an energy increase at high wavenumbers.
To close the SBS balance, we examine the viscous dissipation term D, see figure Figure 8(d ). First, we note that only a small amount of the injected energy (less than 5% for all cases) is absorbed by the dissipation term at the scale of the forcing, κ = 2. The overall effect of the dispersed phase is to shift the energy dissipation towards smaller scales. This constitutes the natural reaction of the system to the increased activity in the dissipative range caused by the surface tension term. This behavior becomes more evident as α increases and progressively enhances dissipation at those small scales where the single-phase dissipation is negligible.
Summarising, the surface tension introduces an alternative path for energy transmission from large towards small scales, as discussed for binary flows in Perlekar (2019). The amount of energy transferred by the surface tension is directly proportional to the total droplets surface area A, as shown in Figure 9(a) where we display the maximum energy transferred via surface tension max ( κmax i=1 S σ (i)) and the total area of the dispersed phase A for the different volume fractions under consideration with a linear fit to the data. This observation reinforces our previous conclusion that the interface transfers energy among different scales by disrupting larger turbulent structures and creating smaller ones, hence affecting the canonical -5/3 slope of the turbulence spectra. Note also that, while in mono-dispersed flows this results in a deviation at a specific spectral frequency (Dodd & Ferrante 2016), in poly-dispersed flows this behaviour is seen at all scales.
We next consider the dynamics of the dispersed phase. We first examine the DSD for all the values of volume fraction studied, see Figure 9(b), where we display the droplet diameters normalized by the single-phase (SP1) Kolmogorov scale η sp . The dashed black line depicts the d −3/2 law by Deane & Stokes (2002) and the solid line the d −10/3 law by Garrett et al. (2000) valid for larger droplets. For small droplets, the -3/2 law is well captured also for marginally resolved droplets (with d/η sp < 6). For droplets larger than the Hinze scale, the -10/3 law is also a very good fit, with increasing accuracy for increasing values of α. Our data are in agreement with the findings by Mukherjee et al. (2019) and explained by higher coalescence probability at higher volume fractions, leading to a bigger population of larger droplets. Interestingly, the Hinze scale turns out to approximately define the transition between the -3/2 and the -10/3 scalings as proposed in Deane & Stokes (2002), although for higher values of α, the onset of the d −10/3 power-law occurs at larger diameters. As the droplet distributions can be, to a good approximation, represented by these 2 laws, it follows that A ∝ α, explaining why A, S σ and α are linearly correlated (see panel a of the same figure).
We now consider the phase-averaged energy budget, introduced in Section 2.4. The different terms of Equation (2.10), production, dissipation and transport by pressure and viscous forces, are shown in Figure 10, normalized by the single-phase dissipation. We first observe that the total production and dissipation Single phase Figure 11: PDF of velocity fluctuations u, vorticity ω and dissipation . All quantities are normalized as standard score.
for α < 0.5. The energy production density P m (green markers in panel a), is higher in the dispersed phase for low volume fractions, while it is comparable to that of the carrier phase for α > 0.1. The energy dissipation rate per unit volume in the dispersed phase ε d (green markers in panel b) is also larger at low volume fractions and monotonically decreases with increasing α. The dissipation in the carrier phase, ε c , also decreases, as it compensates for the energy transport T m from the carrier flow towards the dispersed phase. The viscous transport (see blue markers in panel b of the same figure) is significantly lower than its pressure-induced counterpart, T p m , although they exhibit a similar behavior: they first increase until α = 0.1 and then decrease to reach zero for a binary mixture. Note again, that the total transport term is zero, i.e. the sum of T p and T ν from both phases. The case α = 0.5 deserves a specific mention. In this case, production and dissipation in the two phases are equal, hence the transport term T m = 0. Intuitively, it is not possible to define unambiguously a carrier and dispersed phase in binary mixtures; while the energy is locally transported from one phase to the other, the global average is zero for both pressure and viscous transfer.
We finally analyze the PDF of velocity fluctuations u n = u/σ u , vorticity fluctuations ω n = (|ω i | − |ω i | )/σ ω and energy dissipation ε n = (ε − ε )/σ ε , normalized by their standard deviation. In Figure 11(a), we observe that, while the PDF remains symmetric, the tails of the PDF of the velocity fluctuation strongly deviate from the typical pseudo-Gaussian behavior of single-phase turbulence (Sreenivasan & Antonia 1997;Jimenez 2000;Ishihara et al. 2009). As concerns the vorticity in panel (b), no deviation is observed in the Gaussian core (as defined in Sreenivasan & Antonia 1997). However, the distributions of the multiphase flows strongly depart from the single-phase case in the tails. In particular, the exponentially decaying tails have a higher exponent in the case of emulsions, indicating more events with strong vorticity. Interestingly, while increasing the volume fraction does not influence the value of the exponent, increasing α induces deviations in the distributions already at lower values of ω n . We observe a similar behavior for ε n in panel (c): the intermittency of the single-phase flow is amplified by the presence of the interface. As for the vorticity, departures from the single-phase distributions are observed at lower values of ε when increasing the volume fraction α. As a final general remark, the analysis of the PDFs reveals that strong deviations are induced by the presence of the interface, already at low volume fractions, overall increasing the intermittent behavior of the flow. As no collapse is observed for the normalized variables, it can be inferred that the small-scale statistics are affected by the presence of the interface.

Influence of viscosity ratio
We consider now the influence of the viscosity ratio on the flow turbulence, i.e. cases BEx, V1x and V2x in Table 1. The viscosity ratios analyzed span the range 0.01 < µ d /µ c < 100, while W e L = 42.7 for all cases. Two values of the volume fractions are considered, α = 0.03 (series V1x) and α = 0.1 (series V2x). A render of the two-fluid interface (corresponding to the value of the VOF function φ = 0.5) is shown in Figure 12 for cases V11, BE1 and V14 (from left to right). As µ d /µ c increases, larger droplets appear; at low viscosity ratios we find a significantly higher number of droplets.
We start by examining the Taylor Reynolds number of the emulsion flows for the  Figure 13 show the variation of Re λ versus the viscosity ratio for the two volume fractions considered, α = 0.03 and α = 0.1. As expected, Re λ decreases with the viscosity ratio. Significant variations in Re λ are observed already for small volume fractions, the effects being amplified for α = 0.1. In the insets of the same figure, we can observe that λ (i.e. the local variations of k/ε) does not increase linearly with µ d /µ c , indicating increased velocity fluctuations for the dispersed phase at lower viscosity.
To better quantify the variations of the flow gradients, we show the phase-averaged (see eq. 2.8) enstrophy ω 2 m in Figure 14, normalized by the single-phase values from SP1. The viscosity ratio strongly affects enstrophy in the dispersed phase, while the magnitude in the carrier phase is almost constant. Further, smaller variations can be observed when changing the volume fraction from 0.03 to 0.1. For µ d /µ c 1, the enstrophy in the dispersed phase goes approximately as ω 2 c ∝ −log(µ d /µ c ). As the viscosity of the dispersed phase becomes larger, µ d > µ c , ω d decreases below the average value of the single-phase flow and tends towards zero, as high viscosity dampens velocity fluctuations in the dispersed phase. It is worth noting that, for incompressible flows, the energy dissipation rate can be defined as ε ≡ ν|ω i | 2 ; however, when phase averaging, the two formulations differ for by a term proportional to ∂ ii p.
We now discuss the influence of viscosity ratio on the compensated energy spectra, shown in Figure 15(a) for α = 0.03 and in panel (b) for α = 0.1. Similarly to previous observations for Figure 7, the Hinze scale shows, to a good approximation, the pivoting point, below which energy increases with respect to the single-phase spectra. Differences in the inertial subrange are hardly observable for α = 0.03, while they become more prominent when the volume fraction is increased, see panel (b). Analysis of the data for κ < κ H , reveals that the simulations with a dispersed phase present less energy than the single phase case. In the dissipative range the trend emerges more clearly. As κ > κ H , the less viscous the dispersed phase, the more energy is injected in the smaller scales. As discussed in the previous section, energy reduces at large scales and increases at small scales when increasing the volume fraction α. Figure 16 shows the DSD for all configurations with different viscosity ratios. As for the data in Section 3.1, we also display the -3/2 power-law, which well describes the distribution of small droplets, and the d −10/3 law from Garrett et al. (2000) for larger Single phase   droplets. In this range, d > d H , the -10/3 law is observed only in a limited region of the spectrum. As noted previously, this is most likely due to the low volume fraction considered. The variation of µ d has an influence on large droplets, as higher viscosity in the dispersed phase increases the probability of formation of these large droplets. This was also observed qualitatively in Figure 12 and confirms previous findings (Roccon et al. 2017). We present the SBS energy budget for α = 0.1 in Figure 17. The results for α = 0.03 show similar trends, see Appendix A for the details. Following the same scheme as in the previous section, we depict in panel (a) the energy balance for case V22, when µ b /µ c = 0.1. Similarly to previous observations, the dispersed phase absorbs energy at large scales and redistributes it to small scales, that is the presence of the interface provides an alternative path for energy transfer from small to large wavenumbers and no inverse cascade is observed. The non-linear energy transfer T , see panel (b), displays a weak sensitivity to the viscosity ratio (almost negligible for α = 0.03 as shown in Figure 26 in Appendix A). Thus, the differences in the Re λ and energy spectra discussed above are not associated with an extension of the inertial range. For wavenumbers larger than that of the forcing, the non-linear energy transfer is higher at large scales and lower at small scales than for the single-phase case. Panel (c) in Figure 17 show the energy transport due to the surface tension term, S σ . As µ d /µ c increases, the wavelength where the positive energy transport is maximum shifts to larger scales. This behavior is possibly due to increased coalescence for high µ d /µ c ( as discussed later in this section). We also observe that with decreasing viscosity ratio, µ d /µ c < 1, the curves tend to collapse, as the data for µ d /µ c = 0.1 and µ d /µ c = 0.01 are approximately overlapping. At the injection scale, κ = 2, almost all cases behave similarly. At intermediate wavelengths, the lower the viscosity ratio, the higher the energy absorbed by the surface tension forces. As previously observed, the Hinze scale represents, to a good approximation, the point where the energy transfer towards small scales by the surface tension term S σ is maximum. All these observations apply to the two values of α considered, see also Appendix A. A note should be made on the flow with the highest viscosity ratio: in this case, the energy is not transferred down to the dissipative range. A qualitative explanation is given by the following scenario. When the interface interacts with a sufficiently large vortex in the carrier phase, it tends to deform and, in doing so, absorbs energy through the work of the interfacial stresses. The deformation of the interface induces shear in the dispersed phase which is opposed by viscous forces. A higher viscosity in the dispersed phase will therefore dump larger and more energetic structures, reducing the energy available at small scales.
Finally, panel (d ) of the same figures shows the transfer function of the energy dissipation term D. We observe that simulations with higher viscosity of the dispersed phase dissipate more energy at large scales, hence dumping turbulence in the inertial range, as expected for more viscous flows. This trend is maintained until the dissipative range, where, instead, a lower viscosity ratio produces higher dissipation. This causes the apparently paradoxical situation that, despite there is limited energy transport by the non-linear terms, dissipation is still active at smalls scales because of the energy brought by the interfacial stresses; this may suggest the need of a specific definition of dissipative range for multiphase flows.
Next, we discuss the influence of the viscosity ratio on the phase-averaged energy budget, shown in Figure 18 for α = 0.1. As for the SBS balance, the same analysis for α = 0.03 can be found in Appendix A, as the variation of volume fraction does not significantly chance the underlying physical process. The production density (green symbols in panel a) shows only slight variations with viscosity ratios in the carrier phase, whereas it increases in the dispersed phase when its viscosity increases; in particular, P d < P c when µ d < µ c . A similar trend is observed for the dissipation rate (red symbols in panel a), when the differences between dispersed and carrier phases become more evident. In this case, the dissipation in the dispersed case increases with its viscosity until µ d /µ c = 100, when it decreases because of the lower energy transferred to smaller scales inside the droplets. The transport terms, T µ and T p in panel (b), indicate that Single phase Figure 19: PDF of velocity fluctuations u (panel a), vorticity ω (b) and energy dissipation (c). All quantities are normalized by their standard deviation. The data pertain cases V2x in Table 1, with α = 0.1.
energy is always transferred from the carrier to the dispersed phase. Both terms increase in magnitude when decreasing the viscosity ratio, indicating that energy needs to be supplied to the dispersed phase to sustain turbulence when viscous forces are increasing. The pressure transport is the preferential path for energy transfer from the carrier to the dispersed phase for low and moderate values of µ. For the case with largest viscosity of the dispersed phase, the transfer due to pressure forces becomes lower than that associated to viscous forces. Finally, we consider the effect of the viscosity ratio on the PDFs of velocity, vorticity and dissipation rate, see Figure 19 for α = 0.1, while data for α = 0.03 can be found in Appendix A. A low viscosity in the dispersed phase generates larger velocity fluctuations (see also Figure 14), hence the tails of PDFs are more evident for small values of µ d /µ c in panel (a). Interestingly, when the viscosity ratio increases above unity, velocity fluctuations in the dispersed phase are quenched, the standard deviation decreases and the statistics are closer to those of the single-phase reference case. The distributions of the normalized vorticity are shown in panel (b). As for the velocity fluctuations, a higher viscosity in the dispersed phase decreases the intermittency and the distributions approach the single-phase values. For µ d /µ c < 1, intermittency increases and the tail of the distribution are more evident; nonetheless they can still be fitted with decaying exponentials. As previously observed for varying α, the pseudo-Gaussian part of the vorticity PDF collapse for all cases.
The PDF of the energy dissipation show almost no alteration between cases with different viscosity ratios, while intermittency is strongly increased with respect to the single-phase case. Due to the normalization with σ ε , the curve collapse indicates that variations induced by µ d /µ c of the small-scale dynamics are negligible.
To conclude, the turbulence is significantly affected by variations of the viscosity ratio already at small volume fractions. Higher viscosity of the dispersed phase dampens the small-scale structures because of higher viscous dissipation at all wavelengths. For emulsions with viscosity of the dispersed phase lower than that of the carrier phase, the activity at small scales increases and so does intermittency. The surface tension term S σ significantly contributes to the transfer of energy to the smallest scales in this case.

Influence of Weber number
The influence of the surface tension coefficient, expressed through the large scale W e L number, is examined in this section. As discussed in literature (Komrakova et al. 2015;Roccon et al. 2017;Mukherjee et al. 2019), the combination of volume fraction, surface tension coefficient and energy injection scale, L, has to be accurately chosen because the HIT configuration is very sensitive to coalescence. Furthermore, high W e L may generate an excess of unresolved droplets, significantly affecting the results. Therefore, all the simulations discussed in this section are performed at α = 0.03, while the forcing is maintained at κ 0 = 2. The cases discussed in this section are BE1 and the series W1x with reference to Table 1, covering a significant range of W e L , from 10.6 to 106.5. In Figure 20 we show a render of the flow for different values of W e L . As expected, at low W e L we observe the appearance of large liquid structures due to higher surface tension forces. At high W e L , on the other hand, the dispersed phase undergoes severe fragmentation. The presence of small droplets resulting from fragmentation can be observed in all cases. We start by studying the global behavior of the flow through the integral quantities in Figure 21. Re λ , reported in panel (a), shows an almost linear increment both with the Taylor scale λ and W e L (represented with colors). Unlike previous observations in Section 3.1 where the increase of Re λ was mostly due to local variations of the ratio k/ε, decreasing surface tension also lowers the volume-averaged energy dissipation, as shown in the inset. These findings are in agreement with the results on turbulent emulsions in Rosti et al. (2020). As the viscosity ratio µ d = µ c is constant, the decrease of the dissipation is caused by lower enstrophy levels, as it can be appreciated from the data for the carrier phase in panel (b). The behavior of the enstrophy of the dispersed phase is less intuitive, exhibiting a non-monotonic behavior, and will be addressed later when discussing the phase-averaged energy balance. Figure 22(a) shows the compensated energy spectra at different W e L . As we are varying the surface tension, the Hinze scale varies in each case (see vertical dotted lines of corresponding color). As mentioned before, the Hinze scale defines with good approximation the spectra pivoting point. As previously discussed, energy is reduced at larger scales in the inertial range and increases at smallest scales. With increasing W e L , higher energy is observed at high wavelengths. Figure 22(b) shows the droplet-size-distribution for all the W e L under investigation. As for panel (a) we show the Hinze scale for each case with vertical dotted lines. Again we observe that the -10/3 power-law from Garrett et al. (2000) provides a reasonable description for the largest droplets, d > d H ; as we increase σ, i.e. low W e L , larger droplets may appear, as expected by the increased cohesion forces. In this case, the energy required to breakup large droplets is only available in large eddies. As their turnover time is in the order of T , large droplet breakup becomes a rare event and the distributions are more noisy, so that it is more difficult to identify a clear trend. In addition, by reducing W e L , the distribution becomes more irregular for d < d H as most of the dispersed phase is in large droplets. The effects of W e L can be better described by the scale-by-scale analysis, shown in Figure 23. The complete energy balance is shown for case W11 (W e L = 10.6) in panel (a). Unlike the the cases shown previously for W e L = 42.6, the surface tension energy transfer S σ is more uniform through the different scales and its effects are globally less evident. To deepen the analysis, we display the non-linear energy transfer function T for each case at different W e L in figure Figure 23(b). At the injection wavelength κ = 2, no major differences are observed when varying the surface tension. At small wavelengths, κ > 2, we observe that the energy transfer by the non-linear term increases with W e L , compensating for the effect of the energy absorption from the surface tension. The energy transfer at smaller scales, after the peak, increases with σ, approaching the values of the single-phase flow.
The energy transfer via the interfacial stresses, S σ , is shown in panel (c) of the same figure. The energy is again absorbed at large scales and distributed at small scales. Flows with small W e L absorb more energy at small wavenumbers and the transmission of energy (i.e. positive S σ ) is smeared over a higher range of scales, hence the peak (max(S σ )) is also less evident. For all W e L investigated, the surface tension term S σ transfers energy also within the dissipative range at small scales, where the transport from non-linear terms has become negligible.
The energy dissipation D, panel (d ) of Figure 23, decreases with W e L at large and intermediate scales, as energy is partially absorbed by S σ . The amplitude of the dissipation rates becomes however almost independent of the Weber number at the smallest scales. Further, as previously observed, the presence of the dispersed phase delays the onset of the dissipative range.
The phase-averaged energy balance from simulations with different Weber number is shown in Figure 24. Both production and dissipation (panel a) are found to decrease in the carrier phase when increasing W e L , while the former increases and then decreases in the carrier phase. This can be possibly related to the droplet size distributions: decreasing the droplet size increases the internal dissipation, which may explain the behavior at the lower Weber examined. On the other hand, high deformability decreases the dissipation close to the interface, which may explain the decrease at the largest W e L . For all values of W e L considered, the dispersed phase extracts kinetic energy from the carrier phase, as T c > 0 (panel b). The decrease of surface tension forces results in a monotonic decrease of the viscous transfer and an increase of the pressure transport for the dispersed phase. Consistently, dissipation is always higher in the dispersed phase, while it decreases in the carrier phase when increasing W e L .
The analysis of the PDF for velocity, vorticity and dissipation are finally shown in panels (a), (b) and (c) of Figure 25. Strong variations are induced in all PDFs, showing that indeed a more rigid interface favors the appearance of extreme events. Since a more deformable interface offers lower resistance to the propagation of velocity disturbances from one phase to the other, less modifications of the PDFs with respect to the singlephase case can be expected at higher W e L (see also Rosti et al. 2020). This is indeed observed in all PDFs, where the distributions are seen to approach the single-phase one when increasing W e L . Nevertheless, rare events are still evident also at the largest Weber considered, especially for the energy dissipation. Vorticity shows, again, that the psuedo-Gaussian part of the distribution is identical for all W e L , while the exponentially decaying tails display strong variations. All quantities are normalized by their standard deviation. The data pertain cases W1x, BE1 and SP1 in Table 1.

Conclusions
In this work we discuss how volume fraction, viscosity ratio and Weber number influence HIT in emulsions. The analyses are performed at different levels of details, spanning from phase averaged balances to SBS energy transfer in spectral space. Some observations are common to all configurations and highlight some fundamental physical effects introduced by the dispersed phase. Here, we first consider these different aspects and then discuss the modulation introduced by the variation of material properties.

Spectra and SBS energy balance.
In all simulations with a dispersed phase, the energy decreases at large scales and increases at small scales, corroborating previous findings (Ten Cate et al. 2004;Perlekar et al. 2014;Dodd & Ferrante 2016;Mukherjee et al. 2019;Rosti et al. 2020;Olivieri et al. 2020a). Interestingly, this behaviour applies to both solid and liquid dispersed phases in HIT. Furthermore, the pivoting point of the energy spectra is found to be described, with a good approximation, by the Hinze scale. This has also been observed in binary mixtures (Perlekar et al. 2014) and emulsions (Mukherjee et al. 2019) and is here extended to several operating conditions.
In general, the mechanisms of energy transport are modified as follows: the transfer by the non-linear advection terms decreases, as the surface tension forces absorb energy at large scales. In an emulsion, energy is transferred to small scales also by the surface tension force, well within the dissipative range of the corresponding single-phase flow, forcing the viscous dissipation to be active at even smaller scales. No inverse cascade has been observed in the present simulations.
The general idea, according to which coalescence and breakup are responsible for modifications of the energy spectra seems to only partially explain our observations. In fact, according to this hypothesis, significant deviations should be observed when comparing spectra for different volume fractions. Here, instead, we observe the largest deviations in the energy spectra, in particular at small scales. when varying the viscosity ratio. This issue may be further addressed in future studies if coalescence is inhibited, reduced or controlled numerically.

Effects of the dispersed phase on the dissipative range
The classical "far" dissipative range (κ ∼ κ max ), where both non-linear energy transport and energy dissipation of the SBS budget are zero, is lost when a dispersed phase is introduced. In multiphase flows, despite the non-linear energy transfer vanishes at certain small scales, the energy dissipation does not because energy is brought to these smaller scales by the action of the surface tension. As discussed above, energy dissipation is thus forced to extend towards smaller scales, overall increasing the range of wavelengths where there is activity. In other words, this increased activity at small scale translates also into an extension of the dissipative range, with the non-linear transport substituted by the surface-tension transport.
It is important to understand how the scaling in the inertial range might be affected by these modifications of the dissipation range. From a practical viewpoint, the results in Section 2.4 shows that increasing the mesh resolution does not result in significant alterations of the inertial range, indicating that a relevant analysis of the inertial range dynamics is still possible even in simulations where the surface tension terms are slightly under-resolved at small scale. Nevertheless, resolving the dissipative range is important for a complete discussion of the SBS budget and e.g. the DSD; understanding turbulence at small scales in multiphase flows remains therefore a relevant question also from a fundamental point of view.

Flow intermittency
We have observed that the presence of a dispersed phase increases intermittency, unless the dispersed phase is highly viscous. The probability of detecting rare events increases, mainly for energy dissipation and vorticity, as shown here by the PDF analysis.
In particular, at higher volume fractions and constant µ d /µ c and W e L , the exponent describing the distribution tail exponential decay is independent of the volume fraction α. The onset of the exponential tail (hence the probability of an extreme event) is, on the other hand, affected by α, proving that these events are mostly occurring at the interface. This is a confirmation of the observations in Dodd & Ferrante (2016) on the increased energy dissipation at the interface. The variation of the exponential tail for both energy dissipation rate and vorticity at different µ d /µ c and W e L reveals that intermittency is significantly affected by the fluid properties. In cases with high µ d and low surface tension, the vorticity intermittency is attenuated and similar to the single-phase cases. On the other hand, the dissipation seems to be always affected by the multiphase nature of the flow.

Droplet statistics
In all the conditions analyzed, the droplet-size distributions show both the -3/2 exponential scaling from Deane & Stokes (2002) for the small droplets and the -10/3 from Garrett et al. (2000) for the larger ones, confirming and extending the previous findings of Mukherjee et al. (2019) to a significant number of different configurations. Moreover, employing a VOF approach, and its known mass conserving properties, allows to extend the -3/2 scaling to significantly small droplets.
The power-law d −10/3 well describes the distributions of larger droplets when the volume fraction is below 10%, with only a small loss in accuracy for higher values of α, in agreement with the assumption of negligible coalescence in Garrett et al. (2000). Although this power law was obtained under the assumption of a dilute dispersed phase, recent works based on a diffuse-interface approach report the same scaling in the presence of coalescence ( estimate through accurate sharp-interface simulations a similar exponent, −3, so that it might be difficult to have a clear distinction on the different effects. Finally, we show that the estimate of the Hinze scale as transition point between the two power laws is less accurate for α > 0.1, suggesting different model coefficients may be needed when coalescence is relevant.

Role of the fluid properties
Our analyses demonstrate that the volume fraction α is the parameter that mostly modifies the energy fluxes in the flow; yet, increasing the volume of the dispersed phase does not change the underlying physics. This is notably documented in Section 3.1 where we show that the amount of total interface area determines the energy transport across scales. Moreover, the simulation data reveal that the energy transfer via surface tension forces is enhanced at low viscosity ratios, while high viscosity in the carrier phase inhibits the propagation of vortices through the interface, hence reducing the overall energy transport. Changing the Weber number amounts to modulating the pivoting frequency below which energy transfer through surface tension is directed towards smaller scales. In particular, as the dispersed phase is less deformable, the energy absorption from the dispersed phase occurs at larger scales, and turbulence is progressively reduced. In fact, as the surface tension increases, more energy is required to deform the droplets, an energy which can only be found in large-scale eddies.
To study the role of the viscosity ratio, see Section 3.2, we consider values ranging from 10 −2 (a value typical of bubbles) to 10 2 (typical of droplets).
The analysis reveals that for µ d /µ c 1, Re λ increase significantly, due to the lower viscosity in the dispersed phase. The scale-by-scale energy budget shows that the interfacial and non-linear transport terms are not strongly affected at these low viscosity ratios. For µ d /µ c > 1, on the other hand, the turbulence in the dispersed phase is reduced, which implies a significantly smaller Re λ , below the value of the single-phase case. In these cases, the energy transfer induced by the interfacial stresses is significantly reduced, suggesting that large differences may be found in liquid-gas and gas-liquid emulsions. The droplet-size distribution does not show strong differences, although larger droplets are more likely to be generated by a more viscous dispersed phase. Note, as discussed above, that the viscosity ratio has a significant impact on the flow intermittency.
Finally, we have examined the role of the large-scale Weber number W e L . At low W e L , coalescence is more likely to occur, hence there is a higher probability to find large droplets. Nevertheless, the Hinze scale proves to be an accurate estimation of the transition between the -3/2 and -10/3 for all the cases analyzed. Changing W e L and thus the droplet size distribution also affects the energy transport across scales by the surface tension forces. Specifically, when decreasing W e L the energy injection from interfacial tension moves to larger scales. the single-phase case are small, see panel (a) of the figure. The energy transport due to surface tension (panel c) is attenuated at high viscosity ratios and shifts towards small wavelengths due to increased coalescence (see Section 3.2). Finally, energy dissipation, panel (d ), shows again limited variations due to reduced volume fraction, although, it can be observed again that the small scale energy transfer is unaffected at high viscosity ratios. The phase-averaged energy balance in Figure 27 shows only weak variations with respect to the cases at α = 0.1 in Figure 18. Again, we notice that energy dissipation in the dispersed phase increases at higher µ d , while energy is always transferred from the carrier to the dispersed phase, as for α = 0.1.
We finally present the PDFs of velocity, vorticity and energy dissipation in Figure 28(a,b,c). Again, small variations can be observed with respect to cases at α = 0.1 (Figure 19). For vorticity and energy dissipation, we report lower probability to observe rare events at lower volume fraction, as discussed in Section 3.1.