1. Introduction
Droplet-laden turbulent flows are ubiquitous in nature and industry (Jähne & Haußecker Reference Jähne and Haußecker1998; Karsa Reference Karsa1999; Schramm, Stasiuk & Marangoni Reference Schramm, Stasiuk and Marangoni2003; Kralova & Sjöblom Reference Kralova and Sjöblom2009; Dickinson Reference Dickinson2010). A few examples include the capture of atmospheric CO
$_2$ at the surface of seas and oceans, which is mediated by the entrainment of air bubbles by breaking waves (Merlivat & Memery Reference Merlivat and Memery1983; Deane & Stokes Reference Deane and Stokes2002; Pereira et al. Reference Pereira, Ashton, Sabbaghzadeh, Shutler and Upstill-Goddard2018), or the dynamics of liquid jets and sprays that is of fundamental interest for combustion, cooling, irrigation and firefighting (Mugele & Evans Reference Mugele and Evans1951; Faeth, Hsiang & Wu Reference Faeth, Hsiang and Wu1995; Herrmann Reference Herrmann2011; Canu et al. Reference Canu, Puggelli, Essadki, Duret, Menard, Massot, Reveillon and Demoulin2018; Kooij et al. Reference Kooij, Sijs, Denn, Villermaux and Bonn2018). These flows rarely involve pure fluids: instead, they often include small amounts of impurities that may act as surface-active agents (surfactants). Surfactants are compounds that assemble at the fluid interface and modify the local surface tension. Even small amounts of surfactant can drastically change the flow behaviour, making their presence crucial in many practical scenarios (Koshy, Das & Kumar Reference Koshy, Das and Kumar1988; Dobbs Reference Dobbs1989; Sjoblom Reference Sjoblom2005; Takagi & Matsumoto Reference Takagi and Matsumoto2011).
The interface between the carrier phase and the dispersed phase, i.e. the droplets, serves as a conduit for various physical and chemical exchanges, such as heat, vapour (Scapin, Costa & Brandt Reference Scapin, Costa and Brandt2020; Onofre Ramos et al. Reference Onofre Ramos, Zhao, Bouali and Mura2022), solutes and aerosols (de Leeuw et al. Reference de Leeuw, Andreas, Anguelova, Fairall, Lewis, O'Dowd, Schulz and Schwartz2011). The rate of these exchanges is determined by the product of the interfacial flux and the interfacial area, emphasizing the pivotal role of interfacial characteristics. There has been considerable effort to estimate these exchanges via empirical correlations (Akita & Yoshida Reference Akita and Yoshida1974; Kelly & Kazimi Reference Kelly and Kazimi1982; Delhaye & Bricard Reference Delhaye and Bricard1994) or via population balance equations relying on the droplet size distribution and on droplet breakage and coalescence models (Luo & Svendsen Reference Luo and Svendsen1996; Babinsky & Sojka Reference Babinsky and Sojka2002; Andersson & Andersson Reference Andersson and Andersson2006a,Reference Andersson and Anderssonb; Martínez-Bazán et al. Reference Martínez-Bazán, Rodríguez-Rodrìguez, Deane, Montañes and Lasheras2010; Chan et al. Reference Chan, Dodd, Johnson, Urzay and Moin2018; Chan, Johnson & Moin Reference Chan, Johnson and Moin2021; Gaylo, Hendrickson & Yue Reference Gaylo, Hendrickson and Yue2023). Hence, in this study we aim to answer the question, ‘what is the interfacial area of surfactant-laden droplets in turbulence?’ Indeed, we measure the interfacial area of each droplet and discover the presence of two universal regimes: the interfacial area of small droplets is proportional to the square of their characteristic size, whereas, for large droplets, it is proportional to the cube of their characteristic size. The information on the individual interfacial area, combined with the droplet size distribution, provides an estimate of the total interfacial area available. The two different regimes are directly linked to the shape of the droplets: small droplets are spheroid-like or ellipsoid-like, whereas large droplets take long, filamentous shapes. We find that the length scale separating these two regimes is the Kolmogorov–Hinze scale, defined as the maximum size of a droplet that is not broken apart by turbulent fluctuations; droplet breakage becomes prevalent for droplets larger than the Kolmogorov–Hinze scale. The concept of the Kolmogorov–Hinze scale originates from the works of Kolmogorov (Reference Kolmogorov1949) and Hinze (Reference Hinze1955), who applied Kolmogorov's (Reference Kolmogorov1941) assumptions to droplets in turbulence. Some recent studies, however, have disputed the theoretical framework upon which Hinze's theory is based and, hence, the relevance of the Kolmogorov–Hinze scale: Qi et al. (Reference Qi, Tan, Corbitt, Urbanik, Salibindla and Ni2022) showed that droplets interact with eddies of a range of length scales, rather than solely with eddies of a size similar to the droplet, and Vela-Martín & Avila (Reference Vela-Martín and Avila2022) showed that droplet breakup does, in fact, occur below the Kolmogorov–Hinze scale. Motivated by these studies, we ask, ‘does the Kolmogorov–Hinze framework hold in surfactant-laden flows?’ As we shall see, the answer is yes. That is, our numerical simulations show that the Kolmogorov–Hinze scale can be used as a key parameter to describe the morphology of surfactant-laden droplets in turbulence, and that the values obtained using Hinze's original formulation (Hinze Reference Hinze1955) are in good agreement with a more recent formulation, which uses the work done by the interface to define a length scale for the droplet phase (Crialesi-Esposito, Chibbaro & Brandt Reference Crialesi-Esposito, Chibbaro and Brandt2023). Furthermore, we show that Hinze's theory can be extended to surfactant-laden droplets, provided that a suitable value of the surface tension is selected. In our set-up, surfactant effects on the morphology of the droplets can indeed be well approximated using an averaged value of the surface tension, thereby maintaining the simplicity and efficacy of the Kolmogorov–Hinze framework.
Beyond an average reduction in surface tension, surfactants introduce more intricate dynamics into the flow: a surfactant is an additional phase that is transported by the local flow and by motion and deformation of the interface, and that reduces the value of surface tension according to its local concentration. This can lead to an inhomogeneous value of surface tension over the interface of the droplets, giving rise to Marangoni stresses, i.e. stresses that act tangentially to the interface and originate from surface tension gradients. Marangoni stresses have been shown to be crucial in hindering and preventing coalescence (Dai & Leal Reference Dai and Leal2008; Soligo, Roccon & Soldati Reference Soligo, Roccon and Soldati2019b), reducing the rising velocity of bubbles (Takagi & Matsumoto Reference Takagi and Matsumoto2011; Elghobashi Reference Elghobashi2019), and in the build-up of bubble layers in wall-bounded flows (Lu & Tryggvason Reference Lu and Tryggvason2008; Takagi & Matsumoto Reference Takagi and Matsumoto2011; Tryggvason & Lu Reference Tryggvason and Lu2015; Lu, Muradoglu & Tryggvason Reference Lu, Muradoglu and Tryggvason2017; Ahmed et al. Reference Ahmed, Izbassarov, Costa, Muradoglu and Tammisola2020). An increase in the drag coefficient has been reported when adding surfactant to wall-bounded bubbly flows (Takagi, Ogasawara & Matsumoto Reference Takagi, Ogasawara and Matsumoto2008; Takagi et al. Reference Takagi, Ogasawara, Fukuta and Matsumoto2009; Verschoof et al. Reference Verschoof, Van Der Veen, Sun and Lohse2016): surfactant reduces the size of the droplets and causes a lower drag reduction compared with surfactant-free cases. Our study examines a statistically stationary, homogeneous and isotropic multiphase flow at a moderate Reynolds number (see figure 1). This type of flow does not allow for the build-up of large-scale surfactant gradients commonly found in up-flow and down-flow set-ups, where the velocity difference between the carrier and the dispersed phase generates and maintains a surfactant gradient along the interface of the droplets or bubbles (Takagi et al. Reference Takagi, Ogasawara and Matsumoto2008; Lu et al. Reference Lu, Muradoglu and Tryggvason2017). For this reason, we expect the effect of Marangoni stresses to be localized. In this study, we ask, ‘what is the effect of Marangoni stresses on droplet morphology in homogeneous homogeneous isotropic isotropic turbulence?’ As we shall show, surfactant effects on droplet morphology in our set-up can be summarized as an average surface tension reduction, and the effect of Marangoni stresses can only be appreciated by analysing the local flow dynamics at the interface.
Snapshots of the simulated domain in the cases (a) S05 and (b) S20 showing the interface of the droplets; the scale bar shows the Kolmogorov–Hinze scale for each case.
To investigate the complex dynamics of clean and surfactant-laden droplets in turbulence, we use direct numerical simulations. In recent years there has been a consistent growth in the number of numerical studies on multiphase turbulence, supported as well by the increased availability and capability of high-performance computing infrastructures. Multiphase turbulence is characterized by a wide range of scales, from the smallest, molecular-size interfacial scale, to the smallest flow scales – the Kolmogorov length scale – and up to the large-scale structures of the flow. The separation of scales usually spans over about eight to ten orders of magnitude, while direct numerical simulations on leading-edge high-performance computing systems can simulate about four orders of magnitude (Soligo, Roccon & Soldati Reference Soligo, Roccon and Soldati2021). The common choice is to simulate the larger scales of the multiphase flow, from the large-scale structures down to about the Kolmogorov length scale, and introduce models for the smaller-scale physics (Tryggvason et al. Reference Tryggvason, Thomas, Lu and Aboulhasanzadeh2010, Reference Tryggvason, Dabiri, Aboulhasanzadeh and Lu2013; Soligo et al. Reference Soligo, Roccon and Soldati2021). Several authors have devoted their attention to the development of models for the unresolved scales to be used in direct numerical simulations and large eddy simulations. In their original works Kolmogorov (Reference Kolmogorov1949) and Hinze (Reference Hinze1955) identified the maximum size of a non-breaking droplet in turbulence. Experimental measurements (Garrett, Li & Farmer Reference Garrett, Li and Farmer2000; Deane & Stokes Reference Deane and Stokes2002) later showed the existence of two different regimes in the droplet size distribution, separated by the Kolmogorov–Hinze scale. These works laid the foundations for the understanding of the dynamics of the droplets and for the development of sub-grid-scale models for the interfacial dynamics (Herrmann Reference Herrmann2013; Xiao, Dianat & McGuirk Reference Xiao, Dianat and McGuirk2014; Evrard, Denner & van Wachem Reference Evrard, Denner and van Wachem2019). The power-law scaling exponents for the two regimes measured in experiments have been confirmed as well by numerical simulations: 
$-10/3$ for the breakage-dominated regime (Perlekar et al. Reference Perlekar, Biferale, Sbragaglia, Srivastava and Toschi2012; Skartlien, Sollum & Schumann Reference Skartlien, Sollum and Schumann2013; Deike, Melville & Popinet Reference Deike, Melville and Popinet2016; Mukherjee et al. Reference Mukherjee, Safdari, Shardt, Kenjereš and Van den Akker2019; Rosti et al. Reference Rosti, Ge, Jain, Dodd and Brandt2019b; Soligo, Roccon & Soldati Reference Soligo, Roccon and Soldati2019a; Crialesi-Esposito et al. Reference Crialesi-Esposito, Chibbaro and Brandt2023) and 
$-3/2$ for the coalescence-dominated regime (Rivière et al. Reference Rivière, Mostert, Perrard and Deike2021; Crialesi-Esposito et al. Reference Crialesi-Esposito, Chibbaro and Brandt2023). The droplet size distribution and the population balance equation for the droplets are fundamental tools in the modelling of droplets of similar size to the grid resolution and smaller. Perlekar et al. (Reference Perlekar, Biferale, Sbragaglia, Srivastava and Toschi2012) correlated the instantaneous Weber number of droplets to their deformation, showing that large Weber numbers correspond to strongly deformed droplets. They simulated emulsions at increasing volume fractions and proved the validity of the Kolmogorov–Hinze scale at low volume fractions; they reported deviations from the Kolmogorov–Hinze theory for more concentrated emulsions, possibly due to the effect of coalescence, which was neglected in the original works by Kolmogorov and Hinze. Vela-Martín & Avila (Reference Vela-Martín and Avila2022) showed that drop breakage is a memoryless process, i.e. the relaxation time of the droplet is much lower than its expected lifetime. In the same work, they debated the validity of the Kolmogorov–Hinze scale as an absolute threshold between breaking and non-breaking droplets, arguing that all droplets will eventually break apart, provided there is enough time for breakage. It was shown also that, in the absence of coalescence, the breakage rate depends on the Weber number alone. Gaylo et al. (Reference Gaylo, Hendrickson and Yue2023) investigated the fragmentation of bubbles in statistically stationary homogeneous isotropic turbulence and characterized several fundamental time scales of the bubbles: the relaxation time, the expected lifetime and the time needed for the largest bubble to break down to the Kolmogorov–Hinze scale. It is now well known that the presence of a dispersed phase strongly modifies the dynamics of the flow at all scales: the interface extracts energy at the large flow scales and re-injects it into the flow at much smaller scales, competing with the classic turbulent energy cascade (Mukherjee et al. Reference Mukherjee, Safdari, Shardt, Kenjereš and Van den Akker2019; Perlekar Reference Perlekar2019). This effect is reflected in a deviation from the 
$k^{-5/3}$ power-law scaling of the turbulent kinetic energy spectrum. When there is a considerable slip velocity between the droplet or bubble and the carrier fluid, a completely different scaling of the energy spectrum, 
$k^{-3}$, has been reported in numerical simulations (Roghair et al. Reference Roghair, Mercado, Annaland, Kuipers, Sun and Lohse2011; Pandey, Ramadugu & Perlekar Reference Pandey, Ramadugu and Perlekar2020; Paul et al. Reference Paul, Fraga, Dodd and Lai2022) and experiments (Mercado et al. Reference Mercado, Gomez, Van Gils, Sun and Lohse2010; Prakash et al. Reference Prakash, Mercado, van Wijngaarden, Mancilla, Tagawa, Lohse and Sun2016).
As part of our analysis of droplet morphology, we measure the Euler characteristic of the droplet interfaces. Similarly to the number of droplets in a flow, the Euler characteristic of an interface is an integer-valued topological invariant, and any change in its value requires a splitting or merging of interfaces. Despite its physical significance, the Euler characteristic has only recently been applied to multiphase flows: Dumouchel, Thiesset & Ménard (Reference Dumouchel, Thiesset and Ménard2022) linked the Euler characteristic to the Gaussian curvature of the droplets and used it to parametrize the morphology of liquid droplets undergoing breakup. The Euler characteristic is commonly used in characterizing the sintering of metal powders (DeHoff, Aigeltinger & Craig Reference DeHoff, Aigeltinger and Craig1972; Mendoza et al. Reference Mendoza, Thornton, Savin and Voorhees2006), classifying lung tissues (Boehm et al. Reference Boehm, Fink, Attenberger, Becker, Behr and Reiser2008) and correcting MRI scans of the human brain (Yotter et al. Reference Yotter, Dahnke, Thompson and Gaser2011). In this paper we use the Euler characteristic to count the number of voids and handles on the droplets, demonstrating that the large droplets are made up of extremely interconnected filaments.
The paper is structured as follows. We introduce the numerical method and the computational set-up adopted for the simulations in § 2. Our findings are reported in § 3, where we first focus on the morphology of the dispersed phase, and then on the statistics of the local flow around the droplets. Finally, § 4 summarizes the main results presented in the present work.
2. Numerical model
We solve a system of equations including the momentum (2.1) and mass (2.2) conservation, the volume of fluid (2.3) and surfactant (2.4) transport equations to simulate the dynamics of an ensemble of breaking, coalescing and deforming finite-size droplets in homogenous isotropic turbulence. The two phases, the carrier fluid and the dispersed phase (i.e. the droplets), have the same density 
$\rho$ and dynamic viscosity 
$\mu$. Thus, the governing equations are
\begin{gather} \rho \frac{\partial \boldsymbol{u}}{\partial t} + \rho \boldsymbol{u}\boldsymbol{\cdot} \boldsymbol{\nabla} \boldsymbol{u} ={-}\boldsymbol{\nabla} p +\boldsymbol{\nabla} \boldsymbol{\cdot} [ \mu ( \boldsymbol{\nabla}\boldsymbol{u}+\boldsymbol{\nabla} \boldsymbol{u}^T)] + \boldsymbol{\nabla} \boldsymbol{\cdot} (\boldsymbol{\tau}_{\!c} f_\sigma) + \boldsymbol{f_s}, \end{gather}
\begin{gather}\boldsymbol{\nabla}\boldsymbol{\cdot} \boldsymbol{u}=0, \end{gather}
\begin{gather}\frac{\partial \phi}{\partial t} +\boldsymbol{\nabla} \boldsymbol{\cdot}(\boldsymbol{u}H)=0 , \end{gather}
\begin{gather}\frac{\partial \psi}{\partial t} +\boldsymbol{u}\boldsymbol{\cdot} \boldsymbol{\nabla} \psi=\boldsymbol{\nabla} \boldsymbol{\cdot} (M_\psi \boldsymbol{\nabla} \mu_\psi). \end{gather}
Here 
$\boldsymbol {f_s}$ is the spectral forcing used to sustain turbulence. We use the one-fluid approach, whereby the fluid velocity 
$\boldsymbol {u}$ and pressure 
$p$ are defined in both phases and continuous across the interface; the volume-of-fluid variable 
$\phi$ is used to define the instantaneous position of the interface. The volume of fluid can be understood as a colour function characterizing the local concentration of the dispersed phase: it is equal to 
$\phi =0$ in the carrier phase and to 
$\phi =1$ in the droplet phase. The volume-of-fluid method is an interface capturing method (Prosperetti & Tryggvason Reference Prosperetti and Tryggvason2009) where the concentration of each phase is transported using (2.3) and the interface is implicitly defined as the 
$\phi =0.5$ level. The effect of the interface on the flow is accounted for in the momentum equation via the surface tension forces: the Korteweg tensor 
$\boldsymbol{\tau}_{\!c}=\boldsymbol{\mathsf{I}}-\boldsymbol {n} \otimes \boldsymbol {n}$ (Korteweg Reference Korteweg1901) accounts for the position and shape of the interface, and the surface tension equation of state 
$f_\sigma$ defines the local value of surface tension. In the definition of the Korteweg tensor, 
$\boldsymbol{\mathsf{I}}$ is the identity matrix and 
$\boldsymbol {n}=-\boldsymbol {\nabla }\phi / \|\boldsymbol {\nabla } \phi \|$ is the unit-length, outward-pointing normal to the interface. The local surfactant concentration 
$\psi \in [0,1]$ is expressed as a fraction of the maximum surfactant concentration, which is usually determined by steric hindrance between surfactant molecules. Hence, 
$\psi$ is dimensionless in this formulation. To account for the effect of surfactant on surface tension, we use a modified Langmuir equation of state (Muradoglu & Tryggvason Reference Muradoglu and Tryggvason2008, Reference Muradoglu and Tryggvason2014; Soligo et al. Reference Soligo, Roccon and Soldati2019b) 
$f_\sigma = \max [\sigma _{min},\sigma _0(1+\beta _s \log (1-\psi )) ]$, where 
$\sigma _0$ is the reference surface tension of a clean (i.e. surfactant-free interface), 
$\sigma _{min}$ the minimum surface tension and 
$\beta _s$ the elasticity number. In the original formulation by Bazhlekov, Anderson & Meijer (Reference Bazhlekov, Anderson and Meijer2006), Pawar & Stebe (Reference Pawar and Stebe1996), the Langmuir equation of state provides a good fit at low to moderate surfactant concentration values; however, it fails to account for surfactant saturation dynamics at higher concentrations. Experimental measurements (reviewed by Chang & Franses Reference Chang and Franses1995) showed that, beyond a critical concentration of surfactant, surface tension no longer changes for increasing surfactant concentrations. Hence, to qualitatively account for surfactant saturation dynamics at the interface, we limit the surface tension to be greater than 
$\sigma _{min}$ at all points on the interface. In (2.1), surface tension forces act perpendicular and tangential to the interface: a capillary component (normal to the interface) and a tangential component – Marangoni stresses – proportional to the surface tension gradient. The tangential component is characteristic of surfactant-laden flows, where surface tension changes along the interface according to the local surfactant concentration.
The volume of fluid 
$\phi$ is transported using a simple advection equation (2.3). The multi-dimensional tangent of hyperbola for interface capturing (MTHINC) method (Ii et al. Reference Ii, Sugiyama, Takeuchi, Takagi, Matsumoto and Xiao2012) is used to reconstruct the local volume-of-fluid value 
$\phi$ starting from the cell-local indicator function 
$H$. The transport equation for the cell-local indicator 
$H$ is integrated over a control volume, yielding (2.3) (Ii et al. Reference Ii, Sugiyama, Takeuchi, Takagi, Matsumoto and Xiao2012; Rosti, De Vita & Brandt Reference Rosti, De Vita and Brandt2019a; Rosti et al. Reference Rosti, Ge, Jain, Dodd and Brandt2019b). To compute the surfactant chemical potential, we first calculate a signed-distance function 
$s$ and a smoothed colour function 
$\hat \phi$. A re-distancing equation is solved to compute the signed-distance function over the pseudo-time 
$\tau$,
\begin{equation} \frac{\partial s}{\partial \tau}=\text{sgn}(s_0) (1-\| \boldsymbol{\nabla} s\|),\end{equation}
where 
$\text {sgn}$ is the sign function and the initial guess 
$s_0$ is found as 
$s_0=(2\phi -1)0.75\varDelta$ (Albadawi et al. Reference Albadawi, Donoghue, Robinson, Murray and Delauré2013); this choice guarantees that the zero level of the signed-distance function always corresponds to the interface (Russo & Smereka Reference Russo and Smereka2000; De Vita et al. Reference De Vita, Rosti, Caserta and Brandt2019). The signed-distance function is updated at every time iteration as the volume of fluid is advected. Next, we compute the smoothed colour function as 
$\hat \phi =\tanh ({s}/{3\varDelta })$, which is bounded in 
$-1\le \hat \phi \le 1$ and where the smoothing width is set to three times the grid spacing 
$\varDelta$.
An advection–diffusion equation (2.4) is solved to track the surfactant concentration 
$\psi$ in the entire domain. We use a soluble surfactant: surfactant preferentially collects at the interface between the two fluids, but at the same time, it also dissolves in limited amounts in the bulk of the phases. We use a one-fluid formulation for the surfactant, in which the variable 
$\psi$ defines the surfactant concentration in the entire domain (i.e. in the bulk and at the interface). The adsorption (accumulation of surfactant from the bulk to the interface), desorption (release of surfactant from the interface to the bulk) and diffusion dynamics of the surfactant phase are determined by the chemical potential of the surfactant. The chemical potential is made, in order, of three contributions: a free energy of mixing term, an adsorption term and a bulk-penalty term (Engblom et al. Reference Engblom, Do-Quang, Amberg and Tornberg2013; Yun, Li & Kim Reference Yun, Li and Kim2014; Soligo et al. Reference Soligo, Roccon and Soldati2019b),
\begin{equation} \mu_\psi=\alpha \ln \frac{\psi}{1-\psi}-\beta\frac{(1-\hat\phi^2)^2}{2} +\gamma\frac{\hat\phi^2}{2}.\end{equation}
The first term, the free energy of mixing term, favours a uniform surfactant distribution throughout the entire domain and plays the part of diffusion, with the coefficient 
$\alpha$ controlling the magnitude of the diffusive process. The adsorption term (second term) is a negative contribution to the free energy of the system because the accumulation of surfactant at the interface reduces the total energy of the surfactant configuration. The coefficient 
$\beta$ controls the adsorption dynamics. The last term, the bulk-penalty term, is representative of the cost of free surfactant, i.e. of surfactant dissolved in the bulk of the phases rather than adsorbed at the interface, and the coefficient 
$\gamma$ determines the energy cost of surfactant dissolved in the bulk phases. The adsorption term is maximum in magnitude at the interface (
$\hat \phi =0$), indicating a decrease in the energy of the system, while the bulk-penalty term is maximum in the bulk of the phases (
$\hat \phi =\pm 1$) indicating an increase in the energy. The logarithmic formulation of the free energy of mixing term mandates for a non-constant mobility parameter 
$M_\psi =m \psi (1-\psi )$ (Engblom et al. Reference Engblom, Do-Quang, Amberg and Tornberg2013), with 
$m$ being a numerical coefficient controlling the magnitude of the diffusive-like surfactant dynamics. This choice of the mobility parameter ensures the boundedness of the surfactant concentration, 
$\psi \in [0,1]$.
Finally, we couple the volume-of-fluid method for simulating the interfacial dynamics with a phase-field-based method to track the concentration of a soluble surfactant. The volume-of-fluid method guarantees exact mass conservation of each phase and allows for a sharper interface compared with other diffuse-interface methods. We rely on a method to simulate surfactant dynamics that has been successfully adopted in the past to simulate flows with surfactant-laden interfaces (Engblom et al. Reference Engblom, Do-Quang, Amberg and Tornberg2013; Yun et al. Reference Yun, Li and Kim2014; Soligo et al. Reference Soligo, Roccon and Soldati2019a,Reference Soligo, Roccon and Soldatib; Soligo, Roccon & Soldati Reference Soligo, Roccon and Soldati2020a,Reference Soligo, Roccon and Soldatib) using a phase-field method to model interfacial dynamics. In particular, we employ a volume-of-fluid method to simulate the dynamics of the dispersed and carrier phases, and we use a smoothed colour function, 
$\hat \phi$, to couple the interfacial dynamics (based on the volume of fluid) with the surfactant dynamics (based on a phase-field method). The smoothed volume of fluid field accounts for the interfacial dynamics while at the same time providing the diffuse-interface basis onto which the surfactant model is built. We thus combine the aforementioned strengths of the volume-of-fluid method with the advantages of a formulation of the surfactant phase that accounts for adsorption, desorption and diffusion in a thermodynamically consistent framework.
2.1. Computational method
The system of equations is discretized on a uniform Cartesian grid. The computational grid is staggered: pressure, density, viscosity, volume of fluid and surfactant concentration are defined at the cell centres, and the fluid velocities are stored at the cell faces. Spatial derivatives are approximated using a second-order finite difference scheme, and time advancement is performed via a second-order, explicit Adams–Bashforth scheme. A fractional-step method (Kim & Moin Reference Kim and Moin1985) is adopted to advance the mass and momentum conservation equations in time, with the resulting Poisson equation for the pressure solved via a fast pressure solver. The volume of fluid is transported using a directional splitting method combined with an upwind scheme (Ii et al. Reference Ii, Sugiyama, Takeuchi, Takagi, Matsumoto and Xiao2012; Rosti et al. Reference Rosti, De Vita and Brandt2019a). The same scheme is used for the advective term in the surfactant transport equation.
The surfactant is resolved on a refined grid to capture the steep concentration gradients at the interface and to keep a sharp interfacial profile. The smoothing width of the colour function 
$\hat \phi$ should be large enough to accurately discretize the surfactant profile across the interface, and at the same time, it should be small enough to keep the surfactant profile sharp. Using a refined grid thus allows us to capture the modelled interfacial dynamics while maintaining a thin interfacial surfactant layer. The finer computational grid used for the surfactant transport is still a staggered, uniform, Cartesian grid; linear interpolation is used to interpolate variables from/to the standard grid (for the velocity, pressure, density, viscosity and volume of fluid) to/from the fine grid (for the surfactant concentration). The surfactant transport is carried on the fine grid, with the velocity and smoothed volume-of-fluid fields interpolated to the fine grid. Surface tension forces are instead at first computed on the fine grid and then applied to the momentum conservation equation. Tests have been performed with different grid refinement factors: the pressure jump across the interface, the surfactant concentration value at the interface and the total surfactant concentration show minimal changes compared with the reference case (i.e. unitary refinement factor). For the sake of comparison among the different cases, the smoothing width was kept constant in all cases, while in our numerical simulations, the smoothing width is adapted to the refinement factor, thus allowing for smaller values of the smoothing width and for a thinner surfactant interfacial layer.
We use the in-house code Fujin to perform all the numerical simulations presented here. The code has been used and validated in the past on a variety of different flow configurations (Rosti et al. Reference Rosti, Ge, Jain, Dodd and Brandt2019b; Olivieri et al. Reference Olivieri, Brandt, Rosti and Mazzino2020; Brizzolara et al. Reference Brizzolara, Rosti, Olivieri, Brandt, Holzner and Mazzino2021; Cannon et al. Reference Cannon, Izbassarov, Tammisola, Brandt and Rosti2021; Mazzino & Rosti Reference Mazzino and Rosti2021; Rosti & Takagi Reference Rosti and Takagi2021; Abdelgawad, Cannon & Rosti Reference Abdelgawad, Cannon and Rosti2023; Rosti, Perlekar & Mitra Reference Rosti, Perlekar and Mitra2023). Further validation cases are available on the group's website, https://groups.oist.jp/cffu/code. Specific validation tests for the surfactant model and its implementation are reported in Appendix B.
2.2. Computational set-up
We perform direct numerical simulations in a cubic box of size 
$L$ with periodic boundary conditions in all spatial directions. Homogenous isotropic turbulence is sustained using the force 
$\boldsymbol {f_s}$ in (2.1). We use the spectral forcing scheme proposed by Eswaran & Pope (Reference Eswaran and Pope1988), whereby the flow is forced in a shell of Fourier modes 
${2{\rm \pi} /L_{a} \leq |\boldsymbol {k}| \leq 2{\rm \pi} /L_{b}}$, and the force on each mode evolves randomly in time (Uhlenbeck & Ornstein Reference Uhlenbeck and Ornstein1930) with variance 
$\rho \sigma _L^2$ and relaxation time 
$T_L$. Hence, 
${U_L\equiv \sigma _L^{2/3}T_L^{1/3}(L/2{\rm \pi} )^{1/3}}$ is the characteristic velocity scale of the forcing. We set the forcing Reynolds number 
${Re_L\equiv \rho U_L L/(2 {\rm \pi}\mu )=41.6}$ to give a turbulence intensity that is tractable on our numerical grid. We choose the dimensionless relaxation time 
${T_L^*\equiv 2{\rm \pi} T_L U_L/L=2.08}$ to give variations at the time scale of the large eddy turnover time. To prevent droplets from spanning the entire periodic domain, we set the forcing at a length scale smaller than 
$L$ (Mukherjee et al. Reference Mukherjee, Safdari, Shardt, Kenjereš and Van den Akker2019; Crialesi-Esposito et al. Reference Crialesi-Esposito, Chibbaro and Brandt2023). Hence, the minimum and maximum wavelengths of forcing are set to 
$L_{b}=L/3$ and 
$L_{a}=L/2$, respectively.
The computational domain is discretized using an equispaced Cartesian grid with 
$N=500$ grid points in all directions; to better resolve the sharp surfactant gradients at the interface and keep a sharper surfactant profile across the interface, the surfactant transport equation is resolved on a twice-refined grid, with 
$N_\psi =1000$ grid points in all directions. A refinement factor of 2 for the surfactant concentration grid has been selected as it significantly improves the surfactant profile's sharpness at the interface while keeping the computational cost within reasonable limits. With this refinement factor, the computational cost increases by roughly 
$\sim 25\,\%$ and the storage requirements by 
$\sim 116\,\%$.
We report in table 1 the chosen parameters for all cases reported in this paper. We use one single-phase reference case (
$\langle \phi \rangle =0$), three cases with clean droplets (
$\langle \psi \rangle =0$) and three cases with surfactant-laden droplets (
$\langle \psi \rangle =0.1$). The droplet-laden flows were initialized using fluid velocity and pressure from the single-phase reference case once it had reached a statistically steady state. A single spherical droplet of radius 
$R\simeq 0.288L$ (corresponding to 
$\langle \phi \rangle =0.1$) was initialized at the centre of the computational box. Due to the action of the surrounding turbulent flow, the droplet deforms and breaks apart into smaller droplets. For our surfactant-laden cases, the surfactant is initially distributed in the domain following the equilibrium profile with 
$\psi =0.1$ in the bulk phase, computed by zeroing the gradient of the chemical potential. We focus on a highly soluble surfactant and set the coefficients of the surfactant chemical potential to 
$\alpha ={0.0242 u'_0}^2$, 
$\beta ={0.0121 u'_0}^2$ and 
$\gamma ={0.0121 u'_0}^2$, and the numerical coefficient of the mobility parameter to 
$m=0.0307L_b/u'_0$. For the flows with droplets, we fix the reference surface tension 
$\sigma _0$, allowing us to define a reference Weber number 
$We\equiv \rho {u'_0}^2 L_{b}/\sigma _0$ based on the single-phase root-mean-square velocity 
$u'_0$ and the minimum wavelength of the forcing 
$L_{b}$. We select a moderate-strength surfactant with elasticity number 
$\beta _s=5$ and a relatively high surfactant saturation concentration, yielding a low minimum surface tension, 
$\sigma _{min}=0.1\sigma _0$.
List of simulations performed. Here 
$\langle {\cdot }\rangle$ denotes an average over the domain volume and 
$\langle {\cdot }\rangle _I$ denotes an average over the interface. All simulations have been carried out at constant volume fraction 
$\langle \phi \rangle$; an additional reference case (single phase, 
$\langle \phi \rangle =0$) is performed. We investigate four different Weber numbers (
$We$) and two different values of the mean surfactant concentration 
$\langle \psi \rangle$. The measured values of the Kolmogorov length scale 
$\eta$, the Taylor micro-scale 
$\lambda$, Taylor Reynolds number 
$Re_\lambda$, average surface tension 
$\langle \sigma \rangle _I$, effective Weber number 
$We_e$ and Kolmogorov–Hinze diameters 
$d_H$ and 
$d_{H\sigma }$ are reported. The largest and smallest values of each parameter are shown in bold.
To verify the grid independence of the present results, we performed two additional simulations on a more refined grid, 
$N_f=2N=1000$ grid points: a single-phase (same parameters as case SP) case and a surfactant-free multiphase case (same parameters as case C10). Results on the fine grid, 
$N_f=1000$, are reported in the following and summarized here. The scale-by-scale energy budget confirms that the standard grid, 
$N=500$, is sufficient to capture all relevant turbulence scales. For the multiphase case, we do not report any relevant differences between the statistics of the fine grid and the standard grid cases; we do observe the presence of smaller droplets in the fine grid case, as expected due to the increased resolution. Clearly, when adopting a more refined grid in the numerical simulation of multiphase flows, we do expect differences in the morphology of the dispersed phase, as smaller droplets and interfacial features (for instance, ligaments and sheets) can be resolved. This impacts the droplet size distribution and the simulation of interface breaking and merging. Interface breaking is a relatively fast phenomenon that can be well approximated using a continuum formulation (Eggers Reference Eggers1995, Reference Eggers1997; Soligo et al. Reference Soligo, Roccon and Soldati2019a). It has been shown that grid resolution has a minor effect on the simulation of interface breaking (Herrmann Reference Herrmann2011; Lu & Tryggvason Reference Lu and Tryggvason2018, Reference Lu and Tryggvason2019). Interface merging, on the other hand, is dependent on the grid resolution: the final stages of interface merging depend on physics acting at scales as small as the molecular scale (MacKay & Mason Reference MacKay and Mason1963; Aarts, Schmidt & Lekkerkerker Reference Aarts, Schmidt and Lekkerkerker2004; Chen et al. Reference Chen, Kuhl, Tadmor, Lin and Israelachvili2004; Aarts et al. Reference Aarts, Lekkerkerker, Guo, Wegdam and Bonn2005; Kamp, Villwock & Kraume Reference Kamp, Villwock and Kraume2017; Sreehari et al. Reference Sreehari, Borg, Chubynsky, Sprittles and Reese2019). These small scales are not resolved in continuum simulations of multiphase flows and interface merging occurs whenever the interface–interface distance becomes smaller than the grid spacing. Since coalescence occurs based on an artificial length scale, it has been named numerical coalescence. As interface merging involves physics down to the molecular scale, a complete simulation of coalescence cannot be attained solely by grid refinement (Scardovelli & Zaleski Reference Scardovelli and Zaleski1999; Tryggvason et al. Reference Tryggvason, Dabiri, Aboulhasanzadeh and Lu2013; Soligo et al. Reference Soligo, Roccon and Soldati2019a, Reference Soligo, Roccon and Soldati2021). Numerical coalescence is still an open issue and several approaches and models have been put forward to improve the simulation of interface merging; see Soligo et al. (Reference Soligo, Roccon and Soldati2021) for a general review. In our case, we do not employ any models for interface merging, and as shown in figures 3–6 and 8, grid resolution effects are seen only for statistics on the smallest droplets.
3. Results
In table 1 we report integral quantities from all cases studied. Length scales of the turbulent flow are the Kolmogorov scale 
$\eta \equiv (\mu /\rho )^{3/4}\varepsilon ^{1/4}$ and the Taylor micro-scale 
$\lambda \equiv \sqrt {15\mu /(\rho \varepsilon )}u'$, where 
$u'$ is the root-mean-square velocity of the flow and 
$\varepsilon$ is the mean dissipation rate. The Taylor Reynolds number 
$Re_\lambda \equiv \rho u' \lambda / \mu$ is 178 in the single-phase case, and as was previously observed by Rosti et al. (Reference Rosti, Ge, Jain, Dodd and Brandt2019b) and Crialesi-Esposito et al. (Reference Crialesi-Esposito, Rosti, Chibbaro and Brandt2022), 
$Re_\lambda$ increases slightly when droplets are present. The surface tension averaged over the interface 
$\langle \sigma \rangle _I$ is the same as the reference surface tension 
$\sigma _0$ in the cases with clean droplets. However, it is reduced by more than half in the presence of surfactant. This motivates us to define an effective Weber number, 
$We_{e}\equiv \rho u'^2 L_{b}/\langle \sigma \rangle _I$, to better compare the different cases.
The Kolmogorov–Hinze diameter 
$d_H\equiv 0.725 \langle \sigma \rangle _I^{3/5}\rho ^{-3/5} \varepsilon ^{-2/5}$ is an estimate of the diameter of the largest droplet that does not break up. It is made by balancing surface tension with turbulent velocity fluctuations, using an empirical constant of 
$0.725$ (Hinze Reference Hinze1955). We also use a more recent formulation of the Kolmogorov–Hinze diameter from Crialesi-Esposito et al. (Reference Crialesi-Esposito, Chibbaro and Brandt2023); at large scales, droplets predominantly break up and the interface takes energy from the flow (negative work), whereas at smaller scales, droplets predominantly coalesce and the interface returns energy to the flow (positive work). The length scale at which the work done by the interface is zero is defined as 
$d_{H\sigma }$. The two estimates of the Kolmogorov–Hinze diameter are in fairly good agreement in the cases with and without surfactant.
Figure 2(a) shows the kinetic energy spectra of the turbulent flows. In all cases, the most energetic modes are in the range 
$2 \le k L/(2{\rm \pi} )\le 3$, where turbulent forcing is applied. The single-phase case shows the 
$k^{-5/3}$ scaling predicted by Kolmogorov (Reference Kolmogorov1941), persisting for over a decade of wavenumbers. As it has been previously reported by Perlekar (Reference Perlekar2019), Rosti et al. (Reference Rosti, Ge, Jain, Dodd and Brandt2019b) and Crialesi-Esposito et al. (Reference Crialesi-Esposito, Rosti, Chibbaro and Brandt2022), the cases with droplets show a slight reduction in energy at small wavenumbers and an increase at large wavenumbers. We also see from figure 2 that the Kolmogorov–Hinze scale is well within the inertial range. Hence, we can assume that self-similarity applies to the turbulent velocity fluctuations that dictate droplet deformation and breakup, investigated in the following subsection. The turbulent kinetic energy spectrum computed on the fine grid (
$N_f=1000$, shown with empty markers) superposes onto that computed on the standard grid and further extends into the small scales, thus indicating that the standard grid is sufficiently refined to simulate turbulence accurately. Our simulations are in a statistically steady state, and so the fluid kinetic energy contained in each wavenumber is constant in time, and the energy flux through each wavenumber 
$k$ is constant and equal to the energy injection rate 
$\epsilon$. This is expressed by the equation 
${\mathcal {F}(k)+\varPi (k)+\mathcal {D}(k)+\varSigma (k)=\epsilon }$, where the terms on the left-hand side are the energy flux due to forcing, advection, viscous dissipation and surface tension, respectively. We calculate these terms using the method given in the supplementary information of Abdelgawad et al. (Reference Abdelgawad, Cannon and Rosti2023) (see also chapter 6 of Pope Reference Pope2000). Namely, we take a three-dimensional Fourier transform of each term in the Navier–Stokes equation (2.1), multiply by the fluid velocity and integrate over the region bounded by a sphere of radius 
$k$ in wavenumber space. For dissipation, we choose the region inside the sphere; for the other terms, we choose the region outside the sphere. This way, as 
$k \to \infty$, 
$\mathcal {D}=\epsilon$ and 
$\mathcal {F}=\varPi =\varSigma =0$. Figure 2(b) shows the energy balance for the single-phase flow and a droplet-laden flow. The single-phase case shows the canonical Richardson cascade; energy is injected by forcing at the large scales and carried to smaller scales by advection, where it is dissipated by viscosity. In the droplet-laden case, surface tension also carries energy to smaller scales. Minor differences can be appreciated between the standard grid cases (solid lines) and the fine grid cases (semi-transparent lines). For the single-phase case, all the energy is almost completely transferred into dissipation and a plateau in the viscous dissipation is reached at large wavenumbers. A similar result is also observed for the multiphase case, where a fraction of the energy remains stored at small scales in the surface tension term for both grid resolutions.
(a) Turbulent kinetic energy spectrum 
$E$ for cases SP and C10, made dimensionless using the domain size 
$L$ and the root-mean-square velocity 
$u'$ of each case. Vertical dashed lines show the wavenumber of the Kolmogorov–Hinze scale 
$k_H\equiv 2{\rm \pi} /d_H$ for all the cases with droplets. The solid grey line above reports the Kolmogorov scaling 
$k^{-5/3}$. Data from simulations performed on a fine grid, 
$N_f=1000$, are reported with empty markers. (b) Scale-by-scale energy balance for cases SP and C10. Energy flux due to forcing 
$\mathcal {F}$, viscous dissipation 
$\mathcal {D}$, advection 
$\varPi$ and surface tension 
$\varSigma$ are plotted using dot-dashed, dotted, dashed and solid lines, respectively; semi-transparent lines identify the data from the simulations performed on a fine grid. Vertical dashed lines mark 
$k_{H\sigma }\equiv 2{\rm \pi} /d_{H\sigma }$, the wavenumber at which 
$\varSigma$ is maximum for every droplet case.
3.1. Droplet statistics
The inset of figure 3(a) shows the average number of droplets 
$\mathcal {N}$ in our simulations. To identify and count each droplet, we use a stack-based, six-way flood fill on the computational cells characterized by 
${\phi \ge 0.5}$; the algorithm is a direct extension to three-dimensional space of the two-dimensional four-way flood fill algorithm (Newman & Sproull Reference Newman and Sproull1979). The number of droplets has been counted over several instantaneous snapshots and averaged in time once the simulation has reached a statistically steady state, i.e. once the Taylor–Reynolds number and the number of droplets fluctuate about a constant mean value. We note that clean and surfactant-laden cases at similar values of the effective Weber number, i.e. S05 and C10, S10 and C20, S20 and C40, have approximately the same average number of droplets. This result suggests that the effect of surfactant on the dispersed phase manifests mainly as an average surface tension reduction, with negligible effects from Marangoni stresses. As previously found by Rosti et al. (Reference Rosti, Ge, Jain, Dodd and Brandt2019b) in shear turbulence, we see that the number of droplets is proportional to the Weber number, in our case with a factor 
$\mathcal {N}=434We_e$.
(a) The PDF of the droplet diameters 
$d$, with the dashed and solid lines showing scalings previously found in the coalescence and breakup regimes (Deane & Stokes Reference Deane and Stokes2002), respectively; data from the simulation performed on a fine grid, 
$N_f=1000$, are reported with empty markers. The inset reports the total number of droplets 
$\mathcal {N}$ in each case, with error bars showing the root-mean-square variation in time. The dotted line is a fit 
${\mathcal {N}=434We_e}$. In (b) we calculate the mean and standard deviation of the surfactant concentration 
$\psi$ on the interface of each drop; we then average these values over the ensemble of drops of the same size. The right-hand axis shows the normalized interfacial surface tension resulting from the presence of surfactant.
The average number of droplets 
$\mathcal {N}$ is of the order 
$10^4$. This large sample size allows us to make accurate statistics of the droplets, even when binned by their equivalent diameter 
$d$. We define the equivalent diameter of a droplet
\begin{equation} d\equiv(6V/{\rm \pi})^{1/3},\end{equation}
as the diameter of a sphere with volume 
$V$, where 
$V$ is the volume of said droplet. The droplet size distribution for all cases is shown in figure 3(a). We observe a collapse of all the curves when the droplet size 
$d$ is normalized by the Kolmogorov–Hinze scale in each case. Note that the Kolmogorov–Hinze scale separates two different regimes: the coalescence-dominated regime for 
$d/d_H<1$ and the breakage-dominated regime for 
$d/d_H>1$. The Kolmogorov–Hinze scale (Hinze Reference Hinze1955) is defined as the largest droplet that resists breakage due to turbulent fluctuations. The coalescence-dominated regime characterizes droplets smaller than the Kolmogorov–Hinze scale; breakage is highly unlikely for these droplets, which instead are in a state of constant coalescence. On the other hand, droplets larger than the Kolmogorov–Hinze scale are prone to breaking apart. The droplet size distribution shows a clear power-law behaviour in the breakage- and coalescence-dominated regimes. Considering the rate of breakage, exponents for the two regimes have been obtained: 
$-3/2$ for the coalescence-dominated regime (Riviere et al. Reference Rivière, Ruth, Mostert, Deike and Perrard2022) and 
$-10/3$ for the breakage-dominated regime (Garrett et al. Reference Garrett, Li and Farmer2000). Deane & Stokes (Reference Deane and Stokes2002) measured the bubble size distribution in breaking waves and found good agreement between the experimental measurements and the analytical scalings. Several previous computational works confirmed the same 
$-10/3$ power-law exponent in the breakage-dominated regime (Skartlien et al. Reference Skartlien, Sollum and Schumann2013; Deike et al. Reference Deike, Melville and Popinet2016; Mukherjee et al. Reference Mukherjee, Safdari, Shardt, Kenjereš and Van den Akker2019; Rosti et al. Reference Rosti, Ge, Jain, Dodd and Brandt2019b; Soligo et al. Reference Soligo, Roccon and Soldati2019a; Crialesi-Esposito et al. Reference Crialesi-Esposito, Chibbaro and Brandt2023), while fewer works have captured the 
$-3/2$ power-law exponent for the coalescence-dominated regime (Rivière et al. Reference Rivière, Mostert, Perrard and Deike2021; Crialesi-Esposito et al. Reference Crialesi-Esposito, Chibbaro and Brandt2023). Our results in figure 3(a) show a size distribution that is compatible with the breakage-dominated regime scaling, 
$-10/3$. However, the available data in the breakage-dominated regime spans at most one decade; hence inferring an accurate power-law scaling here is challenging. For the low 
$We_e$ cases, we see some deviation from the 
$-3/2$ scaling just below the Kolmogorov–Hinze scale. However, further into the coalescence-dominated regime, around 
$d/d_H\approx 10^{-1}$, all cases follow the 
$-3/2$ scaling. Furthermore, we show that the macroscopic effect of surfactant on the droplet size distribution is well captured by considering a lower, average surface tension value when computing the Kolmogorov–Hinze scale. A similar result was previously obtained for surfactant-laden flows (Skartlien et al. Reference Skartlien, Sollum and Schumann2013; Soligo et al. Reference Soligo, Roccon and Soldati2019a), although for the breakage-dominated regime alone. Here, we extend the result to droplets smaller than the Kolmogorov–Hinze scale. We also report the droplet size distribution data for the case C10 on a more refined grid (empty markers); a good agreement is observed with data from the same case on the standard grid. Clearly, on the more refined grid, we are able to capture smaller droplets; hence, the droplet size distribution extends slightly further to smaller diameters.
Figure 3(b) shows the dependence of surfactant concentration 
$\psi$ on the equivalent droplet diameter 
$d$. For these values and error bars, the average surfactant concentration and its standard deviation were computed at the interface of each droplet, and then averaged over all the droplets of size 
$d$. This way, the error bars capture not the variation of 
$\psi$ between droplets, but the average variation on each droplet, which governs Marangoni stresses. We see that the average surfactant concentration at the interface is higher than the initial concentration in the bulk phase (
$\psi =0.1$), as from (2.6) it is energetically favourable for the surfactant to assemble on the interface. We observe a trend in the average surface tension value at the interface for increasing values of the Weber number: as the total amount of interfacial area increases (the total number of droplets is roughly proportional to the Weber number), the average surfactant concentration at the interface reduces. For the surfactant parameters considered in this study, there is little dependence of the mean surfactant concentration at the interface on the size of the droplets. The average surfactant concentration for all droplet sizes is about 
$\psi \approx 0.115$, with a slightly lower value for the smallest droplets. This results in an average reduction of the surface tension to approximately 40 % of its clean value. We note that a surface tension reduction of around one-half is typical of real-world surfactant-laden interfaces, such as Tween 80 and NaCl in water (Qazi et al. Reference Qazi, Schlegel, Backus, Bonn, Bonn and Shahidzadeh2020). Error bars show the mean standard deviation of surfactant concentration on a droplet with equivalent diameter 
$d$. The standard deviation on the droplets has instead a mild dependence on the characteristic size of the droplet, showing about a twofold increase between the smallest and the largest droplets (i.e. over a 
$\sim 100\times$ increase in the droplet equivalent diameter). The variation in surface tension is approximately 15 % of the mean surfactant concentration at the interface, corresponding to about 8 % change in surface tension on each droplet.
Next, we look at the shape of the droplets, which provides an indication of the competition between surface tension and turbulence: a droplet in a quiescent fluid takes a spherical shape and any deviation from this shape is to be attributed to the flow. We evaluate the shape of the droplets using two dimensionless parameters: the aspect ratio 
$\sqrt {I_1/I_3}$ and the sphericity 
$d/d_A$.
The aspect ratio is computed as the ratio between the smallest and the largest eigenvalues of the moment of inertia tensor, respectively 
$I_1$ and 
$I_3$, with 
$I_1 \le I_2 \le I_3$. The moment of inertia tensor is computed as
\begin{equation} \boldsymbol{\mathsf{I}}\equiv \int_V \rho (\|\boldsymbol{r}\|^2 \boldsymbol{\mathsf{I}}-\boldsymbol{r}\otimes \boldsymbol{r}) \, {\rm d} V, \end{equation}
where 
$\boldsymbol {r}$ is the vector from the droplet centre of mass to a point inside the droplet and 
$V$ is the volume of the considered droplet. The aspect ratio uses a similar definition of the bubble deformation parameter defined in Bunner & Tryggvason (Reference Bunner and Tryggvason2003), although in their original work, the bubble deformation was defined as the square root of the largest over the smallest eigenvalue and a slightly different formulation of the moment of inertia tensor was used. Here, we choose to define the aspect ratio as the inverse of the bubble deformation found in Bunner & Tryggvason (Reference Bunner and Tryggvason2003), such that the values of the aspect ratio are bounded between 0 (e.g. infinitely long and thin filament or sheet) and 1 (e.g. sphere or cube). To exclude grid resolution effects when calculating aspect ratios, we consider only droplets with a volume of more than 100 grid cells.
The sphericity is defined as the ratio of the volume-equivalent diameter over the surface-equivalent diameter; 
$d/d_A$. We use a different definition of sphericity from that found in the literature (Wadell Reference Wadell1935), such that it is bounded between 0 and 1. The surface-equivalent diameter is defined as the diameter of a sphere having the same surface area 
$A$ of the considered droplet: 
$d_A\equiv \sqrt {A/{\rm \pi} }$. Sphericity is equal to unity for a sphere (the shape with minimal surface area for a given volume) and reduces as the droplet deforms from the spherical shape. The area 
$A$ of each droplet is computed by counting the number of computational cells crossed by the interface and projecting a face of the computational cell onto the local normal 
$\boldsymbol {n}$ to the interface.
We choose these two parameters as they provide very different information, as illustrated in figure 4(a), where we compute the aspect ratio and sphericity for three sample droplets. The aspect ratio is very sensitive to droplet-scale deformations, for instance, stretching along one of the axes, but is relatively unchanged by small-scale perturbations at the interface. Conversely, the sphericity is less sensitive to large, droplet-scale deformations but is very effective in revealing small-scale perturbations of the interface. Qualitatively, the aspect ratio estimates the shape of a box that bounds the droplet. Meanwhile, the sphericity measures the total area of the droplet interface. This is reflected in figure 4(a): taking the spheroidal droplet as a reference (
$\sqrt {I_1/I_3}\approx 1$ and 
$d/d_A\approx 1$), the ellipsoidal droplet shows a negligible change in the sphericity value and a 
$\sim 30\,\%$ reduction in the aspect ratio. The bulgy droplet, being rather compact, has a value of the aspect ratio relatively close to that of the spheroidal droplet (about 10 % difference) but a much smaller value of sphericity (
$\sim 30\,\%$ smaller).
Mean droplet deformation at each equivalent diameter 
$d$ normalized by the Kolmogorov–Hinze scale 
$d_H$. We plot two measures of deformation; (b) the aspect ratio of the droplets 
$\sqrt {I_1/I_3}$ and (c) the sphericity 
$d/d_A$. On panel (b) we shade the region 
$0.5< d/d_H<2$ where the aspect ratio is observed to decrease with equivalent diameter. In both panels (b) and (c) we report data from the simulation performed on a fine grid, 
$N_f=1000$, with empty markers. Panel (a) above shows the values of deformation for (from left to right) a spheroid, an ellipsoid and a bulgy droplet.
We report the aspect ratio 
$\sqrt {I_1/I_3}$ in figure 4(b) as a function of the droplet size normalized by the Kolmogorov–Hinze scale. We observe that making the droplet size dimensionless using the Kolmogorov–Hinze scale yields a collapse of the aspect ratio for all cases onto a single curve. This suggests that even in the presence of surfactants, Hinze's (Reference Hinze1955) assumptions hold, and droplet deformation is a universal function of 
$d/d_H$. We observe three different regimes for the aspect ratio, which can be distinguished based on the value of 
$d/d_H$. The aspect ratio of small droplets, up to approximatively 
$d=0.5d_H$, is roughly constant and close to unity, indicating that the droplets are spherical or only slightly elongated (i.e. compact shape). At this point, no information can yet be inferred on local, small-scale deformations of the interface. Droplets smaller than the Kolmogorov–Hinze scale are characterized by surface tension forces dominating over turbulent fluctuations. The second regime, observed for 
$0.5d_H \lesssim d \lesssim 2 d_H$, is characterized by a sharp reduction in the value of the aspect ratio (down to 
$\sqrt {I_1/I_3}\approx 0.5$, a similar value of an elongated ellipsoid): droplets become more elongated with an overall deformation that increases with their size. This result is coherent with the definition of the Kolmogorov–Hinze scale as the size of the largest (on average) non-breaking droplets: at about the Kolmogorov–Hinze scale, droplets start to deform significantly. Statistics on the sphericity will provide further information on the small-scale deformation within this regime and will be discussed in the following paragraph. Lastly, in larger droplets, 
$d>2d_H$, there is a recovery of the droplet aspect ratio: droplets partially recover from their deformed state, with values of the aspect ratio nearing 
$\sqrt {I_1/I_3}\approx 0.7$. This result indicates that droplets become less deformed overall, and some rotational symmetry is restored. At first, this may be a counter-intuitive result, but sphericity will help to explain this interesting behaviour. We include also data from the refined grid simulation (empty markers) and we report that it closely follows the standard grid data.
We report the sphericity for all cases in figure 4(c). The droplet size is made non-dimensional using the Kolmogorov–Hinze scale for each case. We observe that data for all cases collapse on a single curve when scaled by the Kolmogorov–Hinze scale, further confirming the validity of the Kolmogorov–Hinze scale as a fundamental length scale. We observe two very distinct regimes, separated by the Kolmogorov–Hinze scale. At scales smaller than the Kolmogorov–Hinze scale, droplets have a close-to-unity and slightly decreasing sphericity, which sharply decreases above the Kolmogorov–Hinze scale. Note that, for the smallest droplets, we have values of sphericity larger than one; we would like to remark that these values are not admissible and are due to inaccuracies in the computation of interface normals and area when the droplets are only a few grid cells in volume. The sphericity for droplets smaller than the Kolmogorov–Hinze scale is close to unity and decreases for increasing droplet sizes; this information, coupled with the results from the droplet aspect ratio, indicates that droplets much smaller than Kolmogorov–Hinze scale have a spheroidal shape with limited elongation and almost no small-scale perturbations of the interface. Kolmogorov–Hinze-scale-sized droplets (but still smaller than the Kolmogorov–Hinze scale) show a substantial reduction in the aspect ratio and only a minor decrease in the sphericity: the shape of these droplets is similar to an ellipsoid, as the droplet is stretched (low aspect ratio) but the sphericity is still close to unity (indicating the absence of relevant perturbations of the interface). Conversely, droplets slightly larger than the Kolmogorov–Hinze scale show a reduced aspect ratio and sphericity: these droplets are not only strongly elongated (low aspect ratio), but small-scale perturbations of the interfaces (small humps and dimples) start forming (low sphericity). The trend in sphericity is kept also for droplets much larger than the Kolmogorov–Hinze scale; these droplets show a recovery of the aspect ratio, indicating either the formation of bulgy droplets (see figure 4a) or of convoluted filaments. Both of these shapes are coherent with the two deformation parameters we investigated, i.e. relatively low aspect ratio and low sphericity. Data from the refined grid simulation (empty markers) collapses well onto the standard grid data, confirming grid independence of the results. Due to the higher grid resolution, the fine grid case is able to capture smaller droplets compared with the standard case, however, the general trend is kept.
To distinguish among the two possible shapes, bulgy droplets versus convoluted filaments, we compute the average radius of curvature 
$R$ at the interface of each droplet, reported in figure 5. To compute the radius of curvature 
$R$, we first compute the mean curvature 
$\kappa$ using the divergence of the normal 
$n$ to the interface:
\begin{equation} \kappa=\boldsymbol{\nabla}\boldsymbol{\cdot}\boldsymbol{n}. \end{equation}
We average the curvature 
$\kappa$ over the droplet interface and define the radius of curvature as the inverse of the mean curvature, 
$R\equiv 1/\kappa$. Note that, in (3.3) we do not use a minus sign in the definition of the curvature, as we have an outward-pointing normal 
$\boldsymbol {n}$ and we choose to assign positive values of curvature to convex surfaces, i.e. the curvature is positive if the surface curves away from the normal. Also, the normal is ill-defined for small droplets, so we only calculate the mean curvature of droplets with a volume greater than 10 computational cells. For a surface in three dimensions, such as our interface, the mean curvature is equal to the sum of the two principal curvatures, 
$\kappa =\kappa _1+\kappa _2$. For a sphere, both principal curvatures are equal to the inverse of the sphere's radius; hence, the mean curvature is twice this value and the radius of curvature is 
$R=d/4$. Figure 5 shows that droplets smaller than the Kolmogorov–Hinze scale follow this scaling (grey dashed line). This result further confirms the rather regular shape (spheroid- or ellipsoid-like) of the small droplets. Above the Kolmogorov–Hinze scale instead, we observe a departure from 
$d/4$: the radius of curvature becomes constant, approximately equal to half the Kolmogorov–Hinze scale, and independent of drop size. The good agreement obtained with data from the fine grid simulation (empty markers) indicates that the grid resolution is more than sufficient to capture the local shape of the interface.
The mean radius of curvature 
$R$ of the surfaces of droplets binned by their diameter 
$d$. The dashed line shows 
$R=d/4$, the radius of curvature of a sphere with diameter 
$d$; data from the simulation performed on a fine grid, 
$N_f=1000$, are reported with empty markers. Inset: the same values, plotted in terms of curvature 
$\kappa =1/R$. Error bars show the standard deviation of curvature on the droplets.
To understand this behaviour, we consider a cylinder of radius 
$d_H/2$. On the curved surface of the cylinder, the two principal curvatures are 
$\kappa _1=0$ and 
$\kappa _2=2/d_H$. Neglecting the two flat ends, the mean curvature of the cylinder is thus 
$\kappa =2/d_H$, i.e. the radius of curvature of the cylinder is equal to half the Kolmogorov–Hinze scale, and is independent of its length. This result, combined with the information obtained from the deformation parameters in figure 4, shows that, above the Kolmogorov–Hinze scale droplets take the shape of filaments with a diameter equal to the Kolmogorov–Hinze scale.
The inset of figure 5 shows the mean curvature 
$\kappa$ of the droplets as a function of their equivalent diameter 
$d$. Error bars indicate the standard deviation of the curvature, which was calculated at the interface of each droplet and averaged over the ensemble of droplets of similar size 
$d$. The standard deviation of the curvature of an interface is a measure of its corrugation, and can be used to quantify Plateau–Rayleigh instabilities that lead to droplet breakup (Rayleigh Reference Rayleigh1878; Villermaux & Bossa Reference Villermaux and Bossa2009; Kooij et al. Reference Kooij, Sijs, Denn, Villermaux and Bonn2018). From the inset of figure 5, we see that droplets of all sizes show corrugation, and the standard deviation of 
$\kappa$ is comparable to its average value. Above the Kolmogorov–Hinze scale, we report an increased probability of negative values of the mean curvature, indicative of dimples and saddle points in the surface of the droplet, which, for example, can be due to impinging jets or large Plateau–Rayleigh instabilities on the interface.
To test our hypotheses on the droplet shapes above and below the Kolmogorov–Hinze scale, we consider their surface area. Figure 6 shows the mean surface area of the droplets at each equivalent diameter 
$d$. As we previously noted, the droplets smaller than the Kolmogorov–Hinze scale are almost spherical in shape, and hence, we see that their surface area grows quadratically with their characteristic size. Above the Kolmogorov–Hinze scale, the droplet surface areas instead grow as the cube of their sizes. As a first approximation, we can think of these droplets as long cylinders with the same diameter 
$d_H$ and different lengths 
$l$. Ignoring the two end faces, the surface area of such a cylinder is 
$A_f={\rm \pi} d_H l$ and the volume is 
$V_f = {\rm \pi}d_H^2 l/4$. Substituting the volume for the equivalent diameter (3.1) we get an expression for the length:
\begin{equation} l = \frac{2d^3}{3d_H^2}. \end{equation}
Plugging this into the formula for the cylinder surface area, we obtain 
$A_f=2{\rm \pi} d^3/(3d_H)$, showing that a filament with variable length 
$l$ and constant diameter 
$d_H$ has an interfacial surface area that grows as the cube of its volume-equivalent diameter 
$d$. The interfacial surface area of the droplets above the Kolmogorov–Hinze scale closely follows 
$A_f$, as reported in figure 6. The two scalings are also confirmed by simulation on a more refined grid (empty markers).
Dependence of droplets’ surface area 
$A$ on their equivalent diameter 
$d$. The dashed line shows the surface area of a sphere with diameter 
$d$ and the solid line shows the surface area of a filament with diameter 
$d_H$. Data from the simulation performed on a fine grid, 
$N_f=1000$, are reported with empty markers.
From our observations of curvature and surface area, it appears that droplets above the Kolmogorov–Hinze scale are filamentous in shape. However, their aspect ratio shows the filaments cannot be straight, since this would produce a monotonic reduction in 
$\sqrt {I_1/I_3}$ with 
$d$, which we do not see in figure 4(a). Hence, a natural question is: do the filaments form loops, or are they simply connected? figure 7 shows two droplets from case S20, we see that both droplets are made up of convoluted filaments, and in many places, the filaments do in fact form loops. These complex shapes are typical of large and deformable viscous droplets in both liquid–liquid mixtures (Collins & Knudsen Reference Collins and Knudsen1970; Andersson & Andersson Reference Andersson and Andersson2006b; Li et al. Reference Li, Miller, Wang, Koley and Katz2017; Xue & Katz Reference Xue and Katz2019) and gas–liquid mixtures (Hinze Reference Hinze1955; Villermaux Reference Villermaux2007; Blenkinsopp & Chaplin Reference Blenkinsopp and Chaplin2010; Dumouchel & Blaisot Reference Dumouchel and Blaisot2014; Jain et al. Reference Jain, Prakash, Tomar and Ravikrishna2015; Kooij et al. Reference Kooij, Sijs, Denn, Villermaux and Bonn2018; Chéron et al. Reference Chéron, Brändle De Motta, Blaisot and Menard2022; Qi et al. Reference Qi, Tan, Corbitt, Urbanik, Salibindla and Ni2022; Kant et al. Reference Kant, Pairetti, Saade, Popinet, Zaleski and Lohse2023). Small droplets (smaller than the Kolmogorov–Hinze scale) have instead a more spherical shape, as indicated in figure 4. To answer the question quantitatively, we measure the topology of the interface of each droplet using its Euler characteristic 
$\chi$, which obeys the equation
\begin{equation} 1-\chi/2=h-v,\end{equation}
where 
$h$ is the number of handles and 
$v$ is the number of voids in the drop (see Appendix A for a derivation of this equation). In the insets of figure 8(a) we show renders of example droplets with one and two handles. Analogously to the handle on a teacup, a handle is a loop of the dispersed phase that extends through the carrier phase. The renders in the insets of figure 8(b) show example droplets with one and two voids. A void is a region of the carrier phase entirely enclosed by the dispersed phase. Note that similarly to droplet breakup and coalescence, a change in the number of handles or voids necessitates a merging or splitting of interfaces. Using a method similar to Mendoza et al.'s (Reference Mendoza, Thornton, Savin and Voorhees2006) characterization of dendritic metal samples, we measure the Euler characteristic using simplicial homology, that is, by dividing the interface into simple polygons and counting the number of nodes 
$n$, edges 
$e$ and faces 
$f$ on the interface of each droplet. The Euler characteristic of the interface is then given by the Poincaré formula (Massey Reference Massey1997, p. 26)
\begin{equation} \chi=n-e+f.\end{equation}
Our volume-of-fluid field 
$\phi$ is defined on a cubic grid and we can define the interface as the boundaries between cells where 
$\phi -0.5$ changes sign. Hence, our interface is already divided into square faces and 
$\chi$ can be calculated by counting these faces and their edges and nodes. We also count the number of voids 
$v$ on each droplet. Voids are counted by rerunning the flood fill algorithm, looking for contiguous cells where 
$\phi <0.5$; to do so, we use a stack-based 26-way flood fill, which is an extension to three dimensions of the two-dimensional eight-way flood fill (Newman & Sproull Reference Newman and Sproull1979). Knowing 
$\chi$ and 
$v$ for each drop, we can use (3.5) to obtain the number of handles 
$h$.
Visualization of two droplets extracted from case S20. The left droplet has 56 handles and six voids. The right droplet has two handles and zero voids. The length of the black bar is the Kolmogorov–Hinze scale.
Figure 8(a) shows how the mean number of handles per droplet depends on the droplet size. We see that the largest droplets are very self-connected, having of the order of 
$10^2$ handles. Furthermore, the number of handles in each case collapses to a single line when the equivalent diameter is normalized by the Kolmogorov–Hinze scale 
$d_H$. This universality occurs because surface tension is constantly acting to destroy the handles, so the higher the surface tension, the shorter the lifetime of the handle. The fit 
$h=0.04(d/d_H)^3$ gives the number of handles as a function of the droplet volume. Using our earlier result that the large droplets are filaments with lengths given by (3.4), we can convert this expression into the number of handles per unit length of the filament: 
$h/l=0.06/d_H$. Figure 8(b) shows us, on the other hand, that the number of voids is independent of surface tension, which is reasonable as a void inside a droplet experiences no net surface tension force. The number of voids simply scales with the volume of the droplet. Again we can obtain a fit, 
$v=500(d/L)^3$; substituting the equivalent diameter for the droplet volume 
$V$ (3.1) we find there are roughly 
$v/V\approx 950/L^3$ voids per unit droplet volume in all cases. We suspect the void concentration depends on the droplet coalescence rate and turbulence intensity (
$Re_\lambda$). However, we leave a proper investigation of this dependence to future works.
(a) Number of handles 
$h$ per droplet. The grey line shows the fit 
$h=0.04(d/d_H)^3$. Renders show example droplets with one and two handles (
$\chi =0$ and 
$\chi =-2$, respectively). (b) Number of voids 
$v$ per droplet. The grey line shows the fit 
$v=500(d/L)^3$. We mark the 
$x$ axis to show the Kolmogorov–Hinze scale of each case. We report data from the simulation performed on a fine grid, 
$N_f=1000$, using empty markers. Renders show example droplets with one and two voids (
$\chi =4$ and 
$\chi =6$, respectively). The renders are obtained from droplets we artificially generated for the sole purpose of demonstrating handles and voids.
Data from the fine grid simulation is reported with empty markers for both handles and voids, and a good agreement with data computed on the standard grid is observed. When counting the number of voids in the fine grid case, we excluded voids with a volume smaller than eight computational cells, (i.e. one computational cell on the standard grid). This method allows for a fair comparison of the number of voids in the fine grid and the standard grid cases.
3.2. Flow statistics
We characterize the effects induced by the presence of clean and surfactant-laden interfaces on the local flow statistics using the flow topology parameter (Perry & Chong Reference Perry and Chong1987), which compares the local flow to three different base flows: purely rotational, pure shear and purely extensional flow. The flow topology parameter 
$Q$ is a combination of the rate-of-strain tensor 
$\boldsymbol{\mathsf{S}}\equiv (\boldsymbol {\nabla } \boldsymbol {u}+\boldsymbol {\nabla } \boldsymbol {u}^T)/2$ and the rate-of-rotation tensor 
$\boldsymbol{\varOmega} \equiv (\boldsymbol {\nabla } \boldsymbol {u}-\boldsymbol {\nabla } \boldsymbol {u}^T)/2$, where 
$\boldsymbol {\nabla }\boldsymbol {u}$ is the velocity gradient tensor. The quantities 
$S^2$ and 
$\varOmega ^2$ are defined as 
$S^2=\boldsymbol{\mathsf{S}}:\boldsymbol{\mathsf{S}}$ and 
$\boldsymbol{\varOmega} ^2=\boldsymbol{\varOmega} :\boldsymbol{\varOmega}$, where ‘
$:$’ identifies the dyadic (double-dot) product. The flow topology parameter categorizes the local flow based on three types of base flow: purely rotational flow (
$S^2=0$ and 
$\varOmega ^2 \ne 0$), purely extensional flow (
$S^2\ne 0$ and 
$\varOmega ^2 = 0$), and pure shear (
$S^2=\varOmega ^2$).
\begin{equation} Q=\frac{S^2-\varOmega^2}{S^2+\varOmega^2}= \begin{cases} -1 & \text{purely rotational,}\\ 0 & \text{pure shear,}\\ +1 & \text{purely extensional.} \end{cases} \end{equation} We compare clean and surfactant-laden cases at similar effective Weber numbers in figure 9, together with the single-phase case for reference. The flow topology parameter is computed in three distinct regions: inside the droplets (
$\phi >0.5$), in the carrier phase (
$\phi <0.5$) and at the interface. This way, we can separate the contribution from the different regions of the flow (Dodd & Jofre Reference Dodd and Jofre2019; Rosti et al. Reference Rosti, De Vita and Brandt2019a; Soligo et al. Reference Soligo, Roccon and Soldati2020b), and investigate the effect of Marangoni stresses at the interface for the surfactant-laden cases and of flow confinement on the flowing condition inside the droplets and at the interface.
Flow topology parameter sampled in different regions of the domain: inside the droplets (solid line), in the carrier phase (dashed line) and at the interface (dotted line). The single-phase case is reported for reference (dash-dotted line). Each panel refers to an approximate value of the effective Weber: (a) 
$We_e\approx 10$ (cases C10 and S05), (b) 
$We_e\approx 20$ (cases C20 and S10) and (c) 
$We_e\approx 40$ (cases C40 and S20).
We first consider the flow topology parameter in the carrier phase. The relatively low volume fraction of the dispersed phase reduces the overall impact of the presence of the interface on the outer flow. The flow topology parameter for all clean and surfactant-laden cases well collapses onto the single-phase line, indicating indeed that the presence of a deformable interface and of Marangoni stresses does not introduce any significant modification of the outer flow at the relatively low volume fraction considered. A similar result was reported by Rosti et al. (Reference Rosti, De Vita and Brandt2019a) for clean droplets and for volume fractions of the dispersed phase up to 30 %. In general, the flow topology shows a predominance of a combination of pure shear and extensional flow, with a limited rotational contribution.
The flow topology computed inside the droplets highlights the effect of flow confinement. Surface tension has competing effects: while on the one hand, a large value of surface tension is more effective in decoupling internal and external flow, a low value of surface tension produces many small-size droplets that increase the flow-confinement effect. We observe a reduction in extensional flow and an increase in pure shear for increasing values of the Weber number, which has been attributed to the confinement effect of small droplets (Soligo et al. Reference Soligo, Roccon and Soldati2020b). We do not report changes in the rotational component, which can be attributed to the lack of large-scale coalescence events (Rosti et al. Reference Rosti, De Vita and Brandt2019a). When comparing cases at similar effective Weber numbers, we observe that the addition of surfactant has a similar effect to a reduction in the surface tension value: a reduction in the extensional flow and an increase in pure shear. The difference in the flow topology between clean and surfactant-free droplets reduces as the effective Weber number increases.
It is interesting to note that, at the interface, the addition of surfactant instead has an opposite effect: the presence of surfactant, especially at low values of the effective Weber number, leads to an increase of extensional flow and a decrease in rotational flow – pure shear is unchanged. This difference results from the action of Marangoni stresses, and is indeed more apparent at low values of the effective Weber number, i.e. at high values of surface tension. We observe a shift towards extensional flow for increasing values of the effective Weber number, confirming previous findings (Soligo et al. Reference Soligo, Roccon and Soldati2020b).
We now focus on flow statistics at the interface to better understand the local surface tension effects induced by clean and surfactant-laden interfaces. We compute the alignment of vorticity at the interface of the droplets: to quantify the direction of vorticity, we use the cosine of the angle 
$\theta$; see the sketch in figure 10(a). This quantity is found by taking the scalar product between the interface normal 
$\boldsymbol {n}$ (outward-pointing, unit-length normal) and the unit-length vorticity vector, 
$\boldsymbol {\omega }/\|\boldsymbol \omega \|$. The probability density function (PDF) of the vorticity-interface alignment is reported in figure 10(a) for all the cases we simulated. Similarly to what was found by Mukherjee et al. (Reference Mukherjee, Safdari, Shardt, Kenjereš and Van den Akker2019), the vorticity is mostly orthogonal to the interface normal. At very large values of the surface tension, i.e. 
$We\to 0$, the interface between the droplet and the carrier phases acts similarly to a slip wall: the high surface tension makes the interface close to undeformable but imposes no condition on the tangential component of the flow. We thus expect that at high values of surface tension the vorticity vector is tangential to the interface, i.e. 
$\boldsymbol {n}\boldsymbol {{\cdot }}\boldsymbol \omega /\|\boldsymbol \omega \|=0$. As the interface becomes more deformable (corresponding to a higher Weber number), this condition is relaxed and the PDF widens. We observe that the PDF remains symmetric for all cases; this is an expected result as positive and negative values of the alignment correspond to flow structures at the interface rotating in the anti-clockwise and clockwise directions, respectively, and the two are equally probable. The alignment of vorticity at the interface also highlights the effect of surfactant and, in particular, of Marangoni stresses tangential to the interface: the cases at a similar effective Weber number (namely 
$We_e\approx 10$ C10 and S05; 
$We_e\approx 20$ C20 and S10; 
$We_e\approx 40$ C40 and S20) show notably different distributions, suggesting that in terms of the local flow around the droplets, the effect of surfactant cannot be approximated as only an average surface tension reduction. We indeed note a decoupling among the various cases at approximately the same effective Weber number, with cases C40 (
$We_e=39.3$) and S10 (
$We_e=25.8$) showing the very same distribution in vorticity alignment. The presence of Marangoni stresses promotes the formation of flows tangential to the interface, whose gradients contribute to the vorticity component normal to the interface.
Alignment between unit-length normal to the interface 
$\boldsymbol {n}$ and (a) vorticity at the interface, (b) velocity at the interface. The inset in panel (a) shows how the alignment has been computed: for a generic vector 
$\boldsymbol {a}$, we define the alignment as 
$\boldsymbol {n}\boldsymbol {{\cdot }} \boldsymbol {a}/\|\boldsymbol {a}\|$, which is equal to the cosine of the angle 
$\theta$ in between the two vectors.
The alignment of the fluid velocity at the interface further corroborates the role of Marangoni stresses in modifying the local flow velocity at the interface; figure 10(b) shows the PDF of the scalar product between the unit-length velocity vector and the normal to the interface. The surfactant-laden cases show a more peaked distribution at 
$\boldsymbol {n}\boldsymbol {{\cdot }} \boldsymbol {u}/\|\boldsymbol {u}\|\approx 0.1$, corresponding to a local fluid velocity almost tangential to the interface, with a small outward (i.e. from the droplet phase towards the carrier phase) component. This velocity alignment is clearer for the cases at low Weber number: the magnitude of Marangoni stresses directly depends on the local surface tension, hence, the cases at high Weber number are characterized by weaker Marangoni stresses. Indeed, the case S20 shows similar velocity alignment to the surfactant-free cases, being that the Marangoni stresses are weaker compared with the other surfactant-laden cases. For the clean cases, we observe two separate peaks in the distribution, one at 
$\boldsymbol {n}\boldsymbol {{\cdot }} \boldsymbol {u}/\|\boldsymbol {u}\| = -1$ and one at 
$\boldsymbol {n}\boldsymbol {{\cdot }} \boldsymbol {u}/\|\boldsymbol {u}\|\approx 0.3$. The former corresponds to an inward flow perpendicular to the interface and the latter to a fluid velocity mainly tangential to the interface, although with a larger normal component compared with the surfactant-laden cases. We attribute this reduction in the probability of having flow tangential to the interface to the absence of Marangoni stresses for the surfactant-free cases. A similar distribution is also achieved by case S20, which is characterized by weak Marangoni stresses (due to the low reference surface tension value), further highlighting the role of Marangoni stresses.
So far we have only considered the angle between the flow velocity at the interface and the interface itself; we now proceed to analyse the magnitude of the flow velocity at the interface. The flow velocity is decomposed into two components, a normal component 
$u_n\equiv \boldsymbol {u}\boldsymbol {{\cdot }} \boldsymbol {n}$ aligned with the outward-pointing normal to the interface 
$\boldsymbol {n}$ and a tangential component 
$u_t\equiv \|\boldsymbol {u}-u_n \boldsymbol {n}\|$. The sign of the normal component is important; the interface is advected with the flow so positive 
$u_n$ occurs in places where the interface moves outward in the direction of the carrier phase, and negative 
$u_n$ occurs in places where the interface moves inward in the direction of the dispersed phase. Due to volume conservation of the phases, 
$u_n$ has zero mean. On the other hand, the choice of the tangential direction in the plane is arbitrary: it is taken as the remainder of the subtraction of the normal component from the total velocity. For this reason, we have only positive values of the tangential velocity 
$u_t$.
In figure 11(a) the PDFs of the normal component of the velocity 
$u_n$ show peaks at relatively low positive 
$u_n$ and are negatively skewed in all cases, i.e. extreme negative values (inward fluid velocity) are more probable than extreme positive values (outward fluid velocity). To elucidate the cause of the skewed distributions, we also show histograms of the logarithm of the dissipation 
$\epsilon$ for case C40. The tails of these histograms can be attributed to extreme events and, hence, to intermittency in the flow (Kaneda & Morishita Reference Kaneda and Morishita2012). We see that the flow outside the droplets gives slightly wider tails and, hence, has higher intermittency than that inside the droplets (also seen by Crialesi-Esposito et al. (Reference Crialesi-Esposito, Boffetta, Brandt, Chibbaro and Musacchio2023) for fluid velocity differences inside and outside droplets). Therefore, we attribute the increased likelihood of extreme inward velocities at the interface to increased intermittency of the flow in the bulk phase. Returning to the main panel of figure 11(a), we note that cases with a higher 
$We_e$ show a wider distribution of 
$u_n/u'$. A higher effective Weber number implies a more deformable interface with a lesser damping effect on the normal velocity: extreme (positive and negative) events become more probable as surface tension is reduced. The difference among surfactant-laden and clean cases at similar effective Weber numbers is reported with two-colour lines (rescaled by a factor of two for improved readability). We observe that the surfactant suppresses extreme events, especially at low values of the Weber number (high surface tension), when Marangoni stresses are greatest in magnitude. A similar turbulence suppression effect was seen by Shen, Yue & Triantafyllou (Reference Shen, Yue and Triantafyllou2004) for surfactants in free shear flows, and was attributed to the elasticity of the surfactant-laden interface.
(a) Normal component of the flow velocity at the interface 
$u_n$ and (b) tangential component of the flow velocity at the interface 
$u_t$. Two-colour lines show the difference between the PDF of surfactant-laden and clean pairs of cases with similar 
$We_e$ (these data are magnified by a factor two for better reading). The inset in panel (a) shows histograms of the logarithm of the dissipation 
$\epsilon$ inside (dashed line) and outside (solid line) the droplets for case C40, where 
$m_\epsilon$ and 
$s_\epsilon$ are the mean and standard deviation of 
$\ln \epsilon$, respectively. The inset in panel (b) shows the decomposition of the flow velocity along a direction normal (
$\boldsymbol {n}$) and tangential (
$\boldsymbol {t}$) to the interface.
When considering the tangential component, we observe a reversal of the trend: surfactant, via Marangoni stresses, increases the probability of flow tangential to the interface. Two-colour lines show the difference between surfactant-laden and clean cases at similar 
$We_e$, with the data rescaled by a factor of two for ease of reading. It is clear that the presence of surfactant increases the probability of finding tangential velocities in the range 
$u' \lesssim u_t \lesssim 2u'$. All surfactant-laden cases have approximatively a similar value of the peak of the distribution, which is slightly larger than that of the surfactant-free cases, indicating that surfactant-laden cases have a higher probability of large values of the tangential velocity. Interestingly, as the Weber number increases, the peak shifts to slightly higher values, and the likelihood of large tangential velocity increases as well, as shown by the two-colour lines. This result suggests that while increasing the flow velocity tangential to the interface and suppressing large normal components, Marangoni stresses also have a modulating effect on the tangential component.
4. Conclusions
In this study we perform direct numerical simulations of surfactant-laden droplets in homogeneous isotropic turbulence. The interfacial dynamics is solved using an MTHINC volume-of-fluid method coupled with a phase-field-based approach to simulate surfactant dynamics. By examining droplet morphology and local flow statistics, we can shed light on the interfacial characteristics and dynamics of these complex systems.
The Kolmogorov–Hinze length scale is a fundamental quantity in multiphase flows laden with clean droplets or bubbles. Our numerical results confirm that the Kolmogorov–Hinze framework can be extended to surfactant-laden droplets by using an averaged surface tension value, thus accounting for the presence of surfactant. In the configuration we adopted in this study, the role of Marangoni stresses is minor, thus, the effect of surfactant can be approximated as a simple average surface tension reduction. We indeed observe a collapse of most statistics when plotted as a function of 
$d/d_H$. We compute the droplet size distribution for all clean and surfactant-laden cases and verify that (i) the Kolmogorov–Hinze scale effectively separates the breakage- and coalescence-dominated regimes, and (ii) the power-law scaling for these two regimes can be applied to surfactant-laden droplets. The combined results on the deformation of the droplets, i.e. aspect ratio and sphericity, prove that droplets smaller than the Kolmogorov–Hinze scale have a relatively compact and regular shape (spheroid- or ellipsoid-like shapes), whereas droplets larger than the Kolmogorov–Hinze scale have a coiled, filamentous shape, supporting previous observations of filamentous water drops (Villermaux & Bossa Reference Villermaux and Bossa2009; Jackiw & Ashgriz Reference Jackiw and Ashgriz2021). The filamentous droplets are found to have an average diameter that is independent of the overall droplet size 
$d$ and is equal to the Kolmogorov–Hinze scale, further evidencing the relevance of this length scale.
The very different shapes of large and small droplets have direct implications on the total area of the interface, which is a crucial parameter in determining the overall exchange of species, momentum and energy among the carrier and the dispersed phase. We report the existence of two regimes, separated by the Kolmogorov–Hinze scale: the area of droplets smaller than the Kolmogorov–Hinze scale is proportional to the square of the characteristic size of the droplet, whereas it is proportional to the cube of the characteristic size for droplets larger than the Kolmogorov–Hinze scale. The two different scalings can be directly traced back to the shape of the droplets, spheroid-like below the Kolmogorov–Hinze scale and filamentous above the Kolmogorov–Hinze scale. The large and filamentous droplets are coiled up, as indicated by their aspect ratio. Thus, we investigate their self-connectedness, using the Euler characteristic to count the number of handles and voids on each droplet. To the best of our knowledge, this is the first time the number of handles and voids on a droplet has been measured. We find that the number of handles depends directly on the size of the droplet and its surface tension, as data from all cases collapse on a single curve when normalized by the Kolmogorov–Hinze scale; we also provide a scaling for the linear density of handles, 
$h/l=0.06/d_H$. Conversely, the number of voids depends on the droplet size alone. Our interpretation is that the restoring action of surface tension reduces the lifetime of a handle, whereas the dynamics of a void are unaffected by surface tension. Going further, the Poincaré-Hopf theorem relates the Euler characteristic 
$\chi$ of an interface to the number of topological defects in any tangent vector field (e.g. fluid velocity, alignment of molecules, stresses) at the interface (Maroudas-Sacks et al. Reference Maroudas-Sacks, Garion, Shani-Zerbib, Livshits, Braun and Keren2021) – a well-known example of this for 
$\chi =2$ is the hairy ball theorem. A future investigation may count topological defects on the surface of droplets and relate the number to their ability to resist breakup.
Results from the morphology of the droplets highlight the validity of the Kolmogorov–Hinze scale for both clean and surfactant-laden flows: a rescaled value of the surface tension, accounting for the average surface tension reduction induced by the surfactant, can be effectively used to define the Kolmogorov–Hinze scale. This finding suggests that the effect of surfactant on droplets in homogeneous isotropic turbulence can be mainly summarized as a reduction in surface tension. The lack of a large-scale and time-persisting velocity difference among the carrier and dispersed phase, as found in up-flow and down-flow configurations (Takagi et al. Reference Takagi, Ogasawara and Matsumoto2008; Lu et al. Reference Lu, Muradoglu and Tryggvason2017), prevents the formation of significant, large-scale Marangoni stresses at the interface. Hence, in our simulation set-up, Marangoni stresses play a minor, local role, with negligible effects on the statistics of the droplets. Local flow statistics better show the effect of these tangential stresses: Marangoni stresses modulate the flow at the interface by reducing the velocity component normal to the interface and increasing the tangential component. When computing the flow topology parameter at the interface, we find that Marangoni stresses increase the elongational component and reduce rotational flow at the interface. This result is coherent with the action of Marangoni stresses generating an elongational type of flow with sources corresponding to low-surface-tension regions and sinks to high-surface-tension regions. Inside the droplets, Marangoni stresses reduce elongational flow and increase the pure shear contribution.
In conclusion, we find that a statistically homogeneous and isotropic flow allows for a simplified treatment of the surfactant effects. Results from Hinze's (Reference Hinze1955) theory can thus be applied to surfactant-laden flows, by considering the average surface tension reduction operated by the surfactant. For both clean and surfactant-laden flows, the Kolmogorov–Hinze scale separates two regimes characterized by very different scaling for the surface area of the droplets: a significant increase in the surface area is observed for droplets larger than the Kolmogorov–Hinze scale. This latter result has important implications for environmental and industrial multiphase flows, where the interface serves as a conduit for all species, momentum and energy transfers among the phases.
Acknowledgements
All authors acknowledge the computational resources provided on the Deigo cluster by the Scientific Computing and Data Analysis section of Research Support Division at OIST, and on Oakbridge-CX provided by the Information Technology Center (The University of Tokyo) through the HPCI System Research Project (project ID: hp220100). I.C. would like to thank Julian Katzke (OIST) for guidance with the rendering software Blender, and David O'Connell (OIST) for helpful discussion of the Euler characteristic.
Funding
The research was supported by the Okinawa Institute of Science and Technology Graduate University (OIST) with subsidy funding from the Cabinet Office, Government of Japan.
Declaration of interests
The authors report no conflict of interest.
Data availability statement
The plotted data are available at https://groups.oist.jp/cffu/Cannon2024JFluidMech, and the code used for calculating droplet morphology can be found at https://github.com/marco-rosti/CFF-dropStats.git.
Appendix A. Euler characteristic of a droplet interface
Here we use concepts from algebraic topology to derive (3.5), which relates the Euler characteristic of a droplet interface to the number of voids and handles it contains. The Euler characteristic and genus are well defined for a single surface; however, as seen in the image in figure 8(b), our droplets can have many distinct surfaces. For a general droplet with a number of voids 
$v$, it has an outer surface with Euler characteristic 
$\chi _0$, and the voids create 
$v$ distinct inner surfaces with Euler characteristics 
${\chi _1, \chi _2,\ldots \chi _v}$. We define the Euler characteristic of the droplet's interface as the sum of the Euler characteristics of all the surfaces,
\begin{equation} \chi=\sum_{i=0}^v\chi_i.\end{equation}
Both the inner and outer surfaces of the droplet can have handles. For any distinct orientable surface, the number of handles 
$h_i$ is its genus, and the genus is related to the Euler characteristic by 
$\chi _i=2-2h_i$ (Massey Reference Massey1997, p. 30). Hence, (A1) can be written and rearranged to
\begin{equation} \chi=\sum_{i=0}^v (2-2h_i)=2(1+v)-\sum_{i=0}^v 2h_i . \end{equation}
The number of handles on the drop is the sum of the handles on all surfaces of the droplet, i.e. 
$h=\sum _{i=0}^v h_i$, so we can write
\begin{equation} \chi=2+2v-2h. \end{equation}This can be easily rearranged to obtain (3.5).
Appendix B. Validation of the surfactant model
We report in this appendix the validation tests of the proposed numerical method. Two benchmark tests are presented and compared against existing theories, experiments and simulations.
We verify at first the adsorption dynamics of the surfactant at the interface: we start from a uniform surfactant concentration, equal to the surfactant concentration in the bulk 
$\psi _b$, and let the surfactant adsorb onto the interface (located at 
$x=L/2$ in figure 12a). Initially, the system is out of equilibrium and surfactant diffuses towards the interface to restore the equilibrium; at equilibrium, the chemical potential is equal everywhere. We can thus equate the chemical potential in the bulk (
$\hat \phi =\pm 1$, surfactant concentration 
$\psi _b$) and at the interface (
$\hat \phi =0$, surfactant concentration 
$\psi _0$) to obtain the so-called Langmuir isotherms (Engblom et al. Reference Engblom, Do-Quang, Amberg and Tornberg2013). For given values of the parameters of the chemical potential, 
$\alpha$, 
$\beta$ and 
$\gamma$, the Langmuir isotherm relates the surfactant concentration at the interface 
$\psi _0$ to that in the bulk 
$\psi _b$, in equilibrium conditions
\begin{equation} \psi_0=\frac{\psi_b}{\psi_b+(1-\psi_b) \, {\rm e}^{-({\beta+\gamma})/{2\alpha}}}.\end{equation}
Figure 12(a) shows the set-up we use to test our code against the Langmuir isotherm benchmark: a flat interface is located at 
$x=L/2$, and surfactant is initially uniformly distributed with a concentration equal to 
$\psi _b$. The flow is initially at rest and, due to the absence of any forcing, stays at rest throughout the entire simulation; similarly, the interface does not move from its initial position. At the boundaries 
$x=0$ and 
$x=L$, we impose a far-field value of surfactant, 
$\psi =\psi _b$, thus allowing surfactant to enter the system. In the 
$y$ and 
$z$ directions, we impose periodic boundary conditions. This benchmark is, in principle, a one-dimensional test (
$x$ direction). However, we perform three-dimensional numerical simulations in order to use the very same solver described in § 2; all variables are uniform in the 
$y$ and 
$z$ directions. We resolve the flow and volume-of-fluid equations on a 
$N_x \times N_y \times N_z=100\times 4\times 4$ computational grid, and use a more refined grid for the surfactant transport equations, 
$500\times 4\times 4$. The smoothing width of the interface is 
$3\varDelta$. We use 
$\beta /\alpha =0.741$ and explore a range of values for the energy cost in the bulk 
$\gamma /\alpha =\{0.0741, 0.741, 3.70, 7.41, 11.1, 14.8\}$. Figure 12(b) shows the simulated interfacial surfactant concentrations 
$\psi _0$ once they had reached equilibrium for various values of 
$\psi _b$. The values obtained from our numerical simulations fall on top of the corresponding Langmuir isotherms, proving that the implemented numerical method can correctly capture the surfactant dynamics.
(a) Adsorption of surfactant onto an interface at 
$x=L/2$. The right-hand axis shows the smoothed colour function 
$\hat \phi$ (in red), which has a smoothing width 
$3\varDelta$. The left-hand axis shows the surfactant concentration (in blue) at four time instants during the simulation. The bulk surfactant concentration is 
$\psi _b=0.05$ and the energy cost of surfactant in the bulk is 
$\gamma =11.1\alpha$. (b) The equilibrium interfacial surfactant concentration 
$\psi _0$ for a range of 
$\psi _b$ and 
$\gamma$. Lines show Langmuir isotherms, given by (B1), and markers show the results of our simulations. The value of 
$\gamma$ is represented by shading from light blue to dark blue.
The second benchmark aims to verify the surfactant transport on a moving interface and the computation of surface tension forces by measuring the deformation of a droplet in shear flow. Taylor (Reference Taylor1934) quantified the deformation of a droplet in a shear flow using the deformation parameter
\begin{equation} D\equiv\frac{L-B}{L+B}, \end{equation}
where 
$L$ and 
$B$ are the droplet's largest and smallest principal diameters, respectively. Figure 13(a) shows the simulation set-up that we use to reproduce Taylor's experiment. The effective Capillary number 
$Ca_e\equiv a U \mu /h\langle \sigma \rangle$ describes the ratio of viscous forces to surface tension forces in the system, where 
$h$ is the domain half-height, 
$a=0.4h$ is the initial radius of the droplet, 
$\pm U$ is the fluid velocity of the top and bottom walls respectively, 
$\mu$ is the dynamic viscosity, which is the same for the bulk and droplet phases, and 
$\langle \sigma \rangle$ is the (average) surface tension at the droplet interface. As was previously demonstrated by Soligo et al. (Reference Soligo, Roccon and Soldati2020a), there are negligible differences between the two-dimensional and three-dimensional cases in the limit of small Reynolds and Capillary numbers, as those considered here. Hence, we chose to perform two-dimensional numerical simulations to reduce the computational cost of the benchmark simulations. We select a computational grid with 
$N_x \times N_y \times N_z=600\times 200 \times 4$ points for the flow and volume of fluid variables, and we use a twice more refined grid to discretize the surfactant transport equation. As we use the same three-dimensional solver introduced in § 2, we have to use four grid points in the 
$z$ direction; however, all variables are uniform in the 
$z$ direction.
Deformation of droplets in shear flows. (a) The simulated 2-D domain (partially shown here) with velocity boundary conditions 
$u=U$ at 
$y=2h$, 
$u=-U$ at 
$y=0$, and periodic boundary conditions at 
$x=0$ and 
$x=6h$. We show the location of the droplet interface for a clean and surfactant-laden case. (b) Dependence of the steady-state deformation parameter 
$D$ on the effective Capillary number 
$Ca_e$. Our simulation results are shown with solid markers and values from the literature are shown with empty markers. Clean droplets are marked in red and surfactant-laden droplets are marked in blue. The solid line is the analytical relation from (Taylor Reference Taylor1934) with the confinement correction proposed by Shapira & Haber (Reference Shapira and Haber1990).
We run a number of simulations with various values of 
$Ca_e$ and three values of the initial surfactant concentration in the bulk phase 
$\psi _b\in \{0,0.01,0.02\}$, where 
$\psi _b=0$ is the clean droplet case. In all cases, the Reynolds number of the flow 
$Re\equiv \rho U h/\mu$ is 
$Re=0.1$. The droplets are initially circular and deform in the shearing flow. Figure 13(b) shows our measured values of the deformation parameter 
$D$ once they had reached a steady state. We compare our results (solid markers) with two-dimensional droplets simulated using the boundary integral method (Rallison Reference Rallison1981), experiments of three-dimensional droplets (Guido & Simeone Reference Guido and Simeone1998), three-dimensional droplets simulated using the volume-of-fluid method (Li, Renardy & Renardy Reference Li, Renardy and Renardy2000), three-dimensional droplets simulated using the boundary integral method with insoluble surfactants (Bazhlekov et al. Reference Bazhlekov, Anderson and Meijer2006), three-dimensional droplets simulated using dissipative particle dynamics (Pan, Phan-Thien & Khoo Reference Pan, Phan-Thien and Khoo2014) and surfactant-laden droplets simulated using a phase-field method (Soligo et al. Reference Soligo, Roccon and Soldati2020a). Our results are in good agreement with all of the above and also closely follow the analytical result from Taylor (Reference Taylor1934), with the confinement correction proposed by Shapira & Haber (Reference Shapira and Haber1990), i.e. for droplets and carrier fluid having the same viscosity, the Taylor deformation parameter is equal to
\begin{equation} D=\frac{35}{32} Ca_e\left[1+C_{S H} \frac{3.5}{2}\left(\frac{a}{2h}\right)^3\right], \end{equation}
where 
$C_{SH}=5.6996$ is a numerical coefficient accounting for the confinement due to the top and bottom boundaries (Shapira & Haber Reference Shapira and Haber1990).
In this appendix we have shown that the proposed numerical method can accurately simulate the transport of surfactant over moving and deforming interfaces, as well as the action of surfactant on surface tension.
Open access
