F Thiesset

B Duret

T Ménard

C Dumouchel

J Reveillon

F X Demoulin

Liquid transport in scale space

When a liquid stream is injected into a gaseous atmosphere, it destabilizes and continuously passes through different states characterized by different morphologies. Throughout this process, the flow dynamics may be different depending on the region of the flow and the scales of the involved liquid structures. Exploring this multi-scale, multi-dimensional phenomenon requires some new theoretical tools, some of which need yet to be elaborated. Here, a new analytical framework is proposed on the basis of two-point statistical equations of the liquid volume fraction. This tool, which originates from single phase turbulence, allows notably to decompose the fluxes of liquid in flow-position space and scale space. Direct Numerical Simulations of liquid-gas turbulence decaying in a triply periodic domain are then used to characterize the time and scale evolution of the liquid volume fraction. It is emphasized that two-point statistics of the liquid volume fraction depend explicitly on the geometrical properties of the liquid-gas interface and in particular its surface density. The stretch rate of the liquid-gas interface is further shown to be the equivalent for the liquid volume fraction (a non diffusive scalar) of the scalar dissipation rate. Finally, a decomposition of the transport of liquid in scale space highlights that non-local interactions between non adjacent scales play a significant role.

Introduction

Flows of air and water like fluids are plethoric in nature as in our everyday life. They were consubstantial to the emergence of life on earth and remain essential to its sustainability. Some involve both a liquid and gas phase, which are referred to as twophase flows. Our everyday experience reveals that e.g. a water jet ensuing into quiescent air has strong propensity to break up in liquid fragments, subsequently atomizing into a stream of dispersed droplets in a gaseous atmosphere, called spray. Usage of sprays dates back to the Paleolithic for painting hand stencils and are now employed in a very wide range of situations: irrigation, cosmetics, pharmaceutics, manufacturing, decontamination. Therefore, sprays have strong implications in the economy, health, safety, and well-being of modern societies.

According to [START_REF] Villermaux | Fragmentation[END_REF], there exist three classes of theoretical elaborations for describing the atomization of liquid streams: the sequential cascade models, the aggregation models, and the maximum entropy formalism. Their key benefit is to provide a plausible vision of the mechanisms of rupture and spray formation, leading to some very useful mathematical formulations for the spray drop-size distribution in the far field. The first class of elaborations dates back to Kolmogorov (1941b) and considers the atomization process as a successive reduction of 'mother drops' into a hierarchy of smaller 'daughter drops'. This approach in its original and modern versions [START_REF] Gorokhovski | Analyses of Kolmogorov's model of breakup and its application into Lagrangian computation of liquid sprays under air-blast atomization[END_REF] thus conceptualizes the atomization process as a cascade mechanism, analogous to single fluid turbulence (Kolmogorov 1941a), where the sense of evolution is directed towards small scales. The second class of models is due mainly to the work by Villermaux and co-workers [START_REF] Villermaux | Ligament-mediated spray formation[END_REF]Marmottant & Villermaux 2004). They postulated that droplets emanate from the rupture of detached ligaments which are formed during the primary atomization process. These ligaments are assumed to be virtually constituted of spherical blobs, which successively aggregate yielding a given droplet size distribution in the spray. Hence, the aggregation scenario indicates that the transfer among the different scales is directed towards the large scales. Finally, the last class of theory relies on the Maximum Entropy Formalism [START_REF] Sellens | A simplified prediction of droplet velocity distributions in a spray[END_REF][START_REF] Dumouchel | The maximum entropy formalism and the prediction of liquid spray dropsize distribution[END_REF]). The idea is to consider, again, the liquid as being constituted as a set of elementary bricks and to compute the most probable distribution of disjoint clusters given some constraints that are imposed on the system. In the Maximum Entropy Formalism, the process of cluster formation is done in one step and there is no explicit description of the underlying dynamics or kinematics.

It is worth mentioning that these three classes of models do not apply to the same flow regime or flow region. Indeed, the cascade model is intended to address high Weber numbers flow where the effect of surface tension is small compared to inertial forces. On the contrary, the aggregation scenario applies to the final stage of atomization when break-up of liquid structures is caused by surface tension effects. Consequently, depending on the local Weber number, there could exist some zones in a flow complying more with the cascade scenario and some other regions where the aggregation process is at play. In general situations, it remains very unclear in which particular direction, if one were to exist, the evolution in scale space preferentially occurs [START_REF] Herrmann | An outline of a LES subgrid model for liquid/gas phase interface dynamics[END_REF].

All aforementioned models emphasize the cornerstone notion of constitutive scale (mother and daughter drops, blobs) and scale transfer. However, in both visions, this notion is rather ambiguous and lacks of a handy, measurable definition. If a typical scale were to be unambiguously defined, it could only correspond to the diameter of a sphere (ligament) while the liquid structures are generally far from being spherical (cylindrical) in the primary or secondary atomization zone. A direct consequence is that predictions of such theories can be tested mostly in the very far field of the spray, where droplets are spherical for the 'drop-size' to be defined and measured. Hence, they cannot be validated in the region of the flow where the processes they intend to describe are acting.

Therefore, although the cascade and the aggregation scenario have revealed to be extremely beneficial to provide the mathematical foundations for the spray drop size distribution, tackling the complexity and the multi-scale facets of liquid atomization requires elaborating a new observable:

(i) which applies all along the atomization process from the liquid injection to the spray dispersion, and irrespectively to the flow regime, (ii) which unambiguously define the notion of scale thereby allowing the different mechanisms at play and the effect of some physical parameters (surface tension, density, viscosity, etc) to be scrutinized for each scale, (iii) whose evolution can be explored in the light of a transport equation derived from first principles.

Some observables have been proposed to fulfill these requirements. First, the volumebased scale distribution [START_REF] Dumouchel | Multi-scale analysis of simulated capillary instability[END_REF], which is the 3D extension of the surface based scale distribution [START_REF] Dumouchel | Analysis of two-dimensional liquid spray images: the surface-based scale distribution[END_REF]) has been introduced and has shown great potential for characterizing the multi-scale facets of atomization. Another approach was proposed by [START_REF] Canu | Where does the droplet size distribution come from?[END_REF] who defined the scale distribution on the basis of intrinsic geometrical properties of the liquid-gas interface, namely its mean and Gaussian curvatures. Recent theoretical work (Thiesset et al. 2019b) have shown that the volume based scale distribution can be expressed in terms of the geometrical metric of the liquid gas interface (its surface area, mean and Gaussian curvatures). There thus exists a strong connection between the two latter approaches. The present work follows these lines and aims at proposing a new observable.

Besides two-phase flows, the concept of hierarchized scales among which energy is transferred has been extensively explored in the context of single phase turbulent flows. This idea dates back to [START_REF] Richardson | The supply of energy from and to atmospheric eddies[END_REF] and has deeply inspired von [START_REF] Von Kármán | On the statistical theory of isotropic turbulence[END_REF] and Kolmogorov (1941a) who were first to derive the transport equations for the two-point statistics, (i.e. correlation or structure functions) in the specific context of single phase homogeneous isotropic turbulence. A similar framework has been followed by [START_REF] Yaglom | On the local structure of a temperature field in a turbulent flow[END_REF] who proposed the two-point statistical equation of a passive scalar field (e.g. temperature) at asymptotically high Reynolds number. Considering an analysis of the relevant metrics (e.g. velocity, temperature, concentration) not only at one point but at two points in space allows to mathematically define the notion of scale as the straight distance between the two points considered, thereby allowing to assess the energy budget for each scale of the turbulence spectrum. Two-point statistical equations are thus sometimes referred to as scale-by-scale budgets. A generalisation of such equations to non-isotropic yet homogeneous flows was presented by Monin & Yaglom (1975) (see p. 402). In recent years, general equations for inhomogeneous anisotropic cases have been developed, see e.g. [START_REF] Hill | Exact second-order structure-function relationships[END_REF]; [START_REF] Danaila | Progress in studying small-scale turbulence using exact two-point equations[END_REF][START_REF] Danaila | Yaglom-like equation in axisymmetric anisotropic context[END_REF]. These equations are now referred to as generalised Karman-Howarth-Kolmogorov-Monin (KHKM) equations. As far as the scalar field is concerned, the generalised Yaglom's equation was derived by e.g. [START_REF] Gauding | Generalised scale-by-scale energy-budget equations and large-eddy simulations of anisotropic scalar turbulence at various schmidt numbers[END_REF]; [START_REF] Togni | Physical and scale-by-scale analysis of Rayleigh-Bénard convection[END_REF].

Scale-by-scale budgets have been used extensively in recent years to characterize the statistical behaviour of some canonical flows such as grid or box turbulence [START_REF] Danaila | A generalization of Yaglom's equation which accounts for the large-scale forcing in heated decaying turbulence[END_REF][START_REF] Orlandi | Dependence of the non-stationary form of Yaglom's equation on the Schmidt number[END_REF][START_REF] Antonia | Approach to the 4/5 law in homogeneous isotropic turbulence[END_REF][START_REF] Valente | The energy cascade in grid-generated non-equilibrium decaying turbulence[END_REF][START_REF] Gomes-Fernandes | The energy cascade in near-field non-homogeneous non-isotropic turbulence[END_REF] in presence of a mean scalar gradient [START_REF] Danaila | Effect of gradient production on scalar fluctuations in decaying grid turbulence[END_REF][START_REF] Watanabe | Statistics of a passive scalar in homogeneous turbulence[END_REF], forced isotropic turbulence [START_REF] Moisy | Kolmogorov equation in a fully developed turbulence experiment[END_REF][START_REF] Gotoh | Velocity field statistics in homogeneous steady turbulence obtained using a high-resolution direct numerical simulation[END_REF], homogeneous shear turbulence [START_REF] Casciola | Scale-by-scale budget and similarity laws for shear turbulence[END_REF], fully developed channel flows (Danaila et al. 2001;[START_REF] Marati | Energy cascade and spatial fluxes in wall turbulence[END_REF][START_REF] Saikrishnan | Reynolds number effects on scale energy balance in wall turbulence[END_REF][START_REF] Cimarelli | Paths of energy in turbulent channel flows[END_REF], wall turbulence [START_REF] Cimarelli | Sources and fluxes of scale energy in the overlap layer of wall turbulence[END_REF], rotating turbulence [START_REF] Lamriben | Direct measurements of anisotropic energy transfers in a rotating turbulence experiment[END_REF][START_REF] Campagne | Direct and inverse energy cascades in a forced rotating turbulence experiment[END_REF], turbulent wake flows [START_REF] Thiesset | Scale-by-scale turbulent energy budget in the intermediate wake of two-dimensional generators[END_REF][START_REF] Portela | The turbulence cascade in the near wake of a square prism[END_REF], round jet flows [START_REF] Danaila | Progress in studying small-scale turbulence using exact two-point equations[END_REF]Burattini et al. 2005a;Thiesset et al. 2014a), Rayleigh-Bénard convection [START_REF] Togni | Physical and scale-by-scale analysis of Rayleigh-Bénard convection[END_REF] and detached turbulent flows behind a bump [START_REF] Mollicone | Turbulence dynamics in separated flows: the generalised Kolmogorov equation for inhomogeneous anisotropic conditions[END_REF]. Scale-by-scale budgets have also been generalised to account for density variations [START_REF] Lai | A Kármán-Howarth-Monin equation for variable-density turbulence[END_REF]), viscosity variations (Danaila et al. 2017;[START_REF] Gauding | One-point and two-point statistics of homogeneous isotropic decaying turbulence with variable viscosity[END_REF], and the modulation of small-scales statistics due to the presence of coherent structures (Thiesset et al. 2014b).

The main benefit of scale-by-scale budgets relies first on their capacity to decompose the transport of turbulent quantities in geometrical space (i.e. from one position in the flow to another) and in scale space (e.g. from large to small scales). Therefore, the effect of spatial inhomogeneity can easily be distinguished from what is commonly referred to as a cascade process or more generally a transfer in scale space. Second, it is worth stressing that such generalized two-point equations derive from first principles. They are thus exact and do rely on particular assumptions about the flow. Third, since in Large-Eddy Simulations (LES), the notion of scale is inherent to the filtering operation, scale energy budgets have further shown to be useful to better characterizing the energy transfer between the resolved and subgrid scales [START_REF] Gualtieri | Preservation of statistical properties in large-eddy simulation of shear turbulence[END_REF][START_REF] Cimarelli | Analysis of the Kolmogorov equation for filtered wallturbulent flows[END_REF][START_REF] Cimarelli | Sources and fluxes of scale energy in the overlap layer of wall turbulence[END_REF]. LES models can also be developed on the basis of such a framework [START_REF] Lesieur | New trends in large-eddy simulations of turbulence[END_REF][START_REF] Lévêque | Shear-improved Smagorinsky model for large-eddy simulation of wall-bounded turbulent flows[END_REF][START_REF] Cui | A new subgrid eddy-viscosity model for large-eddy simulation of anisotropic turbulence[END_REF]. Finally, concomitantly with the development of two-point statistical equations, has emerged a very powerful tool based on the notion of self-similarity or self-preservation [START_REF] Von Kármán | On the statistical theory of isotropic turbulence[END_REF]Kolmogorov 1941a;George 1992a;Antonia et al. 2003, for instance). Self-preservation states that two-point statistics can be invariant throughout their time and/or space evolutions, when normalized by an appropriate set of variables, called similarity variables. This hypothesis has led to significant contributions to the study of homogeneous isotropic turbulence [START_REF] Von Kármán | On the statistical theory of isotropic turbulence[END_REF][START_REF] Stewart | Similarity and self-preservation in isotropic turbulence[END_REF]George 1992a;[START_REF] Speziale | The energy decay in self-preserving isotropic turbulence revisited[END_REF]Antonia et al. 2003, among others) and to scalar fluctuations evolving in isotropic turbulence (George 1992b;[START_REF] Chasnov | Similarity states of passive scalar transport in isotropic turbulence[END_REF][START_REF] Gonzalez | The approach to self-preservation of scalar fluctuations decay in isotropic turbulence[END_REF][START_REF] Antonia | Similarity of decaying isotropic turbulence with a passive scalar[END_REF]Antonia et al. 2004)). Self-similarity has also been applied to homogeneous shear turbulence [START_REF] George | The self-preservation of homogeneous shear flow turbulence[END_REF][START_REF] Danaila | Progress in studying small-scale turbulence using exact two-point equations[END_REF], wake flows [START_REF] Bevilaqua | Turbulence memory in self-preserving wakes[END_REF][START_REF] George | The self-preservation of turbulent flows and its relation to initial conditions and coherent structures[END_REF][START_REF] Moser | Self-similarity of time-evolving plane wakes[END_REF][START_REF] Antonia | Approach towards self-preservation of turbulent cylinder and screen wakes[END_REF]Ewing et al. 2007b;[START_REF] Thiesset | Scale-by-scale turbulent energy budget in the intermediate wake of two-dimensional generators[END_REF], turbulent jets [START_REF] Bradbury | The structure of a self-preserving turbulent plane jet[END_REF][START_REF] Antonia | Measurements of dissipation rate and some other characteristics of turbulent plane and circular jets[END_REF][START_REF] George | The self-preservation of turbulent flows and its relation to initial conditions and coherent structures[END_REF]Ewing et al. 2007a;Burattini et al. 2005b;Thiesset et al. 2014a) and turbulent shear-layers [START_REF] Bradshaw | The effect of initial conditions on the development of a free shear layer[END_REF][START_REF] Brown | On density effects and large structure in turbulent mixing layers[END_REF][START_REF] Fiedler | On turbulence structure and mixing mechanism in free turbulent shear flows[END_REF]Hussain & Zedan 1978b,a). It is important to stress that the quest for self-similar solutions has motivated many aspects of research in turbulence since it has the potential of reducing partial differential equations to ordinary differential equations leading to simple closures for the scale-by-scale budgets. In addition, when the most appropriate similarity variables are identified, it deeply helps providing insights into the predominant mechanisms at play in the evolution of the two-point statistics.

Note that another way of assessing the scale content of a given phenomena is by inferring the spectral Fourier components of the fluctuating field. However, this is likely to raise some issues in some cases. First, spectra can be computed only over homogeneous directions with periodic boundary conditions. This is a very stringent restriction that makes spectra inapplicable to all flow situations. Correlation or structure functions do not suffer of such a limitation. Second, as discussed by notably [START_REF] Lucci | Modulation of isotropic turbulence by particles of taylor length-scale size[END_REF] and also [START_REF] Duret | DNS analysis of turbulent mixing in two-phase flows[END_REF], two-phase flows reveal local discontinuities associated with the presence of the interface. This can result in oscillations in the large wavenumbers range which corrupt the physical interpretation of energy spectra. This issue can be alleviated in [START_REF] Lucci | Modulation of isotropic turbulence by particles of taylor length-scale size[END_REF] by considering not the Eulerian wavenumber spectrum but the Lagrangian frequency spectrum of fluid trajectories at the sacrifice that it does not emphasize the scale but the time dynamics of the flow under consideration.

In the present study, the machinery of two-point statistical analysis and self-similarity is adapted to a relevant metric of two-phase flows with the goal of exploring the scale distribution and scale energy transfer of the liquid phase. The scalar of interest is here the liquid volume fraction [START_REF] Hirt | Volume of fluid (VOF) method for the dynamics of free boundaries[END_REF][START_REF] Scardovelli | Direct numerical simulation of free-surface and interfacial flow[END_REF]. In the long run, our objectives in deriving such a new framework is to provide insights into (i) the preferential direction (if any) and amplitude of the transfer among the different scales and flow positions during liquid-atomization, and (ii) the effect of different physical parameters (surface tension, viscosity and density, inflow velocity conditions), and the range of scales over which these parameters have an influence. Only point (i) is addressed for the present work. Although this framework applies irrespectively of the investigated zone within the flow and independently of the injection regime, this very first study focuses only on a turbulent two-phase flow decaying in a triply periodic box.

This paper is organized as follows. First, the scale-by-scale budget for the liquid volume Dark regions correspond to the liquid-phase and white regions to the gas-phase The interface is also displayed as black curves. The two-points x `and x ´used to calculate increments are also represented together with the midpoint X and the separation vector r.

fraction is presented. The derivation is given in great details in the Appendix, so that the reader can easily follow the different calculation steps. In this first section, we put particular emphasis on the physical interpretation of the different terms present in the equations. Then, the numerical database allowing this new framework to be appraised is presented. Specific care is given to the appropriateness of the spatial resolution and statistical convergence. Then, the two-point statistical analysis or the liquid volume fraction flow is provided. Conclusions are drawn in a last section.

Theoretical considerations

The liquid volume fraction function

The liquid volume fraction (abbreviated LVF and noted φ) designates the portion of liquid contained within a given volume. The LVF φ " 1 in the liquid phase and φ " 0 in the gas phase. An example of such a field function is portrayed in Fig. 1.

When the flow is incompressible with no phase change, the transport equation for the LVF reads

B t φ `∇x ¨uφ " 0 (2.1)
where B t " B{Bt and ∇ x is the gradient operator with respect to the coordinate x. The equation for the LVF is therefore that of a non-diffusive passive scalar. The effect of fluid viscosity and density jumps together with surface tension of the liquid-gas interface is implicit in this equation through the role they have on the velocity and pressure fields.

LVF transport in 7 dimensional space

The algebra for deriving the two-point equation of the liquid volume fraction is similar to what was presented in e.g. [START_REF] Hill | Exact second-order structure-function relationships[END_REF]; [START_REF] Danaila | Progress in studying small-scale turbulence using exact two-point equations[END_REF] for the velocity field. Application to such a machinery to the LVF is detailed in Appendix A to C.

A convenient form of the scale-by-scale budget is derived in Appendix C and is given by B t pδφ 1 q 2 `∇r ¨pδuq pδφ 1 q 2 `2pδφ 1 qpδu 1 q ¨∇r pδφq `∇X ¨pσuq pδφ 1 q 2 `2pδφ 1 qpσu 1 q ¨∇X pδφq " 0 (2.2) which can formally be written as [START_REF] Mollicone | Turbulence dynamics in separated flows: the generalised Kolmogorov equation for inhomogeneous anisotropic conditions[END_REF] B t pδφ 1 q 2 `∇r ¨Φr `∇X ¨ΦX " Π r `ΠX (2.3)

The quantity under consideration is here the second-order moment (or structure function) of the increment (the difference) δφ 1 " φ 1 px `q´φ 1 px ´q where φ 1 " φ´φ is the fluctuating component of the LVF-field and the overbar denotes a suitable average. pδφ 1 q 2 is generally interpreted as the (cumulative) 'energy' content of the LVF fluctuations at a given position in the flow and at a given scale, This quantity depends a priori on time, on the position in the flow (more precisely the mid-point X between x `and x ´, see Fig. 1) and on the separation vector r " x `´x ´. In addition to a non-stationary term, Eq. (2.3) reveals two flux terms, in the r-and X-space, respectively. These are expressed as follows Φ r " pδuq pδφ 1 q 2 " pδuqpδφ 1 q 2 `pδu 1 q pδφ 1 q 2

(2.4a)

Φ X " pσuq pδφ 1 q 2 " pσuqpδφ 1 q 2 `pσu 1 q pδφ 1 q 2 (2.4b)
where the symbol σβ " pβ ``β ´q{2. ∇ r and ∇ X denotes the gradient operator w.r.t the r-and X-space, respectively. Therefore, ∇ r ¨Φr and ∇ X ¨ΦX are characteristics of a transfer in the r-and X-space, respectively. As per [START_REF] Valente | The energy cascade in grid-generated non-equilibrium decaying turbulence[END_REF]; [START_REF] Gomes-Fernandes | The energy cascade in near-field non-homogeneous non-isotropic turbulence[END_REF], the term Φ r represents the flux of the fluctuating quantity pδφ 1 q 2 from spherical shells centred at X with radius r " |r| either to spherical shells centred at the same X but with different radius or within the same spherical shell but to a different orientation r{r. This flux operates with a characteristic velocity δu " δu `δu 1 , i.e. the sum of mean and fluctuating velocity increments. On the other hand, the term Φ X represents the flux of the fluctuating quantity pδφ 1 q 2 from spherical shells centred at X with radius r " |r| to a spherical shell of same radius centered at another point in flow-position space. The latter is driven by a characteristic velocity σu " σu `σu 1 , i.e. the sum of mean and fluctuating velocity half sum. The RHS of Eq. ( 2.3) further contains two terms which are interpreted as some production terms due to gradient of mean quantities in the r-and X-space, respectively Π r " ´2pδu 1 qpδφ 1 q ¨∇r pδφq (2.5a)

Π X " ´2pσu 1 qpδφ 1 q ¨∇X pδφq (2.5b)
Π r and Π X can further be seen as representative of the 'energy' transfer between the mean and fluctuating fields (see Appendix C).

Eq. ( 2.3) reveals that pδφ 1 q 2 changes in time due to (i) a production mechanism associated with statistical inhomogeneity of mean quantities, and (ii) redistribution through fluxes in both scale-space r and position space X. Contrary to the well known equation for the scalar field such as temperature or pollutant concentration (i.e. the generalised Yaglom equation, see e.g. [START_REF] Danaila | A generalization of Yaglom's equation which accounts for the large-scale forcing in heated decaying turbulence[END_REF]; [START_REF] Gauding | Generalised scale-by-scale energy-budget equations and large-eddy simulations of anisotropic scalar turbulence at various schmidt numbers[END_REF]; [START_REF] Togni | Physical and scale-by-scale analysis of Rayleigh-Bénard convection[END_REF]), the LVF field is a non diffusive scalar. Therefore, the diffusion term and the scalar dissipation rate do not appear in Eq. (2.3). This is a very noticeable peculiarity of the LVF, since it suggests that there is no dissipation mechanism of fluctuations of the liquid volume fraction.

Local and non-local interactions

As per [START_REF] Mollicone | Turbulence dynamics in separated flows: the generalised Kolmogorov equation for inhomogeneous anisotropic conditions[END_REF], it is convenient to introduce the flux convection velocities in either the r-or X-space w r " pδuq pδφ 1 q

2 pδφ 1 q 2 (2.6a) w X " pσuq pδφ 1 q 2 pδφ 1 q 2 (2.6b)
and define m as either r or X, so that w m denotes either w r or w X , which then allows to write Eq. ( 2.3) in a Lagrangian form, viz.

B t pδφ 1 q 2 `wm ¨∇m pδφ 1 q 2 " Π m ´pδφ 1 q 2 ∇ m ¨wm (2.7)
Note that Eq. (2.7) is not equivalent to Eq. ( 2.3) except if the X and r space are independent. In case of statistically homogeneous flows (which will be the case hereafter), gradients with respect to X are zero and Eq. (2.7) applies and is equivalent to Eq. ( 2.3).

The main advantage of writing the scale-by-scale budget in the form provided by Eq. (2.7) is that the transfer term is decomposed into a transport term (i.e. the rightmost term on LHS of Eq. (2.7)) with convection velocity w m whose sign indicates the direction of the flux (negative/positive when directed towards smaller/larger scales), and a sink/source term accounting for the non solenoidal character of w m (the RHS of Eq. (2.7)). Variations of pδφ 1 q 2 along the trajectories of the field w m then appear to be due to the combined effect of production Π m and source term associated with the non solenoidal character of w m .

To further elaborate on the origin of the source term pδφ 1 q 2 ∇ m ¨wm , it is worth applying an Helmholtz decomposition to w m denoting either w r or w X . The velocity w m can then be decomposed as the sum of a divergence-free (solenoidal) and a curl-free (irrotational) component. This writes as

w m " ´∇m P m `∇m ^Am (2.8)
where P m and A m are the scalar and vector potential, respectively in either the r or X space. The source term f m " ∇ m ¨wm can thus be expressed in the form of a Poisson equation, viz. ∇ m ¨wm " ´∇m ¨p∇ m P m q with m " r or X. Using Green's function for Poisson's equation on an unbounded space and assuming that f m decays faster than |m| ´1 at large m, the solution for the scalar potential is given by

P m " ż m 1 G m pm 1 ´mqf m pm 1 qdm 1 (2.9)
G m pm 1 ´mq " t4π|m 1 ´m|u ´1, the Newtonian potential, describes the response of the system at a position m in the manifold m " r or X to a point source located at a position m 1 . Writing the sink/source term in the form of a Poisson's equation therefore suggests that a given scale (or position) interact with all others scales (positions) through a scalar potential field. The sink/source term in Eq. (2.7) therefore contains information about the non-local character of the interactions between different scales or positions. A similar reasoning was elaborated as regards to the non-local effect of pressure in incompressible turbulence by e.g. [START_REF] Tsinober | An informal conceptual introduction to turbulence[END_REF], see p. 28. In other words, Eq. (2.7) can be seen as a way of decomposing the transfer terms into local and non-local processes. The rightmost term on LHS of Eq. (2.7) characterizes local interactions between two neighbouring scales r or positions X. Such interactions are acting through a local transport process. On the other hand, the rightmost term on RHS of Eq. (2.7) emphasizes that there exists a potential field P m which allow separated scales (or positions) to interact.

Further comments

It is worth recalling that Eq. (2.7) is exact (no assumptions were invoked except Eq. (2.1), see Appendix A). It applies to the entire flow field, from the injection to the spray dispersion zone and irrespectively of the flow configuration or regime, as long as the flow remains incompressible with no phase change. The notion of scale is explicit and unambiguously defined as the straight distance between the two-points considered. Consequently, it allows exploring notably the direction and amplitude of the transfer among the different scales, the effect of different physical parameters (surface tension, viscosity and density, inflow velocity conditions), and the range of scales over which these parameters have an influence. Consequently, transposing the two-point statistical analysis to φ appears particularly promising for providing insights into the complex physics of liquid transport and atomization.

However, Eq. ( 2.3) or (2.7) have argument list pX, r, tq, i.e. a 6+1 dimensional problem. This is the price to pay for the transport of liquid to be explored along the geometrical and scale space. For instance, this makes its graphical representation particularly tenuous, and one generally discuss the contribution of the different terms of Eq. ( 2.3) or (2.7) on some submanifolds (Mollicone et al. 2018, and references therein). This complexity can however be reduced given the possible symmetries of the problem and the averaging operators which are used. For instance, statistical homogeneity in a particular direction will cancel out the X gradient in this direction. Isotropy, i.e. invariance by rotation of the r separation vector will allow to write equations as a function of the separation modulus r " |r| only. Note also that this framework is only an exploratory tool, but is incapable of predictions as the sequential cascade [START_REF] Gorokhovski | Analyses of Kolmogorov's model of breakup and its application into Lagrangian computation of liquid sprays under air-blast atomization[END_REF]Kolmogorov 1941a), the aggregation [START_REF] Villermaux | Ligament-mediated spray formation[END_REF]Marmottant & Villermaux 2004) or the maximum entropy [START_REF] Sellens | A simplified prediction of droplet velocity distributions in a spray[END_REF][START_REF] Dumouchel | The maximum entropy formalism and the prediction of liquid spray dropsize distribution[END_REF]) phenomenologies can provide. The reason is that Eq. (2.3) is not closed, i.e. the transport of second-order statistics depends on higher (third-) order terms.

Numerical database description

Numerical setup

The theory is appraised using Direct Numerical Simulations. We use the High-Performance-Computing code ARCHER developed at the CORIA laboratory. It was one of the first code worldwide, undertaking the simulation of liquid-jet atomization under a realistic diesel injection configuration [START_REF] Ménard | Coupling level set/VOF/ghost fluid methods: Validation and application to 3D simulation of the primary break-up of a liquid jet[END_REF]. It solves on a Cartesian mesh the one-fluid formulation of the incompressible Navier-Stokes equation, viz.

B t ρu `∇ ¨pρu b uq " ´∇p `∇ ¨p2µDq `f `γκδ s n (3.1)
where p is the pressure field, D the strain rate tensor, f a source term, µ the dynamic viscosity, ρ the density, γ the surface tension, n the unit normal vector to the liquid-gas interface, κ its mean curvature and δ s is the Dirac function characterizing the locations of the liquid gas interface. For solving Eq. (3.1), the convective term is written in conservative form and solved using the improved Rudman (1998) technique presented in [START_REF] Vaudor | A consistent mass and momentum flux computation method for two phase flows. application to atomization process[END_REF]. The latter allows mass and momentum to be transported in a consistent manner thereby enabling flows with large liquid/gas density ratios to be The simulated configuration is that of a turbulent two-phase flows in a triply periodic box similar to [START_REF] Duret | DNS analysis of turbulent mixing in two-phase flows[END_REF][START_REF] Duret | Improving primary atomization modeling through DNS of two-phase flows[END_REF]. This academic configuration, widely explored in the context of single phase flows, is here extended to a two-phase case. It has already proved to be extremely beneficial for better characterizing some features of primary atomization. In the particular context of two-point statistics, this configuration was chosen because it is statistically homogeneous thereby allowing to explore the dynamics of liquid transport in scale space only.

The physical parameters are the same as [START_REF] Duret | DNS analysis of turbulent mixing in two-phase flows[END_REF][START_REF] Duret | Improving primary atomization modeling through DNS of two-phase flows[END_REF]: the density of liquid and gas phases are set to ρ L " 753 kg.m ´3 and ρ G " 25 kg.m ´3, respectively, leading to a density ratio of about 30. The dynamic viscosity for the liquid and gas phase are µ L " 5.65 ˆ10 ´4 Pa.s and µ G " 1.879 ˆ10 ´5 Pa.s respectively. Therefore, the kinematic viscosity is the same for both the liquid and gas phase. The surface tension is γ " 0.0135 N.m ´1 and the box is filled with φ 0 " 10% of liquid. The domain is L " 1.5 ˆ10 ´4 m wide in each direction and 256 3 grid points are used for discretization. Several validation tests for the appropriateness of this resolution was undertaken by [START_REF] Duret | DNS analysis of turbulent mixing in two-phase flows[END_REF][START_REF] Duret | Improving primary atomization modeling through DNS of two-phase flows[END_REF] and will be further discussed below.

The actual initialisation of such flow is somehow arbitrary. Possible options are to simulate an actual injector, but the cost of such simulation is first very important for turbulent liquid-gas flow and secondly will lead to inhomogeneous anisotropic flow that are more difficult to analyse. Another possibility is to create a turbulent single phase flows and then to plug on it some liquid elements. This approach has been tested but leads to shear stress across the surface that are unlikely. This leads to violent atomisation at small scale with non-physical production of small droplets that can hardly be treated numerically. Thus, to avoid these drawbacks, the strategy we employ is to start from a flow with large enough liquid structures and reduced turbulent intensity, then to apply a numerical forcing procedure to reach a given level of turbulence intensity (see [START_REF] Duret | DNS analysis of turbulent mixing in two-phase flows[END_REF][START_REF] Duret | Improving primary atomization modeling through DNS of two-phase flows[END_REF]). This ensures that the flow reaches a sort of equilibrium between kinetic energy and surface tension/viscous forces. This procedure, although time consuming, allows for a realistic initial turbulent liquid gas flow. The counter part of such an approach is the CPU cost that constrain our simulation to reduced mesh and Reynolds number by comparison with those achieve for single phase flow.

Here a value for the targeted turbulent kinetic energy q 2 " u 2 `v2 `w2 " 7.2 m 2 .s ´2 was chosen as in [START_REF] Duret | DNS analysis of turbulent mixing in two-phase flows[END_REF][START_REF] Duret | Improving primary atomization modeling through DNS of two-phase flows[END_REF]. This leads to a turbulent Reynolds number of 194, the Taylor microscale Reynolds number is 25, the turbulent Weber number W e t " ρ L q 2 L i {γ is 8.1. The ratio of the Taylor microscale and Kolmogorov length scale to the mesh cell size is 20.0 and 2.05 respectively. The forcing is kept activated for several tens of integral time scales. Finally the forcing is turned off and the turbulence is let free to decay. The analysis of two-point statistics is carried out during the decay period only. Typical snapshots of the decaying liquid-gas turbulence simulations are portrayed in Fig. 2. These are separated by about 5 integral time scale T i " L i { b q 2 (here T i is estimated at t " 0 and L i is the integral length-scale of the velocity field, inferred by integrating the velocity correlation function). At initial time (t " 0), i.e. when forcing is turned off, the liquid structures appear highly tortuous, composed of elongated ligaments, drops and bubbles. Then, the turbulence intensity starts to decrease as shown by the decrease of fluctuations of the velocity magnitude. Meanwhile, liquid structures progressively coalesce, and surface tension also starts to dominate over inertial forces. Consequently, as turbulence decays, liquid structures grow in size.

Mesh convergence

Unlike single-phase turbulence simulations, determining a priori the smallest scales to be resolved in DNS of turbulent liquid-gas flows remains yet unfeasible (Gorokhovski & Saveliev 2003; [START_REF] Scardovelli | Direct numerical simulation of free-surface and interfacial flow[END_REF]. In interface capturing methods such as the one considered here, break-up of liquid structures occur automatically when two interface segments are within the same computational grid cell. Therefore, numerical results are inherently grid dependent. As stated by [START_REF] Gorokhovski | Analyses of Kolmogorov's model of breakup and its application into Lagrangian computation of liquid sprays under air-blast atomization[END_REF], it is however expected that when resolution is sufficient, the error in liquid volume (9 the cell volume) becomes negligible compared to the volume of typical liquid structures.

Here, mesh convergence was assessed on the basis of two metrics, one related to the liquid gas interface, viz. the liquid-gas surface density (the total surface area divided by the simulation volume, noted Σ) and one related to the small scale features of the velocity field, i.e. the norm of the strain rate tensor |D|. Results are presented in Fig. 3 where N x " N y " N z , the number of points in x, y and z directions is varied between 32 and 256. One observes that both metric nicely converge to a plateau when resolution increases. |D| and Σ are here computed with the forcing active and time-averaged over about 20 integral time-scales. Σ varies of less than a percent while |D| increases of less than 7% when N x is increased from 192 to 256. A resolution of 256 3 is therefore considered as being appropriate. A 512 3 simulation might be necessary for confirming this, however a more stringent test will be discussed later in section 4.3 concluding that a 256 3 mesh size is nicely tailored. Hereafter, unless explicitly specified, only the results for this particular mesh size will be discussed. In addition, it is worth stressing that this test of mesh convergence pertain to the case where forcing is active. When forcing is turned off, velocity gradients decrease in magnitude while liquid structures grow in size (see Fig. 2) meaning that the numerical resolution is enhanced all along the decay.

Anisotropy

The flow is homogeneous in all three directions. Therefore, two-point statistics are spatially averaged. Given the periodic boundary conditions, the X-terms in Eq. (2.3) vanish, together with the r-production term (see Appendix B). Consequently, in a triply periodic turbulent two-phase flow, the budget for the total LVF field is the same as the fluctuating field (see Appendix C). Further, to improve statistical relevance of our analysis, we have performed 28 simulations of the decaying flow field starting with different initial conditions, all separated by at least one integral time. In addition to the spatial average, two-point statistics were thus further averaged over these 28 independent realizations. Increments are calculated for separations in the range ´N x{2 ď pr x , r y , r z q ď N x{2. Spatially and ensemble averaged statistics are denoted x'y R,E .

Despite this quite large number of samples, two-point statistics were found to remain substantially anisotropic. Isotropy implies that iso-contours of xpδφq 2 y R,E in different Results in different planes are shown, full lines: pr x , r y , 0q plane, dashed lines: pr x , 0, r z q plane, dash-dotted lines: p0, r y , r z q plane. planes passing by the origin, should form concentric circles. Here, Fig. 4 emphasizes that second-order statistics of the liquid volume fraction strongly depend on the orientation of the separation vector r. Although not shown, it also appeared that the degree of anisotropy persists as turbulence decays. To a large extent, this is the result of the lack of statistical convergence of two-point statistics when described in a 4D pr, tq manifold. Indeed, as seen in Fig. 2, the number of liquid structures over which statistics are calculated is somewhat limited. For instance, it is observed that the domain contains about a dozen of liquid structures in Fig. 2(a) finally coalescing to few larger drops in Fig. 2(d). Therefore, even with 28 statistically independent runs, there are not sufficient samples for liquid structures to be isotropically oriented in space and hence for two-point statistics to loose their dependence to the orientation of the vector r.

Statistical convergence

Because of this substantial anisotropy which results from the lack of statistical convergence, we decided to further average two-point statistics over all orientations of the separation vector r. This is achieved by applying a sphere average operation, noted x'y S , as described in Appendix B. Therefore, the two-point budget of the LVF field in decaying turbulence reduces to rotation in r-space). For instance, Thiesset et al. (2019a) provide the two-point budget for a time-evolving liquid-gas shear layer.

B t A pδφq 2 E R,E,S
• r r (δφ) 2 R,E,∂S (b) r = 2dx r = 22dx r = 42dx r = 62dx r = 82dx r = 102dx r = 122dx
In order to check the statistical convergence of our results, we plot xpδφq 2 y R,E,S and 3 r xpδuq ¨r r pδφq 2 y R,E,BS as a function of the number of simulation runs for different separations r (Fig. 5). One observes that while second-order statistics (Fig. 5 (a)) converge rather quickly (« 0.1% after 28 runs), convergence of third-order statistics (Fig. 5 (b)) is much slower and it was necessary to run 28 simulations for the transfer term to be within 5%. This is considered as being sufficient for concluding that statistics have converged.

Running 28 simulations with the ARCHER code corresponds to about 100,000 CPU hours and calculating two-point statistics over about 100 snapshots for each of the 28 runs, represent to about 500 000 CPU hours. In total, the database ends up with 10TB of files. For this reason, it was particularly onerous to address another set of physical parameters (surface tension, total liquid fraction, kinetic energy, etc) and this aspect is left for future investigations.

Results

We now aim at evaluating the potential of two-point statistical analysis of the LVF field using the numerical database described above. This section addresses three main aspects. An analysis of asymptotic scaling of xpδφq 2 y E,R,S when the separation r tends to either large or small values is first performed. This will allow us to emphasize the existence of a characteristic length scale which is closely related to the liquid-gas interface surface density. Next, a similarity analysis is undertaken in order to infer the similarity variables which ought to be used in the scaling of the scale transfer. Finally, light will be shed on the contribution of the different terms of the scale-by-scale budget during the decay of homogeneous liquid-gas turbulence.

Asymptoptic scalings

We start by analysing the asymptotic behaviour of xpδφq 2 y E,R,S in the limit of small and large separations. To this end, it is worth recalling that the LVF field is rather similar to a so-called asymmetric random telegraphic signal, which features sharp transitions between a value of 1 (corresponding here to the liquid) and 0 (for the gas phase).

Consider a 1D asymmetric random telegraphic signal φ T pxq between 0 and 1, where the probability ppφ T pxq " 1q " φ 0 and ppφ T pxq " 0q " 1 ´φ0 . In our case, φ 0 " 0.1, i.e. the box is filled with 10% of liquid. The autocorrelation of such a signal, i.e. the average of φ T pxqφ T px `rq, can be calculated analytically [START_REF] Machlup | Noise in semiconductors: spectrum of a two-parameter random signal[END_REF][START_REF] Fitzhugh | Statistical properties of the asymmetric random telegraph signal, with applications to single-channel analysis[END_REF], and is given by

R φ T prq " φ 0 " φ 0 `p1 ´φ0 q e ´r L ‰ (4.1)
where L is a characteristic length scale which relates to the random jump frequency, the probability of transiting from a value of 1 to 0 and vice-versa. The mean squared value of φ T is R φ T pr " 0q " φ 0 and at large separations, the auto-correlation tends to R φ T pr Ñ 8q " φ 2 0 . The fact that the auto-correlation of a telegraphic signal does not tend to 0 when r Ñ 8 precludes defining an integral length scale. The second-order moment of the increment δφ T , can be expressed in terms of the autocorrelation function, i.e. pδφ T q 2 prq " 2 pR φ T p0q ´Rφ T prqq " 2 pφ 0 ´Rφ T q (4.2)

where again the overbar denotes a suitable average (here a spatial average over the x axis). Therefore, for a telegraphic signal, we have

pδφ T q 2 prq " 2φ 0 p1 ´φ0 q " 1 ´e´r L ‰ (4.3)
Consequently, at large separations pδφ T q 2 prq tends to a constant equal to 2φ 0 p1 ´φ0 q.

The limit towards small separations is lim

rÑ0 pδφ T q 2 prq " 2φ 0 p1 ´φ0 q r L `O ˆr2 L 2 ˙(4.4)
Eq. ( 4.4) indicates that the second-order structure function of a telegraphic signal φ T is proportional to r when r is small compared to L. Eq. (4.4) allows calculating the characteristic length scale L as

L " 2φ 0 p1 ´φ0 q « lim rÑ0 pδφ T q 2 prq r ff ´1 (4.5)
The comparison between xpδφq 2 y E,R,S estimated from the DNS and pδφ T q 2 is portrayed in Fig. 6(a). Results are shown for an arbitrary time t " 4.23T i but similar behaviours were observed for all t. In Fig. 6, xpδφq 2 y E,R,S and pδφ T q 2 are normalized by 2φ 0 p1 ´φ0 q.

At large scales, we then observe that both curves tend to an asymptotic value of 1 consistently with the aforementioned considerations. As a consequence, the limit at large r can be used to compute the total volume occupied by the liquid φ 0 . At small separations (r{L ď 0.1), xpδφq 2 y E,R,S approaches pδφ T q 2 and appears to be nicely proportional to r{L. Therefore, the asymptotic scaling at large and small scales is in perfect agreement with an asymmetric telegraphic signal. It is further worth noting that, at small scales, the scaling xpδφq 2 y E,R,S 9 r{L extends up to scales r « L{2. The fact that at small scales, xpδφq 2 y E,R,S 9 r{L precludes defining a Taylor microscale for the LVF field since the latter assumes a r 2 scaling of second-order structure functions at small scales. Therefore, by analogy with the asymmetric random telegraphic signal, it appears that neither an integral length nor a Taylor length-scale can be defined and that the only length-scale which one can define is L as given by Eq (4.5).

As mentioned earlier, L relates to random jump frequency. Hence, there should exist a link between L and a length scale based on the liquid-gas surface area. [START_REF] Anez | Eulerian-Lagrangian spray atomization model coupled with interface capturing method for diesel injectors[END_REF] proposed a slight modification of the interface length-scale firstly introduced by Lebas

10 -1 10 0 r/L 10 -1 10 0 (δφ) 2 E,R,S / 2φ 0 (1 -φ 0 ) (a) (δφ) 2 E,R,S (
L Σ " 6φ 0 p1 ´φ0 q Σ (4.6)
where Σ is the surface density (the surface area divided by total volume) of the liquid-gas interface. L Σ can be interpreted as a generalization of Sauter Mean Diameter of nonspherical objects, hence the pre-factor 6. A rather similar definition for L Σ was given by [START_REF] Peters | A spectral closure for premixed turbulent combustion in the flamelet regime[END_REF], see Eq. (5.2)) for describing turbulent premixed flames using a level-set formulation. By analysing the DNS data, we found out that

Σ " 3 lim rÑ0 xpδφq 2 y E,R r (4.7)
Eq. (4.7) can most likely be proven analytically knowing that Σ " x|∇φ|y R . However, we were not able to end up with a rigorous demonstration for general cases. In the scientific community of porous media or composite materials, the phase indicator function is often used as a scalar representing the different phases composing the material. It thus has some degrees of ressemblance with the LVF field. For instance, Eq. (4.1) is often used to characterize the structure of particular type of porous media, known as the Debye random media [START_REF] Debye | Scattering by an inhomogeneous solid[END_REF][START_REF] Debye | Scattering by an inhomogeneous solid. ii. the correlation function and its application[END_REF]. By scrutinizing this litterature, we discovered that [START_REF] Yeong | Reconstructing random media[END_REF] derived the same expression as Eq. (4.7) (see their Eq. ( 8)), but again their methodology does address very general cases. [START_REF] Debye | Scattering by an inhomogeneous solid. ii. the correlation function and its application[END_REF] proposed a similar relation as Eq. (4.7), except that they found a prefactor of 8/3 instead of 3. At this stage, we thus have no other solutions than numerically validating Eq. (4.7).

Comparison of Σ as estimated directly by computing the surface area of the zero levelset and the values obtained from Eq. (4.7) is shown in Fig. 6(b). The two curves are almost undistinguishable from each other indicating that, although not demonstrated mathematically, Eq. (4.7) is correct. We have also tested that Eq. (4.7) remains valid when φ 0 is varied (see Appendix D). Consequently, the values for L Σ estimated either from the surface density or from the limit of xpδφq 2 y E,R at small scales are also in perfect agreement (see Fig. 6(b)). Eq. (4.7) yields L " L Σ . In other words, the only relevant length-scale that emerges for the LVF field is L Σ as defined by either Eq. (4.5) or Eq. (4.6).

To sum up, the asymmetric telegraphic signal is a nice surrogate of the LVF field for analysing the asymptotic scaling of xpδφq 2 y E,R at either large and small scales. We have proven that xpδφq 2 y E,R contains informations about the geometrical nature of the liquid-gas interface (its surface area) and the total liquid fraction φ 0 . Nevertheless, the distribution pδφ T q 2 departs significantly from xpδφq 2 y E,R at intermediate scales which indicates that the LVF field is not a purely random field but has certain features representative of the morphology of the liquid structures. Comparing xpδφq 2 y E,R and pδφ T q 2 could possibly be exploited in future investigations to shed light on the peculiar geometrical facets of the liquid-gas interface.

Similarity analysis

Finding which similarity variables are the most appropriate for normalizing two-points statistics provides insights into the prominent mechanisms at play in a given region of flow and over a given range of scales. Potentially, self-similarity can further reduce scaleby-scale budgets (i.e. partial differential equations) to ordinary differential equations, thereby leading to predictions of one-point statistics and simple closures of the scale-byscale budget.

In the light of the methodology developed for single phase turbulence, we apply a similarity analysis to the budget Eq. (3.2). To this end, we start by writing second-and crossed moments of LVF increments in the form

A pδφq 2 E R,E,S " epξqE ˚ptq (4.8a) 3 r A pδuq ¨r r pδφq 2 E R,E,BS " gpξqG ˚ptq (4.8b)
with ξ " r{λptq where λptq is the similarity length scale which depends on time but is a priori unknown. E ˚ptq and G ˚ptq are two normalizing quantities which pertain to second-and crossed moments, respectively. As was done by e.g. George (1992a); [START_REF] Antonia | Similarity of energy structure functions in decaying homogeneous isotropic turbulence[END_REF]; Burattini et al. (2005b), second-order structure functions are generally normalized by their value when the separation tend to large values. Here, we proceed similarly and assume

E ˚ptq " lim rÑ8 A pδφq 2 E R,E,S
" 2φ 0 p1 ´φ0 q (4.9)

It is therefore assumed that E ˚" const. Injecting Eqs. (4.8) into Eq. (3.2), yields

B t epξqE ˚ptq " ´gpξqG ˚ptq (4.10)
Assuming E ˚" const., we have

E ˚ptqtB ξ epξqutB t ξptqu " ´gpξqG ˚ptq (4.11)
recalling that B t ξptq " ξB t log λ, we end up with rE ˚Bt log λs ξB ξ epξq " gpξq rG ˚ptqs (4.12)

Self-similarity implies that the terms which depend only on time (the terms within brackets) should behave similarly w.r.t time, viz.

G ˚ptq 9 B t log λptq (4.13) However, it is not revealed which length-scale λptq should be used for normalising the

-Transfer r / ∂ t log(L Σ ) (b) (b) (b) (b) (b) (b) (b) (b) (b) (b) (b) (b) (b) (b) (b) (b) (b) (b) (b) (b) (b)
Figure 7: Assessement of the degree of self-similarity of (a) second-and (b) crossed moments when the separation r is normalized by L Σ . The different colours from dark to light correspond to increasing time during the decay. separation r. Given the asymptotic scaling described in the previous section, it worth testing L Σ as the similarity length scale. Then, it is expected that

G ˚ptq 9 ´Bt log Σ " K Σ (4.14)
If L Σ were to be used as a similarity length-scale, then the transfer term should scale with the stretch rate K Σ of the liquid-gas interface.

The adequacy of the proposed similarity variables to normalize the scale distribution and scale transfer of the LVF is presented in Fig. 7(a) and (b), respectively. The scale distributions xpδφq 2 y R,E,S collapse very nicely when the separation r is normalized by L Σ . However, self-similarity is not complete, i.e. it does not span the entire range of scales, but appears to be accurately satisfied up to r « L Σ {2. This means that L Σ (and thus Σ) is the relevant quantity for characterizing the time evolution of xpδφq 2 y R,E,S , at least at small scales. As far as the transfer term in Eq. (3.2) is concerned, i.e. Transfer r " ´3{rxpδuq ¨r{r pδφq 2 y R,E,BS (Fig. 7(b)), although the scatter is significant especially at large scales (r Á L Σ {2), it appears that again L Σ and thus K Σ is a relevant similarity variable. In other words, time variations of xpδφq 2 y R,E,S for all scales r À L Σ {2 are uniquely determined by the stretch rate. Consequently, in analogy with the Kolmogorov or Batchelor scaling of second-order structure functions for either the velocity or scalar field, the stretch rate plays here the same role as the kinetic energy or scalar dissipation rate. This statement can be further illustrated by taking the limit at small separations of the scale-by-scale budget. Plugging Eq. (4.7) into Eq. (3.2) yields

´3 r A pδuq ¨r r pδφq 2 E R,E,BS " 2φ 0 p1 ´φ0 q K Σ r L Σ , (4.15)
for all r ! L Σ . Eq. (4.15) can be viewed as the equivalent for a non diffusive scalar such as the LVF, of the Kolmogorov (1941a) four-fifth law or the [START_REF] Yaglom | On the local structure of a temperature field in a turbulent flow[END_REF] fourthird law. It suggests that for scales sufficiently small compared to L Σ , the transfer term scales with r{L Σ and is proportional to the stretch rate. This leads us to the startling conclusion that the equivalent for a non diffusive scalar of the scalar dissipation rate is the stretch rate of the interface separating φ " 0 and φ " 1.

The relevance of the stretch rate as a similarity variable for the energy transfer can be illustrated as in Fig. 8. Consider the case where an interface is present between the two-points x `and x ´, meaning that δφ ‰ 0. Two of such cases are represented in Fig. 8 and are labelled (1) and ( 2). Here, pδuq ¨rpδφq 2 is non zero if pδuq ¨r is non zero, i.e. there exists a velocity gradient along the direction r. When pδuq ¨r ą 0 (Case 1), the liquid structure undergoes a stretching mechanism along the direction r leading to an increase in size in this particular direction (this is illustrated by the green patch in Fig. 8). In contrary, when pδuq ¨r ă 0 (Case 2), the liquid structure is submitted to a compression effect and its size tends to decrease (as illustrated by the green patch in Fig. 8). Therefore, this simple sketch illustrates that the mechanism behind the scale energy transfer is closely related to the stretch/compression rate acting at some scale r. The same reasoning was often stated for the turbulent cascade as being caused by vortex stretching (Lundgren 1982, for instance).

We note that the turbulence characteristics (i.e. the kinetic energy, the integral and Taylor length-scale) do not appear as plausible similarity variables in the scale transport of the LVF field. This proves again that the effect of turbulence is implicit in this theoretical framework. Turbulence shows itself implicitly through its effect on the liquid gas interface and more specifically on the surface density and stretch rate.

LVF Budget and flux decomposition

We now turn our attention to the contribution of each term of the scale-by-scale budget Eq. (3.2). We start by verifying that the budget given by Eq. (3.2) is accurately closed. Results are presented in Figs. 9(a-c) for three arbitrary time values during the decay. There were not noticeable differences for other time values. It appears that the time derivative term is perfectly balanced by the transfer term. Because Eq. (3.2) is exact, this represents a very stringent test of the accuracy of the simulation results. In addition to the mesh convergence presented in section 3.2, it thus reinforces that the spatial resolution used for the present simulations is adequate. Figs. 9(a-c) also present the profiles of the transfer term as calculated from the similarity solution gpξq " ´ξB ξ epξq. The latter appears to be well suited for describing the scale transfer at small scales (r À L Σ {2), i.e. over the range of scales where similarity applies.

It is worth noting that in our case, the r-transfer term in Eq. (3.2), i.e. 9 ´xpδuq r pδφq 2 y R,E,BS is negative irrespectively of r. This indicates that, on average, liquid structures are growing in size. This is consistent with Figs. 2, where it was observed that initially stretched structures (mostly ligaments, see Fig. 2(a)) relax to a more 'spherical' objects (Fig. 2(d)) due to coalescence and surface tension effects. Therefore, for the LVF to be conserved, liquid structures which are longitudinally decreasing in size, concomitantly expands in the ligament transverse directions. Averaging over all 4.16). (g-i) Decomposition of the r-Advection term using a potential and solenoidal velocity Eq. (4.17). Results are presented for t " 4.69, 6.57, 10.33T i . directions r leads to an increase in size, viz. on average the transfer of LVF is directed towards large scales. This is a very peculiar difference with a diffusive passive scalar, for which the transfer term has opposite sign (i.e. is directed towards small scales) and diffusion is the driving mechanism leading, on average, to an increase in size of turbulent scales [START_REF] Danaila | A generalization of Yaglom's equation which accounts for the large-scale forcing in heated decaying turbulence[END_REF][START_REF] Antonia | Similarity of decaying isotropic turbulence with a passive scalar[END_REF]Gauding et al. 2014, for instance). In case of turbulent liquid-gas flows, the LVF does not diffuse, and thus surface tension and coalescence are the only mechanisms responsible for the relaxation of elongated ligaments into larger (on average) liquid structures.

We now turn our attention to the decomposition of the transfer term into local and non-local interactions. Spherically averaging Eq. (2.7) allows writing the transfer term in the form

´Apδuq ¨r r pδφq 2 E E,R,BS loooooooooooooomoooooooooooooon Transfer r " ´Aw r ¨∇ @ pδφq 2 D E,R E S looooooooooooooomooooooooooooooon Advection r ´A@ pδφq 2 D E,R ∇ ¨wr E S looooooooooooomooooooooooooon Non Solenoidal r (4.16)
where w r is given by Eq. (2.6a). As demonstrated in section 2.3, the decomposition Eq. (4.16) is a way of decomposing the transfer term into a local purely convective term and a source term associated with non local interactions in scale space. The contribution of these two mechanisms w.r.t the separation r is presented in Figs. 9(d-f) for the same three time values. Globally, the amplitude of the local and non-local processes have the same order of magnitude. The non-local process peaks at larger scales than the advection mechanism. At small scales (r À 0.3L Σ ), the advection dominates slightly over the nonlocal process. The contrary is observed at rather large scales r Á 0.3L Σ . The two-point analysis indicates a noticeable contribution of non local interactions in scale space, i.e. two separated scales can interact although they are very different in size. We explain this observation by the effect of surface tension. For illustrating this statement, consider the case of the Plateau-Rayleigh instability of a liquid ligament. Surface tension forces depend on the local mean curvature. Along the ligament axis, pressure is thus high at the crests (where curvature is negative) and low at the troughs (positive curvature). This effect is known as the capillary bridge. Surface tension forces will then tend to amplify the contraction (expansion) effect at the neck (crests), ultimately leading to pinch-off, break-up, followed by the retraction of the neck and the formation of one or several droplets. Before break-up, the size of neck and that of the forming droplet can be very different. However, they strongly interact through the capillary bridge, this interaction being mediated by the pressure field, an essentially nonlocal process. Therefore, this simple example suggests that, because of surface tension effects, non-local interactions can play an essential role in the scale transport of the LVF.

We further apply an Helmholtz decomposition to the convection velocity w r , thereby allowing to rewrite the r-Advection term in the form

´Aw r ¨∇ @ pδφq 2 D E,R E S looooooooooooooomooooooooooooooon Advection r " ´Aw pot r ¨∇ @ pδφq 2 D E,R E S loooooooooooooooomoooooooooooooooon Potential r ´Aw sol r ¨∇ @ pδφq 2 D E,R E S loooooooooooooooomoooooooooooooooon Solenoidal r (4.17)
where w pot r and w sol r are the potential (curl-free) and solenoidal (divergence-free) components, respectively. For doing this, we solve the Poisson equation ∇ r ¨wr " ´∇r ¨p∇ r P r q (4.18) and compute w pot r " ´∇r P r and w sol r " w r ´wpot r . Eq. (4.18) is solved using the MINRES algorithm implemented in the Scipy Library of Python. Given the symetry of the problem, zero normal gradient boundary conditions were used. Results are presented in Figs. 9(g-i). The overall contribution of the potential and solenoidal velocity to the advection term is rather similar. At small scales (r À 0.3L Σ ), they contribute almost identically while at large scales (r Á 0.3L Σ ), the potential velocity predominates. It would be speculative to elaborate further on the physical interpretation of the contributions of the potential and rotational velocities. Our intuition is that this decomposition contains informations about the anisotropic transfer between different orientations of the vector r.

Here, the spherical average hides the dependence on the orientation vector and anisotropy is mostly the result of the lack of statistical convergence. For providing further insights into the relation of these two terms to anisotropy, one should explore a flow which features a preferential anisotropy direction. In this context, homogeneous shear liquid-gas turbulence as simulated by [START_REF] Rosti | Emulsions in homogeneous shear turbulence[END_REF][START_REF] Rosti | Droplets in homogeneous shear turbulence[END_REF]) can be instructive.

Conclusion

Two-phase flows reveal complex multi-scale features which require specific theoretical treatments. Here, we propose using a two-point statistical analysis and adapt it to the specific context of turbulent liquid-gas flows. It relies on the transport equation for the liquid volume fraction which is written at two-points arbitrary separated by a distance r. By doing so, the transport of liquid is shown to depend on a transfer and a production process which act together and concomitantly in physical space (i.e. from one position in the flow to the other) and in scale space (e.g. from large to small or vice-versa). The theory further points out the co-existence of local and non-local interactions between different scales and positions. This framework is exact as long as the flow remains incompressible with no-phase change. It is therefore a nicely tailored tool to appraise the appropriateness of a given numerical simulation method and resolution. The notion of scale is explicit and unambiguously defined. It applies to the entire flow field, from the injection to the spray dispersion zone and irrespectively of the flow configuration or regime. The effect of different physical parameters (surface tension, viscosity and density, inflow velocity conditions), although implicit in the equations through their effect on the velocity field, can be probed as a function of the scale and flow position. Consequently, transposing the two-point statistical analysis to the LVF appears particularly promising for providing insights into the complex physics of liquid transport and atomization. It is important recalling that this framework, unlike the sequential cascade, the aggregation or the maximum entropy models, is only an exploratory tool, but is incapable of predictions. The reason is that Eq. (2.3) is not closed. We have however great hope that in the long run, it will spark new ideas for modelling in particular the subgrid scale transport of liquid in a LES context. Here, we invoke the assumption of self-preservation to end up with reasonable predictions of the evolution of the LVF field in scale space.

Light is shed on the dynamics of scale transport in the particular configuration of decaying liquid-gas turbulence. Attention has been paid to the appropriateness of the numerical database by undertaking a careful analysis of the convergence of the results with respect to the mesh resolution and the number of samples used to calculate the statistics. We have performed 28 simulations starting with statistically independent initial conditions. This ensures that spatially, ensemble and spherically averaged statistics have converged. However, a substantial degree of anisotropy was found, which precluded using the isotropic version of the scale-by-scale budgets. The computational resources for these 28 simulations were shown to be so onerous that only one set of physical parameters was studied. Other specific investigations will be carried out in next future to confirm the robustness of our conclusions when varying the surface tension, the density ratio, the turbulence intensity and the total liquid fraction in the domain. Throughout this paper, the results were derived from rather general considerations. It is therefore likely that our conclusions will remain unchanged, at least in this simple particular flow configuration. There remains however an open question on whether our findings will be affected by statistical inhomogeneity (e.g. velocity shear), i.e. a 'non-equilibrium turbulence' configuration (Vassilicos 2015;[START_REF] Valente | The energy cascade in grid-generated non-equilibrium decaying turbulence[END_REF].

One important outcome on the present study is that the LVF field highlights a nice similitude with an asymmetric random telegraphic signal. Their respective scale distribution was found to share the same scaling at both large and small scales. It was however observed that they differ significantly at intermediate scales. The telegraphic signal was also particularly insightful to demonstrate the mathematical impossibility of defining neither an integral length-scale nor a Taylor microscale. There exists only one relevant length-scale which is proved to be L Σ . The scale distribution pδφq 2 is thus a metric of the interface geometrical properties and in particular, its surface area. In future work, the difference between the scale distribution of the LVF field and that of a telegraphic signal can possibly be explored in order to extract some informations about the peculiar geometry of liquid structures. Alternatively, it could be worth analysing or modelling the LVF field in the light of the vast scientific literature pertaining to porous media and/or composite materials, where correlation functions of the phase indicator and their relations to the surface density and volume fraction have been extensively explored [START_REF] Debye | Scattering by an inhomogeneous solid[END_REF][START_REF] Debye | Scattering by an inhomogeneous solid. ii. the correlation function and its application[END_REF][START_REF] Yeong | Reconstructing random media[END_REF].

A second important result has emerged from the application of a similarity analysis to The spatial average commutes with r and t differential operators and with the ensemble average but not with X differential operators. Given any vector q, the divergence theorem implies

x∇ X ¨qy R " 1 V ¡ R p∇ X ¨qq dX " S V ¨1 S £ BR q ¨ndS '" S V xq ¨ny BR (B 3)
where S is the area of the surface BR bounding R, dS the differential of surface area and n the unit normal vector to the surface BR pointing outwards. Therefore, the spatial average of Eq. (A 7) writes Let us denote the spherical average of β over a sphere S of radius r by xβy S , defined by [START_REF] Hill | Exact second-order structure-function relationships[END_REF][START_REF] Nie | A note on third-order structure functions in turbulence[END_REF] xβy S pX, r, tq " ˆ3 4πr 3 ˙¡ S β pr, X, tq dr (B 6)

Again, the spherical average commutes with X and t differential operators and with the ensemble average but not with r differential operators. Similarly, for any vector q x∇ r ¨qy S " ˆ3 4πr 3 ˙¡ S p∇ r ¨qq dr " 3 r

$ & % 1 4πr 2 £ BS q ¨r |r| dS r
, .

-

" 3 r A q ¨r r E BS (B 7)
where x'y BS denotes the average over the shell of the sphere (i.e. the term between braces in Eq. (B 7)). Thus, the spherical average of Eq. (A 7) writes

Figure 1 :

 1 Figure1: Example of LVF field. Dark regions correspond to the liquid-phase and white regions to the gas-phase The interface is also displayed as black curves. The two-points x `and x ´used to calculate increments are also represented together with the midpoint X and the separation vector r.

Figure 2 :

 2 Figure2: Typical snapshots of decaying liquid-gas turbulence. The liquid-gas interface (the grey surface) corresponds to the zero-level set value. Also portrayed is the velocity magnitude. (a) t " 0, (b) t « 5T i , (c) t « 10T i , (d) t « 15T i . From (a) to (d), the turbulent Weber number W e t " ρ L q 2 L i {γ is equal to 8.1, 2.4, 1.1 and 0.7

Figure 3 :

 3 Figure 3: Mesh convergence of the liquid-gas surface density and velocity strain tensor.

Figure 4 :

 4 Figure 4: Contours of xpδφq 2 y R,E .Results in different planes are shown, full lines: pr x , r y , 0q plane, dashed lines: pr x , 0, r z q plane, dash-dotted lines: p0, r y , r z q plane.

  pr, tq, i.e. a two-dimensional problem. In Eq. (3.2), the divergence theorem was used to express the spherically average divergence of the r-flux in terms of the normal flux at the sphere boundary BS (see Appendix B). Note that Eq. (3.2) is a simplification of Eq. (2.3) which applies only to decaying homogeneous flows, the flow under consideration here. For other flow instances, one should refer to the general equation Eq. (2.3) and proceed to the relevant simplifications depending on the flow symmetries (invariance w.r.t time, invariance by translation in X-space or invariance by

Figure 5 :

 5 Figure 5: Statistical convergence of second (a) and third order (b) statistics for different separation r, as a function of the number of independent simulations used to calculate ensemble averages. Results are for t « 7.5T i .

Figure 6 :

 6 Figure6: (a) Comparison of xpδφq 2 y at t " 4.23T i with the second-order structure function of a telegraphic signal with same L. The asymptotic behaviour at small scales (i.e. r{L) and large scales (" 1) are also shown as dashdotted lines. (b) Comparison of Σ and L Σ (lines) estimated directly (from the surface area of the zero level-set) with those obtained from the limit of xpδφq 2 y E,R,S when r Ñ 0 as given by Eqs. (4.6) and (4.7) (symbols)

  (a) (a) (a) (a) (a) (a) (a) (a) (a) (a) (a) (a) (a) (a) (a) (a) (a) (a) (a) (a)

  Schematic representation of pδuqpδφq 2 .

Figure 9 :

 9 Figure 9: (a-c) Contribution of the different terms of the budget Eq. (3.2). The similarity solution gpξq " ´ξB ξ epξq is also displayed as dotted lines. (d-f) Decomposition of the r-transfer term into a local and a non-local transfer term Eq. (4.16). (g-i) Decomposition of the r-Advection term using a potential and solenoidal velocity Eq. (4.17). Results are presented for t " 4.69, 6.57, 10.33T i .

  ensemble average does not change the problem complexity and statistics have argument list pX, r, tq, i.e. a 7D problem. B.2. Spatial average: X-space average Let us denote the spatial average of β over a domain R of volume V in X-space by xβy R , i.e. xβy R pr, R, tq "

  have argument list pR, r, tq. Triply periodic boundary conditions are sometimes used in direct numerical simulations of turbulent flows. For such cases, the average value of the flux at the boundary BR cancel out, viz. problem to a 4D problem, i.e. function of pr, tq. B.3. Spherical average: r-space average

  to a 5D problem with argument list pX, r, tq

: Email address for correspondence: thiesset@coria.fr

Acknowledgements

We are thankful to A. Poux, research engineer at CORIA laboratory who optimized the post-processing code for calculating the terms of LVF budget and who wrote the Poisson solver. Computations have been carried out in CRIANN (Centre Régional Informatique et d'Applications Numériques de Normandie) under the project 2018002 and 2003008. FT also benefited from the financial support from the INSIS institute of the CNRS and the CORIA laboratory which are gratefully acknowledged.

the LVF scale-by-scale budget. In particular, it was shown that second-order structure functions of the LVF field scale very nicely with L Σ implying that the stretch rate of the liquid-gas interface is the relevant similarity variable for the r-transfer term. Self-preservation is however not complete since it holds only at small scales. This suggests that another set of similarity variables, which possibly pertain to the turbulence characteristics, should be used for the large scales. This remains an open question. We have further derived the equivalent of the Kolmogorov four-fifth law or the Yaglom fourthird law to the LVF field. In this relation, K Σ plays the same as the kinetic energy (or scalar) dissipation rate for the small scales of fluid turbulence (or scalar mixing). This is a very nice finding since it has important implications for modelling purposes. Indeed, as for the turbulent cascade in single phase flows, this suggests using the concept of a constant rate of 'energy' transfer between the different scales of the LVF field, this constant being the stretch rate. This further explains the adequacy of the scale similarity LES model for computing the subgrid transport of the LVF field, see e.g. [START_REF] Chesnel | Subgrid analysis of liquid jet atomization[END_REF]. Our study thus proves that the challenge in LES modelling of the LVF subgrid transport is thus pushed further on the modelling of the stretch rate which is a priori unknown. This reinforces the pertinence of LES models with specific treatment of the surface density such as the Eulerian-Lagrangian Spray Atomization (ELSA) model [START_REF] Anez | Eulerian-Lagrangian spray atomization model coupled with interface capturing method for diesel injectors[END_REF]. As stated earlier, the effect of the physical parameters such as surface tension, viscosity and density do not appear explicitly in the two-point equations of the LVF. However, these parameters affect the stretch rate in a way that remains to be explored. In addition, it remains unclear whether this analogy between the scalar dissipation rate and the stretch rate remains valid in presence of e.g. velocity shear. This point should be addressed in future work by considering e.g. homogeneous shear liquid-gas turbulence [START_REF] Rosti | Emulsions in homogeneous shear turbulence[END_REF][START_REF] Rosti | Droplets in homogeneous shear turbulence[END_REF].

Finally, exploring the contribution of the different terms of the scale-by-scale budget emphasizes that the similarity assumption is reasonable for providing a simple closure for the transfer term. Second, it is observed that the flux of LVF in scale space is driven by both local and non-local interactions. Our intuition is that surface tension forces and especially the formation of capillary bridge is the main ingredient explaining that nonadjacent scales can interact. Further work is however needed to confirm this statement. In this perspective, the case of e.g. the Plateau-Rayleigh instability could be insightful.

Declaration of Interests

The authors report no conflict of interest.

Appendix A. Unaveraged equations

Two-point statistical equations are obtained following the lines of e.g. [START_REF] Hill | Exact second-order structure-function relationships[END_REF]; [START_REF] Danaila | Progress in studying small-scale turbulence using exact two-point equations[END_REF]. Eq. (2.1) is written at a point x `and a point x ´separated arbitrarily by a distance r (see Fig. 1)

where superscript plus (minus) denote the quantity at point x `(x ´). Taking the difference between Eq. (A 1a) and Eq. (A 1b) yields

By noting that for any quantity β, we have

´" 0, Eq. (A 2) can be rewritten as

where δβ " β `´β ´" βpx `, tq´βpx ´, tq is the increment (the difference) of the quantity β between the two points x `and x ´. Let X being the midpoint and r the separation vector between x `and x ´, viz.

X "

Using the chain rule, the gradient with respect to X and r are related to

Injecting Eqs. (A 5a) and (A 5b) into Eq. (A 3) leads to

B t pδφq `∇r ¨pδuqpδφq `∇X ¨pσuqpδφq " 0 (A 6)

where σβ " 1 2 pβ ``β ´q is the average of the quantity β between the two points x `and x ´. Multiplying Eq. (A 6) by δφ leads to B t pδφq 2 `∇r ¨pδuq pδφq 2 `∇X ¨pσuq pδφq 2 " 0 (A 7)

Appendix B. Averaged equations

Because the flows under consideration are turbulent, it is worth deriving the averaged forms of Eq. (A 7) using different statistical operators, i.e. ensemble, spatial and spherical averaging procedures. This allows possibly reducing the problem complexity and also express some terms of Eq. (A 7) is an easier way.

B.1. Ensemble average

Because the ensemble average is a summation, it commutes with differential operators. Therefore, the ensemble average of Eq. (A 7) writes

Appendix C. Equation for the mean and fluctuating components

Let us decompose φ and u into a mean (the type of average might depend on the symmetry of the flow, it is here noted by the overbar) and a fluctuating component, viz.

By definition, the averaged fluctuations u 1 , φ 1 " 0. The equations for `δφ ˘2 and pδφ 1 q 2 are obtained multiplying Eq. (A 6) by δφ and δφ 1 , respectively, and then proceed to average.

C.1. Equation for the mean component

Taking each term of Eq. (A 6) individually, multiplying them by δφ and proceeding to average yields,

so that the transport equation for `δφ ˘2 writes B t `δφ ˘2 `∇r ¨pδuq `δφ ˘2 `2pδφq∇ r ¨pδu 1 qpδφ 1 q `∇X ¨pσuq `δφ ˘2 `2pδφq∇ X ¨pσu 1 qpδφ 1 q " 0 (C 3)

C.2. Equation for the fluctuating component

Taking each term of Eq. (A 6) individually, multiplying them by δφ 1 and proceeding to average yields,

pδφ 1 q∇ X ¨pσu `σu 1 q `δφ `δφ 1 ˘" pδφ 1 q∇ X ¨ pσuqpδφ 1 q `pσu 1 qpδφ 1 q `pσu 1 qpδφq ( (C 4c) Therefore the budget for pδφ 1 q 2 writes B t pδφ 1 q 2 `∇r ¨´pδuqpδφ 1 q 2 `pδu 1 q pδφ 1 q 2 ¯`2pδφ 1 qpδu 1 q ¨∇r pδφq `∇X ¨´pσuqpδφ 1 q 2 `pσu 1 q pδφ 1 q 2 ¯`2pδφ 1 qpσu 1 q ¨∇X pδφq " 0 (C 5)

In order to highlight the relation between the equation for the mean Eq. (C 3) and that of the fluctuating counterpart, Eq. (C 5) can be rewritten by noting that triple correlations 2 " `δu `δu 1 ˘`δφ `δφ 1 ˘2

" pδuq `δφ ˘2 `pδuqpδφ 1 q 2 `pδu 1 q pδφ 1 q 2 `2pδφq pδu 1 qpδφ 1 q (C 6a) pσuq pδφq 2 " `σu `σu 1 ˘`δφ `δφ 1 ˘2

" pσuq `δφ ˘2 `pσuqpδφ 1 q 2 `pσu 1 q pδφ 1 q 2 `2pδφq pσu 1 qpδφ 1 q (C 6b) so that Eq. (C 5) writes B t `δφ 1 ˘2 `∇r ¨pδuq pδφq 2 ´∇r ¨pδuq `δφ ˘2 ´2pδφq∇ r ¨pδu 1 qpδφ 1 q `∇X ¨pσuq pδφq 2 ´∇X ¨pσuq `δφ ˘2 ´2pδφq∇ X ¨pσu 1 qpδφ 1 q " 0 (C 7)

It can then easily be shown that Eq. (C 3) + Eq. (C 7) = Eq. (B 1). It is also worth noting that the X-and r-gradient terms appearing in Eq. (C 3) are also present in Eq. (C 7) but with opposite sign. This indicates that what is lost for the mean is gained for the fluctuations, et vice-versa.

Appendix D. Further validation of Eq. (4.7)

This section aims at providing another test of the validity of Eq. (4.7). For this purpose we filled the periodic box with 8 spherical drops (or bubbles) placed at each quarter of the domain. The radius of each droplet (bubble) was adjusted so that the amount of liquid in the domain φ 0 was varied between 10% and 90%. For this test, the computational domain is made of 128 3 grid points. The surface density can then be readily calculated:

In Fig. 10, we plot the values of Σ as given by Eq. (D 1) for φ 0 in the range 10% ď φ 0 ď 90%. Also displayed, the values obtained from either the surface area of the zero level-set or Eq. (4.7). The agreement is again almost perfect. This confirms that Eq. (4.7) is valid.