<figure><img src="image_1.png" /></figure>

HAL
open science

# Isotropic turbulence of variable-density incompressible flows

L. Reynier, B. Di Pierro, F. Alizard

► To cite this version:

L. Reynier, B. Di Pierro, F. Alizard. Isotropic turbulence of variable-density incompressible flows. Physical Review Fluids, 2024, 9 (5), pp.054604. 10.1103/physrevfluids.9.054604. hal-04572051

HAL Id: hal-04572051

https://hal.science/hal-04572051

Submitted on 9 May 2024

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

Public Domain
# Isotropic turbulence of variable-density incompressible flows

L. Reynier,\* B. Di Pierro, and F. Alizard

Université Claude Bernard Lyon 1, LMFA, UMR5509, 69622 Villeurbanne, France

In the present study, the effects of density variations on structures developing in an isotropic incompressible turbulent flow are investigated. Statistical analyses are carried out on datasets obtained from direct numerical simulations of forced turbulence. The discretized variable-density incompressible Navier-Stokes equations are time-advanced with a Fourier-Fourier spectral solver coupled with a semi-implicit second-order in time Runge-Kutta scheme. Turbulence is forced using an extension of the Lundgren method to the variable-density equations including mass diffusion effects. Numerical evidence shows that the introduction of a variable-density field into a turbulent field modifies the coherent structures and the energy spectrum in the inertial range. The analysis of probability density functions of velocity gradients and Lagrangian acceleration suggests an increase in time- and space-intermittency beyond a threshold density ratio associated with the two-fluid mixture. These modifications are not captured by the classical scaling laws of skewness and flatness factors given in the literature for the constant-density flow case. The energy spectra preserve the Kolmogorov slope while exhibiting an energy level alteration within the smallest scales of the inertial range. This region corresponds to the range of modes where the energy levels of the Rayleigh-Taylor instability criterion are the highest, giving some physical arguments that the aforementioned structural modifications may be attributed to Rayleigh-Taylor-like instabilities.

Keywords: Variable-density flows; isotropic turbulence; incompressible Navier-Stokes equations; direct numerical simulation

## I. INTRODUCTION

Flows characterized by large spatial density variations play a crucial role in various industrial and environmental applications, such as releases of pollutants to the environment. It is therefore of fundamental importance to investigate transport properties of gas mixtures into the atmosphere. In this context, the study of isotropic turbulent structures that may emerge from such flows can provide a preliminary understanding of how to more efficiently predict the dispersion of pollutants in the atmosphere.

In this vein, intensive pioneering work carried out by Livescu and Ristorcelli [1], Livescu et al. [2], Rao et al. [3], and Nomura and Elghobashi [4], among a few others, was devoted to exploring the properties of variable-density turbulence and its interactions with other physical phenomena. These authors deeply explore the energy transfer and the mixing process of variable-density turbulence in buoyancy-driven flows, and nonpremixed flames flow configurations. Mohaghar et al. [5] focus on the interactions between variable-density isotropic turbulence and shock waves in a mixing process, and Tian et al. [6] investigate the density effects on the structure and dynamics of post-shock turbulence.

Although the above findings constitute an important step forward, the understanding of the small-scale statistics of variable-density isotropic turbulence is still lacking. Motivated to overcoming this blind spot, the main objective of the present work is to address this problem using direct numerical simulations (DNS) and statistical analyses.

Following the initial research by Kolmogorov [7], which states that small-scale statistics of turbulence are unaffected by forcing and boundary conditions at larger scales, numerous numerical studies devoted to constant density incompressible flows have been conducted based on triply-periodic DNSs [8–12]. Within this framework, computational costs are reduced, making possible the investigation of small-scales turbulence. Recently it also highlights the potential universality of turbulence small-scale statistics. Especially, recent progress in supercomputing facilities has led to a rapid increase in the resolution level of DNSs (e.g. see Ishihara et al. [13]). This achievement has enabled statistical analyses of small-scale isotropic turbulence in constant-density flows. Sreenivasan and Antonia [14] and more recently Dubrulle [15] undertake a comprehensive review of the recent progress made. Notably, Kerr [8], Jiménez et al. [10], and Ishihara et al. [13], conduct studies on the one-point statistics and energy spectra of small-scale turbulence. Their works clearly show the effect of the microscale Reynolds number on temporal and spatial intermittency which are also strongly connected to the statistics of velocity gradients.

As mentioned above, the behavior of small-scale turbulence is mainly investigated through statistical analysis. Since the turbulent fields have random features, statistics are computed using the time average of datasets to reduce

<footnote>\* loic@loicreynier.fr</footnote>
# Isotropic turbulence of variable-density incompressible flows

L. Reynier,\* B. Di Pierro, and F. Alizard

Université Claude Bernard Lyon 1, LMFA, UMR5509, 69622 Villeurbanne, France

In the present study, the effects of density variations on structures developing in an isotropic incompressible turbulent flow are investigated. Statistical analyses are carried out on datasets obtained from direct numerical simulations of forced turbulence. The discretized variable-density incompressible Navier-Stokes equations are time-advanced with a Fourier-Fourier spectral solver coupled with a semi-implicit second-order in time Runge-Kutta scheme. Turbulence is forced using an extension of the Lundgren method to the variable-density equations including mass diffusion effects. Numerical evidence shows that the introduction of a variable-density field into a turbulent field modifies the coherent structures and the energy spectrum in the inertial range. The analysis of probability density functions of velocity gradients and Lagrangian acceleration suggests an increase in time- and space-intermittency beyond a threshold density ratio associated with the two-fluid mixture. These modifications are not captured by the classical scaling laws of skewness and flatness factors given in the literature for the constant-density flow case. The energy spectra preserve the Kolmogorov slope while exhibiting an energy level alteration within the smallest scales of the inertial range. This region corresponds to the range of modes where the energy levels of the Rayleigh-Taylor instability criterion are the highest, giving some physical arguments that the aforementioned structural modifications may be attributed to Rayleigh-Taylor-like instabilities.

Keywords: Variable-density flows; isotropic turbulence; incompressible Navier-Stokes equations; direct numerical simulation

## I. INTRODUCTION

Flows characterized by large spatial density variations play a crucial role in various industrial and environmental applications, such as releases of pollutants to the environment. It is therefore of fundamental importance to investigate transport properties of gas mixtures into the atmosphere. In this context, the study of isotropic turbulent structures that may emerge from such flows can provide a preliminary understanding of how to more efficiently predict the dispersion of pollutants in the atmosphere.

In this vein, intensive pioneering work carried out by Livescu and Ristorcelli [1], Livescu et al. [2], Rao et al. [3], and Nomura and Elghobashi [4], among a few others, was devoted to exploring the properties of variable-density turbulence and its interactions with other physical phenomena. These authors deeply explore the energy transfer and the mixing process of variable-density turbulence in buoyancy-driven flows, and nonpremixed flames flow configurations. Mohaghar et al. [5] focus on the interactions between variable-density isotropic turbulence and shock waves in a mixing process, and Tian et al. [6] investigate the density effects on the structure and dynamics of post-shock turbulence.

Although the above findings constitute an important step forward, the understanding of the small-scale statistics of variable-density isotropic turbulence is still lacking. Motivated to overcoming this blind spot, the main objective of the present work is to address this problem using direct numerical simulations (DNS) and statistical analyses.

Following the initial research by Kolmogorov [7], which states that small-scale statistics of turbulence are unaffected by forcing and boundary conditions at larger scales, numerous numerical studies devoted to constant density incompressible flows have been conducted based on triply-periodic DNSs [8–12]. Within this framework, computational costs are reduced, making possible the investigation of small-scales turbulence. Recently it also highlights the potential universality of turbulence small-scale statistics. Especially, recent progress in supercomputing facilities has led to a rapid increase in the resolution level of DNSs (e.g. see Ishihara et al. [13]). This achievement has enabled statistical analyses of small-scale isotropic turbulence in constant-density flows. Sreenivasan and Antonia [14] and more recently Dubrulle [15] undertake a comprehensive review of the recent progress made. Notably, Kerr [8], Jiménez et al. [10], and Ishihara et al. [13], conduct studies on the one-point statistics and energy spectra of small-scale turbulence. Their works clearly show the effect of the microscale Reynolds number on temporal and spatial intermittency which are also strongly connected to the statistics of velocity gradients.

As mentioned above, the behavior of small-scale turbulence is mainly investigated through statistical analysis. Since the turbulent fields have random features, statistics are computed using the time average of datasets to reduce

<footnote>\* loic@loicreynier.fr</footnote>
statistical variance. It then requires that processes are statistically stationary in time. However, isotropic turbulence is not stationary, and energy must be injected into the flow continuously, which is the so-called forced isotropic turbulence. This technique allows achieving higher Reynolds number and longer statistics than its possible with isotropic decaying turbulence. This has been widely used since the early works of Siggia and Patterson [16], Overholt and Pope [17], Kerr [8] and Eswaran and Pope [18]. The forcing technique developed by the aforementioned authors prevents the decay of turbulence by injecting energy into the large turbulent structures through a narrow region in the spectral space. These forcing methods rely on the cascade notion of Kolmogorov [7] which states that the small scales of turbulence are independent of the behavior of the large-scale motions. More recently, Lundgren [19] introduces a forcing term in the physical space, where the corresponding source term takes the form $f = qu$, which is analogous to the kinetic energy production term that appears in shear flows. One of the advantages of this approach is its simple extension to variable-density — compressible or incompressible — flows [20, 21]. Carroll and Blanquart [20] and Janin et al. [22] propose modifications to the latter method to reduce statistics oscillations and control integral length scale. Despite the benefits, it also adds some extra control parameters not present in our initial dataset. Especially, for inhomogeneous density flows, these improvements would require additional control parameters which should be a function of the Schmidt Sc numbers, density ratio $s$, and initial condition $\rho_0(x)$ for the density field. Hence, the present study will rely upon the Lundgren's method that will be extended to incompressible equations with mass diffusion.

While Lundgren's method enables simulations of forced variable-density isotropic turbulence, the development of specific numerical methods [23, 24] and analyses of mass diffusion effect on turbulence small-scale statistics are still lightly explored topics. From the above discussion, it seems that a comprehensive understanding of the development of isotropic turbulent structures in incompressible flow within inhomogeneous density fields is still lacking.

The purpose of this paper is to extend the studies of small-scale flow structures that emerge in forced isotropic turbulence to the variable-density flow case. In the next section, following a brief presentation of the system of equations and numerical schemes, the forcing method based upon Lundgren's technique is introduced. The following section presents two studies of variable-density flow simulations. The first study offers a statistical analysis of small-scale structures by examining the probability density functions (PDF) and moments of velocity field gradients, as well as their dependence on the microscale Reynolds number. The second study focuses on providing some physical insights from a structural perspective through the use of Fourier spectral analyses and visualization of flow fields to draw a complete picture. Finally, the conclusions are summarized in the last section and prospects are given.

## II. GOVERNING EQUATIONS AND NUMERICAL METHODS

### A. Mathematical model

We consider the isotropic turbulent motion of a viscous fluid in an inhomogeneous medium that takes places in a bounded cubic domain $\Omega$ of volume $V$ — where the Cartesian coordinate system is defined by the $x$, $y$ and $z$ axes — and in time interval $t \in [0, t_f]$. The turbulence fields are assumed to be periodic in each direction of the Cartesian coordinates within the periodic box $\Omega$ of size $2\pi$, so it can be expressed as a Fourier series with both the minimum wavenumber $k_{min}$ and the wavenumber increment being 1. The mathematical model for the variable-density incompressible Navier-Stokes equations (VDINSEs) used in the present contribution is detailed by Frank-Kamenetskii [25], Kazhikhov and Smagulov [26], Antontsev et al. [27], and Guillén-González et al. [28]. Within this framework, mass diffusion is modeled according to Fick's diffusion law and large density variations are considered — where the Boussinesq approximation is no longer verified — while being sufficiently smooth to be accurately projected onto a spectral basis. In regard to the above-mentioned references, we introduce the mean density $\rho(x, t)$ and the mean-volume velocities $u = (u, v, w)^T(x, t)$. The dimensionless equations of motion read

$$ \frac{\partial u}{\partial t} + (u \cdot \nabla)u + \frac{\nabla p}{\rho} = \zeta(\rho, u) + f, \qquad (1a) $$

$$ \frac{\partial \rho}{\partial t} + (u \cdot \nabla)\rho = \frac{1}{Re Sc} \nabla^2 \rho, \qquad (1b) $$

$$ \nabla \cdot u = 0, \qquad (1c) $$

$$ \zeta(\rho, u) = \frac{1}{\rho Re} \nabla^2 u + \frac{1}{\rho Re Sc} ((u \cdot \nabla)\nabla\rho + (\nabla\rho \cdot \nabla)u), \qquad (1d) $$

where $\zeta(\rho, u)$ represents the — viscous and mass — momentum diffusion terms. In the momentum equation (1a), $p$ is a potential function analogous to the pressure and $f$ represents an external volumetric body force. In equations
statistical variance. It then requires that processes are statistically stationary in time. However, isotropic turbulence is not stationary, and energy must be injected into the flow continuously, which is the so-called forced isotropic turbulence. This technique allows achieving higher Reynolds number and longer statistics than its possible with isotropic decaying turbulence. This has been widely used since the early works of Siggia and Patterson [16], Overholt and Pope [17], Kerr [8] and Eswaran and Pope [18]. The forcing technique developed by the aforementioned authors prevents the decay of turbulence by injecting energy into the large turbulent structures through a narrow region in the spectral space. These forcing methods rely on the cascade notion of Kolmogorov [7] which states that the small scales of turbulence are independent of the behavior of the large-scale motions. More recently, Lundgren [19] introduces a forcing term in the physical space, where the corresponding source term takes the form $f = qu$, which is analogous to the kinetic energy production term that appears in shear flows. One of the advantages of this approach is its simple extension to variable-density — compressible or incompressible — flows [20, 21]. Carroll and Blanquart [20] and Janin et al. [22] propose modifications to the latter method to reduce statistics oscillations and control integral length scale. Despite the benefits, it also adds some extra control parameters not present in our initial dataset. Especially, for inhomogeneous density flows, these improvements would require additional control parameters which should be a function of the Schmidt Sc numbers, density ratio $s$, and initial condition $\rho_0(x)$ for the density field. Hence, the present study will rely upon the Lundgren's method that will be extended to incompressible equations with mass diffusion.

While Lundgren's method enables simulations of forced variable-density isotropic turbulence, the development of specific numerical methods [23, 24] and analyses of mass diffusion effect on turbulence small-scale statistics are still lightly explored topics. From the above discussion, it seems that a comprehensive understanding of the development of isotropic turbulent structures in incompressible flow within inhomogeneous density fields is still lacking.

The purpose of this paper is to extend the studies of small-scale flow structures that emerge in forced isotropic turbulence to the variable-density flow case. In the next section, following a brief presentation of the system of equations and numerical schemes, the forcing method based upon Lundgren's technique is introduced. The following section presents two studies of variable-density flow simulations. The first study offers a statistical analysis of small-scale structures by examining the probability density functions (PDF) and moments of velocity field gradients, as well as their dependence on the microscale Reynolds number. The second study focuses on providing some physical insights from a structural perspective through the use of Fourier spectral analyses and visualization of flow fields to draw a complete picture. Finally, the conclusions are summarized in the last section and prospects are given.

## II. GOVERNING EQUATIONS AND NUMERICAL METHODS

### A. Mathematical model

We consider the isotropic turbulent motion of a viscous fluid in an inhomogeneous medium that takes places in a bounded cubic domain $\Omega$ of volume $V$ — where the Cartesian coordinate system is defined by the $x$, $y$ and $z$ axes — and in time interval $t \in [0, t_f]$. The turbulence fields are assumed to be periodic in each direction of the Cartesian coordinates within the periodic box $\Omega$ of size $2\pi$, so it can be expressed as a Fourier series with both the minimum wavenumber $k_{min}$ and the wavenumber increment being 1. The mathematical model for the variable-density incompressible Navier-Stokes equations (VDINSEs) used in the present contribution is detailed by Frank-Kamenetskii [25], Kazhikhov and Smagulov [26], Antontsev et al. [27], and Guillén-González et al. [28]. Within this framework, mass diffusion is modeled according to Fick's diffusion law and large density variations are considered — where the Boussinesq approximation is no longer verified — while being sufficiently smooth to be accurately projected onto a spectral basis. In regard to the above-mentioned references, we introduce the mean density $\rho(x, t)$ and the mean-volume velocities $u = (u, v, w)^T(x, t)$. The dimensionless equations of motion read

$$ \frac{\partial u}{\partial t} + (u \cdot \nabla)u + \frac{\nabla p}{\rho} = \zeta(\rho, u) + f, \qquad (1a) $$

$$ \frac{\partial \rho}{\partial t} + (u \cdot \nabla)\rho = \frac{1}{Re Sc} \nabla^2 \rho, \qquad (1b) $$

$$ \nabla \cdot u = 0, \qquad (1c) $$

$$ \zeta(\rho, u) = \frac{1}{\rho Re} \nabla^2 u + \frac{1}{\rho Re Sc} ((u \cdot \nabla)\nabla\rho + (\nabla\rho \cdot \nabla)u), \qquad (1d) $$

where $\zeta(\rho, u)$ represents the — viscous and mass — momentum diffusion terms. In the momentum equation (1a), $p$ is a potential function analogous to the pressure and $f$ represents an external volumetric body force. In equations
(1b) and (1d), Re and Sc are the Reynolds and Schmidt numbers respectively, defined as

$$Re = \frac{\rho_c v_c L_c}{\eta}, \quad Sc = \frac{\eta}{\rho_c D}, \qquad (2)$$

with η and D the dynamic viscosity and mass diffusivity of the fluid, respectively, and ρ_c, v_c, L_c being characteristic density, velocity and length scales.

It is essential to note that the incompressible formulation arises from the selection of the mean-volume velocity **u** for representing the velocity field of the variable-density flow. Conventionally, the Navier–Stokes equations are written in terms of the barycentric mean-mass velocity **u**, which does not lead to a divergence-free velocity field. Equations (1a) and (1b) are derived by initially formulating the Navier–Stokes equations on the mass-averaged velocity and subsequently substituting the latter with the mean-volume velocity using Fick’s diffusion law

$$\bar{u} = u - \frac{D}{\rho} \nabla \rho, \qquad (3)$$

which relates the velocity $\bar{u}$ with nonzero divergence to the divergence-free velocity **u**. Additional details are provided by Guillén-González et al. [28].

The initial conditions of the system (1) are

$$u(x, 0) = u_0(x), \quad p(x, 0) = p_0(x), \quad \rho(x, 0) = \rho_0(x), \quad \forall x \in \Omega \qquad (4)$$

where $u_0$ and $p_0$ are the initial turbulent velocity and pressure fields obtained from a preliminary constant-density forced isotropic turbulence simulation and $\rho_0$ is the smooth initial density profile. We then introduce an additional control parameter which characterizes the density ratio:

$$s = \frac{\max \rho_0(x)}{\min \rho_0(x)}. \qquad (5)$$

## B. Numerical methods

VDINSE (1) are solved using a Fourier–Fourier spectral numerical scheme. Turbulence fields are expressed as Fourier series where aliasing errors are removed by phase-shifting that keeps all the Fourier modes satisfying $|k| < k_{max} = N/3$, where k is the wavenumber and N the number of grid points along each spatial Cartesian direction. The time-stepping procedure is a second-order Runge–Kutta method for the nonlinear terms with semi-implicit Crank–Nicholson treatment for the momentum diffusion terms (RK2CN). This is based on the work of Reynier et al. [24], which is an extension of Peyret [29], Bell and Marcus [30], Tadjeran [31], Di Pierro and Abid [32], and which is briefly outlined below:

$$\frac{\mathbf{u}^{n+1} - \mathbf{u}^n}{\delta t} = -[(\mathbf{u} \cdot \nabla)\mathbf{u}]^{n+1/2} - \frac{\nabla p^{n+1/2}}{\rho^{n+1/2}} + \frac{1}{2}(\zeta(\rho^{n+1/2}, \mathbf{u}^{n+1}) + \zeta(\rho^{n+1/2}, \mathbf{u}^n)) + \mathbf{f}^n, \qquad (6a)$$

$$\frac{\rho^{n+1} - \rho^n}{\delta t} = -[(\mathbf{u} \cdot \nabla)\rho]^{n+1/2} + \frac{1}{2} \frac{1}{Re Sc} (\nabla^2 \rho^{n+1} + \nabla^2 \rho^n), \qquad (6b)$$

$$\nabla \cdot \mathbf{u}^{n+1} = 0. \qquad (6c)$$

The n superscript symbolizes the solution at time $t^n = n\delta t$ and the n + 1/2 superscript refers to the second-order sub-step solution at time $t^{n+1/2} = (n + \frac{1}{2})\delta t$, where $\delta t$ is the integration time-step. Advection terms $[(\mathbf{u} \cdot \nabla)\mathbf{u}]^{n+1/2}$ and $[(\mathbf{u} \cdot \nabla)\rho]^{n+1/2}$ are estimated through a Runge–Kutta 2 (RK2) scheme. The momentum diffusion term $\zeta(\rho, \mathbf{u})$ is computed with the sub-step density $\rho^{n+1/2} = \frac{1}{2}(\rho^{n+1} + \rho^n)$ following Tadjeran [31] for stability considerations. In order to preserve the second-order accuracy in time and to ensure incompressibility, the pressure — and the update pressure $p^{n+1/2}$ — is computed with the projection technique proposed by Reynier et al. [24], by solving

$$\nabla \cdot \left( \frac{1}{\rho^{n+1/2}} \nabla p^{n+1/2} \right) = f(\mathbf{u}) \qquad (7)$$

through the inversion of the pressure operator. This is performed using a Generalized Minimal Residual (GMRES) algorithm preconditioned with the inverse of the discrete Laplacian operator. The semi-implicit treatment of the
(1b) and (1d), Re and Sc are the Reynolds and Schmidt numbers respectively, defined as

$$Re = \frac{\rho_c v_c L_c}{\eta}, \quad Sc = \frac{\eta}{\rho_c D}, \qquad (2)$$

with $\eta$ and $D$ the dynamic viscosity and mass diffusivity of the fluid, respectively, and $\rho_c$, $v_c$, $L_c$ being characteristic density, velocity and length scales.

It is essential to note that the incompressible formulation arises from the selection of the mean-volume velocity $\mathbf{u}$ for representing the velocity field of the variable-density flow. Conventionally, the Navier–Stokes equations are written in terms of the barycentric mean-mass velocity $\bar{\mathbf{u}}$, which does not lead to a divergence-free velocity field. Equations (1a) and (1b) are derived by initially formulating the Navier–Stokes equations on the mass-averaged velocity and subsequently substituting the latter with the mean-volume velocity using Fick’s diffusion law

$$\bar{\mathbf{u}} = \mathbf{u} - \frac{D}{\rho} \nabla \rho, \qquad (3)$$

which relates the velocity $\bar{\mathbf{u}}$ with nonzero divergence to the divergence-free velocity $\mathbf{u}$. Additional details are provided by Guillén-González et al. [28].

The initial conditions of the system (1) are

$$\mathbf{u}(\mathbf{x}, 0) = \mathbf{u}_0(\mathbf{x}), \quad p(\mathbf{x}, 0) = p_0(\mathbf{x}), \quad \rho(\mathbf{x}, 0) = \rho_0(\mathbf{x}), \quad \forall \mathbf{x} \in \Omega \qquad (4)$$

where $\mathbf{u}_0$ and $p_0$ are the initial turbulent velocity and pressure fields obtained from a preliminary constant-density forced isotropic turbulence simulation and $\rho_0$ is the smooth initial density profile. We then introduce an additional control parameter which characterizes the density ratio:

$$s = \frac{\max \rho_0(\mathbf{x})}{\min \rho_0(\mathbf{x})}. \qquad (5)$$

## B. Numerical methods

VDINSE (1) are solved using a Fourier–Fourier spectral numerical scheme. Turbulence fields are expressed as Fourier series where aliasing errors are removed by phase-shifting that keeps all the Fourier modes satisfying $|k| < k_{max} = N/3$, where $k$ is the wavenumber and $N$ the number of grid points along each spatial Cartesian direction. The time-stepping procedure is a second-order Runge–Kutta method for the nonlinear terms with semi-implicit Crank–Nicholson treatment for the momentum diffusion terms (RK2CN). This is based on the work of Reynier et al. [24], which is an extension of Peyret [29], Bell and Marcus [30], Tadjeran [31], Di Pierro and Abid [32], and which is briefly outlined below:

$$\frac{\mathbf{u}^{n+1} - \mathbf{u}^n}{\delta t} = -[(\mathbf{u} \cdot \nabla)\mathbf{u}]^{n+1/2} - \frac{\nabla p^{n+1/2}}{\rho^{n+1/2}} + \frac{1}{2}(\zeta(\rho^{n+1/2}, \mathbf{u}^{n+1}) + \zeta(\rho^{n+1/2}, \mathbf{u}^n)) + \mathbf{f}^n, \qquad (6a)$$

$$\frac{\rho^{n+1} - \rho^n}{\delta t} = -[(\mathbf{u} \cdot \nabla)\rho]^{n+1/2} + \frac{1}{2} \frac{1}{Re Sc} (\nabla^2 \rho^{n+1} + \nabla^2 \rho^n), \qquad (6b)$$

$$\nabla \cdot \mathbf{u}^{n+1} = 0. \qquad (6c)$$

The $n$ superscript symbolizes the solution at time $t^n = n\delta t$ and the $n + 1/2$ superscript refers to the second-order sub-step solution at time $t^{n+1/2} = (n + 1/2)\delta t$, where $\delta t$ is the integration time-step. Advection terms $[(\mathbf{u} \cdot \nabla)\mathbf{u}]^{n+1/2}$ and $[(\mathbf{u} \cdot \nabla)\rho]^{n+1/2}$ are estimated through a Runge–Kutta 2 (RK2) scheme. The momentum diffusion term $\zeta(\rho, \mathbf{u})$ is computed with the sub-step density $\rho^{n+1/2} = \frac{1}{2}(\rho^{n+1} + \rho^n)$ following Tadjeran [31] for stability considerations. In order to preserve the second-order accuracy in time and to ensure incompressibility, the pressure — and the update pressure $p^{n+1/2}$ — is computed with the projection technique proposed by Reynier et al. [24], by solving

$$\nabla \cdot \left( \frac{1}{\rho^{n+1/2}} \nabla p^{n+1/2} \right) = f(\mathbf{u}) \qquad (7)$$

through the inversion of the pressure operator. This is performed using a Generalized Minimal Residual (GMRES) algorithm preconditioned with the inverse of the discrete Laplacian operator. The semi-implicit treatment of the
equations of motion also requires the inversion of the momentum diffusion operator. These computations are performed with a GMRES method preconditioned with the inverse operator associated with constant-density flow as suggested by Reynier et al. [24]. These preconditioning techniques greatly improve the convergence properties of the iterative solver by reducing the condition number of the aforementioned operators.

Finally, to prevent the emergence of steep structures that are challenging to capture with the spectral scheme, the density field is filtered after each time step when deemed necessary with a spectral Gaussian filter

$$ \rho_f^n = \mathcal{F}^{-1}(\hat{\rho}^n(\mathbf{k}) \exp(-\sigma(\mathbf{k}/k_f)^p)) \quad (8) $$

where $\rho_f^n$ is the filtered density field after the nth time-advancement. In eq. (8) and following, the spectral representation of $\varphi^n$ is noted $\hat{\varphi}^n$ and the Fourier transform of a given function is symbolized $\mathcal{F}(\cdot)$. The filter parameters are set to $k_f = 0.55k_{max}$, $\sigma = 18$ and $p = 16$, hereafter. Specifically, the filter is applied at the beginning of the simulations, spanning 100 iterations. This is done to prevent the emergence of problematic steep structures following the injection of the density field, as elaborated upon in section III A.

### C. Turbulence forcing

Stationary isotropic turbulence is studied numerically by injecting energy through the forcing term $f$ in the momentum equation (1a). We recall that for the constant-density case, the rate of change of the turbulent kinetic energy per unit of mass reads

$$ \int_{\Omega} \frac{\partial}{\partial t} (\frac{1}{2} \mathbf{u}^2) \, dV = \int_{\Omega} \frac{1}{Re} \mathbf{u} \cdot \nabla^2 \mathbf{u} \, dV + \int_{\Omega} \mathbf{f} \cdot \mathbf{u} \, dV. \quad (9) $$

The forcing term is adjusted to achieve a statistically stationary state which guarantees that the integrated kinetic energy is conserved. In that regard, we should verify

$$ \int_{\Omega} \mathbf{f} \cdot \mathbf{u} \, dV = - \int_{\Omega} \frac{1}{Re} \mathbf{u} \cdot \nabla^2 \mathbf{u} \, dV. \quad (10) $$

Lundgren [19] proposes a linear method in the physical space where $f = qu$ injects energy proportionally to the velocity field fluctuations. This approach injects energy over all scales of the flow unlike the spectral forcing methods — such as the method of Kerr [8] — which excite only a narrow region of the spectral space. As the largest scales are subjects to the largest fluctuations, these scales are the most impacted ones by the source term and the smallest scales are essentially unaffected [19, 20]. This method has the advantage of being conveniently extendable to variable-density simulations. The forcing coefficient $q$ is computed from eq. (9) to verify the stationarity condition of eq. (10) and reads

$$ q(\mathbf{u}) = - \frac{1}{Re} \frac{\int_{\Omega} \mathbf{u} \cdot \nabla^2 \mathbf{u} \, dV}{\int_{\Omega} \mathbf{u} \cdot \mathbf{u} \, dV}. \quad (11) $$

The forcing term being proportional to the velocity field, the latter must be initialized. To ensure that small structures are not affected by the initial condition, we adopt an initial condition $u_0$ similar to the forcing term of Kerr's method which forces the Fourier modes associated to the smallest wavenumbers to a constant amplitude:

$$ u_0 = \mathcal{F}^{-1}(\hat{f}), \quad \hat{f}(\mathbf{k}) = f_0, \quad \|\mathbf{k}\|_2 = 1. \quad (12) $$

with $f_0 = \frac{1}{4}$ so that the initial velocity amplitude is unity. This initial velocity condition is finally projected into the divergence-free space.

For variable-density flows, the integrated turbulent kinetic energy per unit of volume $\frac{1}{2}\rho u^2$ is now considered. The latter is driven by

$$ \begin{aligned} \frac{d}{dt} (\frac{1}{2} \rho u^2) &= \rho \frac{d}{dt} (\frac{1}{2} u^2) + \frac{1}{2} u^2 \frac{d\rho}{dt} \\ &= -u \cdot \nabla p + \rho u \cdot \zeta(\rho, u) + \rho u \cdot \mathbf{f} + \frac{1}{2 Re Sc} u^2 \nabla^2 \rho. \end{aligned} \quad (13) $$

Following the same line of thought as in the constant-density case, the forcing term should verify

$$ \int_{\Omega} \rho \mathbf{f} \cdot \mathbf{u} \, dV = - \int_{\Omega} \rho \mathbf{u} \cdot \zeta(\rho, \mathbf{u}) \, dV - \int_{\Omega} \frac{1}{2 Re Sc} \mathbf{u}^2 \nabla^2 \rho \, dV. \quad (14) $$
equations of motion also requires the inversion of the momentum diffusion operator. These computations are performed with a GMRES method preconditioned with the inverse operator associated with constant-density flow as suggested by Reynier et al. [24]. These preconditioning techniques greatly improve the convergence properties of the iterative solver by reducing the condition number of the aforementioned operators.

Finally, to prevent the emergence of steep structures that are challenging to capture with the spectral scheme, the density field is filtered after each time step when deemed necessary with a spectral Gaussian filter

$$ \rho_f^n = \mathcal{F}^{-1}(\hat{\rho}^n(\mathbf{k}) \exp(-\sigma(\mathbf{k}/k_f)^p)) \quad (8) $$

where $\rho_f^n$ is the filtered density field after the nth time-advancement. In eq. (8) and following, the spectral representation of $\varphi^n$ is noted $\hat{\varphi}^n$ and the Fourier transform of a given function is symbolized $\mathcal{F}(\cdot)$. The filter parameters are set to $k_f = 0.55k_{max}$, $\sigma = 18$ and $p = 16$, hereafter. Specifically, the filter is applied at the beginning of the simulations, spanning 100 iterations. This is done to prevent the emergence of problematic steep structures following the injection of the density field, as elaborated upon in section III A.

### C. Turbulence forcing

Stationary isotropic turbulence is studied numerically by injecting energy through the forcing term $f$ in the momentum equation (1a). We recall that for the constant-density case, the rate of change of the turbulent kinetic energy per unit of mass reads

$$ \int_{\Omega} \frac{\partial}{\partial t} (\frac{1}{2} \mathbf{u}^2) \, dV = \int_{\Omega} \frac{1}{Re} \mathbf{u} \cdot \nabla^2 \mathbf{u} \, dV + \int_{\Omega} \mathbf{f} \cdot \mathbf{u} \, dV. \quad (9) $$

The forcing term is adjusted to achieve a statistically stationary state which guarantees that the integrated kinetic energy is conserved. In that regard, we should verify

$$ \int_{\Omega} \mathbf{f} \cdot \mathbf{u} \, dV = - \int_{\Omega} \frac{1}{Re} \mathbf{u} \cdot \nabla^2 \mathbf{u} \, dV. \quad (10) $$

Lundgren [19] proposes a linear method in the physical space where $f = qu$ injects energy proportionally to the velocity field fluctuations. This approach injects energy over all scales of the flow unlike the spectral forcing methods — such as the method of Kerr [8] — which excite only a narrow region of the spectral space. As the largest scales are subjects to the largest fluctuations, these scales are the most impacted ones by the source term and the smallest scales are essentially unaffected [19, 20]. This method has the advantage of being conveniently extendable to variable-density simulations. The forcing coefficient $q$ is computed from eq. (9) to verify the stationarity condition of eq. (10) and reads

$$ q(\mathbf{u}) = - \frac{1}{Re} \frac{\int_{\Omega} \mathbf{u} \cdot \nabla^2 \mathbf{u} \, dV}{\int_{\Omega} \mathbf{u} \cdot \mathbf{u} \, dV}. \quad (11) $$

The forcing term being proportional to the velocity field, the latter must be initialized. To ensure that small structures are not affected by the initial condition, we adopt an initial condition $u_0$ similar to the forcing term of Kerr's method which forces the Fourier modes associated to the smallest wavenumbers to a constant amplitude:

$$ u_0 = \mathcal{F}^{-1}(\hat{f}), \quad \hat{f}(\mathbf{k}) = f_0, \quad \|\mathbf{k}\|_2 = 1. \quad (12) $$

with $f_0 = \frac{1}{4}$ so that the initial velocity amplitude is unity. This initial velocity condition is finally projected into the divergence-free space.

For variable-density flows, the integrated turbulent kinetic energy per unit of volume $\frac{1}{2}\rho u^2$ is now considered. The latter is driven by

$$ \begin{aligned} \frac{d}{dt} (\frac{1}{2} \rho u^2) &= \rho \frac{d}{dt} (\frac{1}{2} u^2) + \frac{1}{2} u^2 \frac{d\rho}{dt} \\ &= -u \cdot \nabla p + \rho u \cdot \zeta(\rho, u) + \rho u \cdot \mathbf{f} + \frac{1}{2 Re Sc} u^2 \nabla^2 \rho. \end{aligned} \quad (13) $$

Following the same line of thought as in the constant-density case, the forcing term should verify

$$ \int_{\Omega} \rho \mathbf{f} \cdot \mathbf{u} \, dV = - \int_{\Omega} \rho \mathbf{u} \cdot \zeta(\rho, \mathbf{u}) \, dV - \int_{\Omega} \frac{1}{2 Re Sc} \mathbf{u}^2 \nabla^2 \rho \, dV. \quad (14) $$
<table><caption>TABLE I: DNS parameters of variable-density forced turbulence simulations, and turbulence characteristics at the final time $t = t_f$ and averaged $⟨·⟩_Γ$ over the study time interval $Γ$. $N$ represents both the number of Fourier modes (prior to aliasing) and the number of mesh points in each direction of physical space. To enhance readability, the notation $⟨·⟩_Γ$ is omitted in the text. Final velocity and pressure fields from the constant-density simulation 256-1 are used as initial conditions along with the sphere density profile from eq. (18) for the variable-density simulations.</caption><thead><tr><th>Run</th><th>N</th><th>Re</th><th>s</th><th>d</th><th>Sc</th><th>t<sub>f</sub></th><th>Γ</th><th>⟨Re<sub>λ</sub>⟩<sub>Γ</sub></th><th>⟨T⟩<sub>Γ</sub></th><th>⟨L⟩<sub>Γ</sub></th><th>⟨λ⟩<sub>Γ</sub></th><th>10<sup>3</sup>⟨η⟩<sub>Γ</sub></th><th>k<sub>max</sub>⟨η⟩<sub>Γ</sub></th><th>⟨λ⟩<sub>Γ</sub>/⟨η⟩<sub>Γ</sub></th><th>−⟨S⟩<sub>Γ</sub></th><th>⟨F⟩<sub>Γ</sub></th></tr></thead><tbody><tr><td>256-0.25</td><td>256</td><td>10<sup>3</sup></td><td>0.25</td><td>10</td><td>1</td><td>6.0</td><td>[2, 6]</td><td>249</td><td>1.43</td><td>0.81</td><td>0.20</td><td>6.41</td><td>1.1</td><td>31.0</td><td>0.51</td><td>4.92</td></tr><tr><td>256-0.5</td><td>256</td><td>10<sup>3</sup></td><td>0.5</td><td>10</td><td>1</td><td>10.0</td><td>[2, 6]</td><td>246</td><td>1.43</td><td>0.81</td><td>0.20</td><td>6.41</td><td>1.1</td><td>31.0</td><td>0.51</td><td>4.92</td></tr><tr><td>256-1</td><td>256</td><td>10<sup>3</sup></td><td>1</td><td>-</td><td>1</td><td>70.0</td><td>[40, 70]</td><td>214</td><td>1.44</td><td>0.78</td><td>0.17</td><td>5.85</td><td>1.0</td><td>28.8</td><td>0.48</td><td>4.97</td></tr><tr><td>256-2</td><td>256</td><td>10<sup>3</sup></td><td>2</td><td>10</td><td>1</td><td>6.0</td><td>[2, 6]</td><td>172</td><td>1.42</td><td>0.81</td><td>0.14</td><td>5.39</td><td>0.92</td><td>25.9</td><td>0.43</td><td>4.82</td></tr><tr><td>256-4</td><td>256</td><td>10<sup>3</sup></td><td>4</td><td>10</td><td>1</td><td>10.0</td><td>[2, 6]</td><td>143</td><td>1.39</td><td>0.80</td><td>0.12</td><td>4.86</td><td>0.83</td><td>23.7</td><td>0.30</td><td>4.37</td></tr><tr><td>256-6</td><td>256</td><td>10<sup>3</sup></td><td>6</td><td>10</td><td>1</td><td>6.0</td><td>[2, 6]</td><td>120</td><td>1.39</td><td>0.80</td><td>0.09</td><td>4.19</td><td>0.72</td><td>20.0</td><td>0.18</td><td>4.03</td></tr><tr><td>256-8</td><td>256</td><td>10<sup>3</sup></td><td>8</td><td>10</td><td>1</td><td>10.0</td><td>[2, 6]</td><td>108</td><td>1.36</td><td>0.80</td><td>0.08</td><td>4.19</td><td>0.72</td><td>20.5</td><td>0.18</td><td>4.03</td></tr></tbody></table>

<table><caption>TABLE II: DNS parameters and turbulence characteristics at the final time $t = t_f$ of the constant-density HIT simulation employed to generate the turbulent initial conditions for the variable-density HIT simulations.</caption><thead><tr><th>Run</th><th>N</th><th>Re</th><th>t<sub>f</sub></th><th>Re<sub>λ</sub></th><th>T</th><th>L</th><th>λ</th><th>10<sup>3</sup>η</th><th>k<sub>max</sub>η</th><th>λ/η</th><th>-S</th><th>F</th></tr></thead><tbody><tr><td>256-1</td><td>256</td><td>10<sup>3</sup></td><td>40</td><td>210</td><td>1.56</td><td>0.89</td><td>0.19</td><td>6.54</td><td>1.1</td><td>28.6</td><td>0.48</td><td>5.36</td></tr></tbody></table>

Within this framework, the forcing coefficient $q$ reads

$$q(\rho, \mathbf{u}) = - \frac{\int_{\Omega} \mathbf{u} \cdot (\rho \zeta(\rho, \mathbf{u}) + (2 \mathrm{Re} \, \mathrm{Sc})^{-1} \mathbf{u} \nabla^2 \rho) \, dV}{\int_{\Omega} \rho \mathbf{u} \cdot \mathbf{u} \, dV}, \quad (15)$$

which is now a function of both $\rho$ and $\mathbf{u}$.

### III. RESULTS

#### A. Simulation methodology

Table I summarizes simulation parameters and some key turbulent parameters values characterizing the simulation runs. The integral length scale $L$ and the Taylor microscale $\lambda$ are computed from the one-dimensional energy spectrum $E(k)$ such as

$$L = \frac{3\pi}{4} \int_{0}^{k_{\max}} k^{-1} E(k) \, dk \bigg/ \int_{0}^{k_{\max}} E(k) \, dk, \quad (16)$$

and

$$\lambda = \left( 5 \int_{0}^{k_{\max}} E(k) \, dk \bigg/ \int_{0}^{k_{\max}} k^2 E(k) \, dk \right)^{\frac{1}{2}}. \quad (17)$$

The microscale Reynolds number $Re_λ$ is defined as $Re_λ = (λ/L) Re$. The eddy turnover time $T$ is given by $T = L/u'$, where $\frac{3}{2}u'^2$ is the total turbulent kinetic energy per unit of mass.

A preliminary simulation is carried out with constant-density using Lundgren's method described by eqs. (10) to (12). Turbulent parameters of this specific run can be found in table I under the label 256-1. The (turbulent) resulting velocity field is used as the initial condition for variable-density simulations. Table II presents the turbulent characteristics of this initial condition and fig. 1 displays its energy spectra.

In this work, the density field is initialized with a smooth profile $ρ_0(x)$ in which a sphere of volume $\frac{1}{2}V$ is located at the center of the computational box

$$\rho_0(r) = 1 + \frac{s-1}{2} \left( 1 - \tanh \frac{r-r_0}{d} \right) \quad (18)$$
<table><caption>TABLE I: DNS parameters of variable-density forced turbulence simulations, and turbulence characteristics at the final time $t = t_f$ and averaged $⟨·⟩_Γ$ over the study time interval $Γ$. $N$ represents both the number of Fourier modes (prior to aliasing) and the number of mesh points in each direction of physical space. To enhance readability, the notation $⟨·⟩_Γ$ is omitted in the text. Final velocity and pressure fields from the constant-density simulation 256-1 are used as initial conditions along with the sphere density profile from eq. (18) for the variable-density simulations.</caption><thead><tr><th>Run</th><th>N</th><th>Re</th><th>s</th><th>d</th><th>Sc</th><th>t<sub>f</sub></th><th>Γ</th><th>⟨Re<sub>λ</sub>⟩<sub>Γ</sub></th><th>⟨T⟩<sub>Γ</sub></th><th>⟨L⟩<sub>Γ</sub></th><th>⟨λ⟩<sub>Γ</sub></th><th>10<sup>3</sup>⟨η⟩<sub>Γ</sub></th><th>k<sub>max</sub>⟨η⟩<sub>Γ</sub></th><th>⟨λ⟩<sub>Γ</sub>/⟨η⟩<sub>Γ</sub></th><th>−⟨S⟩<sub>Γ</sub></th><th>⟨F⟩<sub>Γ</sub></th></tr></thead><tbody><tr><td>256-0.25</td><td>256</td><td>10<sup>3</sup></td><td>0.25</td><td>10</td><td>1</td><td>6.0</td><td>[2, 6]</td><td>249</td><td>1.43</td><td>0.81</td><td>0.20</td><td>6.41</td><td>1.1</td><td>31.0</td><td>0.51</td><td>4.92</td></tr><tr><td>256-0.5</td><td>256</td><td>10<sup>3</sup></td><td>0.5</td><td>10</td><td>1</td><td>10.0</td><td>[2, 6]</td><td>246</td><td>1.43</td><td>0.81</td><td>0.20</td><td>6.41</td><td>1.1</td><td>31.0</td><td>0.51</td><td>4.92</td></tr><tr><td>256-1</td><td>256</td><td>10<sup>3</sup></td><td>1</td><td>–</td><td>1</td><td>70.0</td><td>[40, 70]</td><td>214</td><td>1.44</td><td>0.78</td><td>0.17</td><td>5.85</td><td>1.0</td><td>28.8</td><td>0.48</td><td>4.97</td></tr><tr><td>256-2</td><td>256</td><td>10<sup>3</sup></td><td>2</td><td>10</td><td>1</td><td>6.0</td><td>[2, 6]</td><td>172</td><td>1.42</td><td>0.81</td><td>0.14</td><td>5.39</td><td>0.92</td><td>25.9</td><td>0.43</td><td>4.82</td></tr><tr><td>256-4</td><td>256</td><td>10<sup>3</sup></td><td>4</td><td>10</td><td>1</td><td>10.0</td><td>[2, 6]</td><td>143</td><td>1.39</td><td>0.80</td><td>0.12</td><td>4.86</td><td>0.83</td><td>23.7</td><td>0.30</td><td>4.37</td></tr><tr><td>256-6</td><td>256</td><td>10<sup>3</sup></td><td>6</td><td>10</td><td>1</td><td>6.0</td><td>[2, 6]</td><td>120</td><td>1.39</td><td>0.80</td><td>0.09</td><td>4.19</td><td>0.72</td><td>20.0</td><td>0.18</td><td>4.03</td></tr><tr><td>256-8</td><td>256</td><td>10<sup>3</sup></td><td>8</td><td>10</td><td>1</td><td>10.0</td><td>[2, 6]</td><td>108</td><td>1.36</td><td>0.80</td><td>0.08</td><td>4.19</td><td>0.72</td><td>20.5</td><td>0.18</td><td>4.03</td></tr></tbody></table>

<table><caption>TABLE II: DNS parameters and turbulence characteristics at the final time $t = t_f$ of the constant-density HIT simulation employed to generate the turbulent initial conditions for the variable-density HIT simulations.</caption><thead><tr><th>Run</th><th>N</th><th>Re</th><th>t<sub>f</sub></th><th>Re<sub>λ</sub></th><th>T</th><th>L</th><th>λ</th><th>10<sup>3</sup>η</th><th>k<sub>max</sub>η</th><th>λ/η</th><th>-S</th><th>F</th></tr></thead><tbody><tr><td>256-1</td><td>256</td><td>10<sup>3</sup></td><td>40</td><td>210</td><td>1.56</td><td>0.89</td><td>0.19</td><td>6.54</td><td>1.1</td><td>28.6</td><td>0.48</td><td>5.36</td></tr></tbody></table>

Within this framework, the forcing coefficient $q$ reads

$$q(\rho, \mathbf{u}) = - \frac{\int_{\Omega} \mathbf{u} \cdot (\rho \zeta(\rho, \mathbf{u}) + (2 \mathrm{Re} \, \mathrm{Sc})^{-1} \mathbf{u} \nabla^2 \rho) \, dV}{\int_{\Omega} \rho \mathbf{u} \cdot \mathbf{u} \, dV}, \quad (15)$$

which is now a function of both $\rho$ and $\mathbf{u}$.

### III. RESULTS

#### A. Simulation methodology

Table I summarizes simulation parameters and some key turbulent parameters values characterizing the simulation runs. The integral length scale $L$ and the Taylor microscale $\lambda$ are computed from the one-dimensional energy spectrum $E(k)$ such as

$$L = \frac{3\pi}{4} \int_{0}^{k_{\max}} k^{-1} E(k) \, dk \bigg/ \int_{0}^{k_{\max}} E(k) \, dk, \quad (16)$$

and

$$\lambda = \left( 5 \int_{0}^{k_{\max}} E(k) \, dk \bigg/ \int_{0}^{k_{\max}} k^2 E(k) \, dk \right)^{\frac{1}{2}}. \quad (17)$$

The microscale Reynolds number $Re_λ$ is defined as $Re_λ = (λ/L) Re$. The eddy turnover time $T$ is given by $T = L/u'$, where $\frac{3}{2}u'^2$ is the total turbulent kinetic energy per unit of mass.

A preliminary simulation is carried out with constant-density using Lundgren's method described by eqs. (10) to (12). Turbulent parameters of this specific run can be found in table I under the label 256-1. The (turbulent) resulting velocity field is used as the initial condition for variable-density simulations. Table II presents the turbulent characteristics of this initial condition and fig. 1 displays its energy spectra.

In this work, the density field is initialized with a smooth profile $ρ_0(x)$ in which a sphere of volume $\frac{1}{2}V$ is located at the center of the computational box

$$\rho_0(r) = 1 + \frac{s-1}{2} \left( 1 - \tanh \frac{r - r_0}{d} \right) \quad (18)$$
<figure><img src="image_1.png" /><figcaption>FIG. 1: Energy spectra at final time $t = t_f$ of constant-density HIT run 256-1 used as initial condition for the variable-density HIT simulations. The vertical line delimits the mode associated with the aliasing filter. Simulation parameters for the constant-density HIT DNS are listed in table II page 5.</figcaption></figure>

where $r = \|x\|_2$ and $r_0 = (3\pi^2)^{1/3}$. The parameter $d$ is the length of the transition region between the interior of the sphere and the surrounding fluid, and is fixed to 10 mesh cells. One may recall that the parameter $s$ represents the density ratio between the fluid filling the sphere and the fluid surrounding it.

In simpler terms, simulations are initialized by the injection of an inhomogeneous density field into a turbulent flow. The abrupt injection of the density field is mitigated by applying the filter described in eq. (8) during the initial 100 iterations of the simulation. This filtering process spans a time interval equivalent to one-tenth of an eddy turnover time.

One can specify that using a turbulent field as initial condition is necessary since injecting the variable-density field before forced turbulence is fully developed will result in complete mixing of the mass field before reaching a statistically stationary turbulent state.

In addition, it is also essential to select a specific time window for the study that encompasses mass diffusion effects while ensuring that the turbulent fields can reach an almost statistically stationary state.

For illustration purposes, the run 256-4 is now considered. Figure 2 shows the initial density field and its time evolution.

Data statistics are plotted in fig. 3 where the time evolution of integrated kinetic energy per unit of volume, isotropy degree [33] (see appendix B), microscale Reynolds number, density field standard deviation in space, velocity longitudinal-derivative skewness (the expression of the skewness is recalled in the next section), and eddy turnover time are displayed. The standard deviation of the density field $\sigma(\rho)$ is computed as the square root of the average of the squared deviations from the mean of the density field. On one hand, the figure shows that the adjusted forcing in time is able to guarantee integrated kinetic energy conservation and isotropy. On the other hand, the figure shows that we are able to select a specific time interval of approximately 2 eddy turnover times, where the turbulent Reynolds number as well as the skewness factor keep an almost constant value while mass diffusion effects are still significant. In the following, this time interval used for the statistical study is noted $\Gamma$.

## B. Density field mixing

The introduction of an inhomogeneous density field into a turbulent flow initiates its mixing, facilitated by the turbulent eddies. PDFs can be employed to characterize the state of mixing and its evolution. fig. 4 displays the PDFs of the density field at various selected times for different simulations with distinct initial density ratio $s$. The PDFs initially manifest as a symmetric double peak, corresponding to two initially pure fluids connected by a transition region. As the fluids begin to mix during the simulation, the PDFs become asymmetrical. A progressive stretching of the bell curve towards the left is evident, indicating a positive asymmetry. This implies that the denser fluid component mixes more rapidly than the less dense component. This behavior is also observed in turbulence simulations driven by buoyancy effects conducted by Livescu and Ristorcelli [1].

It is also noteworthy that the differences between the various density ratios $s$, are minimal, Especially concerning the asymmetry, which does not intensify with the increasing in $s$. This asymmetry is substantiated by the skewness
<figure><img src="image_2.png" /><figcaption>FIG. 1: Energy spectra at final time $t = t_f$ of constant-density HIT run 256-1 used as initial condition for the variable-density HIT simulations. The vertical line delimits the mode associated with the aliasing filter. Simulation parameters for the constant-density HIT DNS are listed in table II page 5.</figcaption></figure>

where $r = \|x\|_2$ and $r_0 = (3\pi^2)^{1/3}$. The parameter $d$ is the length of the transition region between the interior of the sphere and the surrounding fluid, and is fixed to 10 mesh cells. One may recall that the parameter $s$ represents the density ratio between the fluid filling the sphere and the fluid surrounding it.

In simpler terms, simulations are initialized by the injection of an inhomogeneous density field into a turbulent flow. The abrupt injection of the density field is mitigated by applying the filter described in eq. (8) during the initial 100 iterations of the simulation. This filtering process spans a time interval equivalent to one-tenth of an eddy turnover time.

One can specify that using a turbulent field as initial condition is necessary since injecting the variable-density field before forced turbulence is fully developed will result in complete mixing of the mass field before reaching a statistically stationary turbulent state.

In addition, it is also essential to select a specific time window for the study that encompasses mass diffusion effects while ensuring that the turbulent fields can reach an almost statistically stationary state.

For illustration purposes, the run 256-4 is now considered. Figure 2 shows the initial density field and its time evolution.

Data statistics are plotted in fig. 3 where the time evolution of integrated kinetic energy per unit of volume, isotropy degree [33] (see appendix B), microscale Reynolds number, density field standard deviation in space, velocity longitudinal-derivative skewness (the expression of the skewness is recalled in the next section), and eddy turnover time are displayed. The standard deviation of the density field $\sigma(\rho)$ is computed as the square root of the average of the squared deviations from the mean of the density field. On one hand, the figure shows that the adjusted forcing in time is able to guarantee integrated kinetic energy conservation and isotropy. On the other hand, the figure shows that we are able to select a specific time interval of approximately 2 eddy turnover times, where the turbulent Reynolds number as well as the skewness factor keep an almost constant value while mass diffusion effects are still significant. In the following, this time interval used for the statistical study is noted $\Gamma$.

## B. Density field mixing

The introduction of an inhomogeneous density field into a turbulent flow initiates its mixing, facilitated by the turbulent eddies. PDFs can be employed to characterize the state of mixing and its evolution. fig. 4 displays the PDFs of the density field at various selected times for different simulations with distinct initial density ratio $s$. The PDFs initially manifest as a symmetric double peak, corresponding to two initially pure fluids connected by a transition region. As the fluids begin to mix during the simulation, the PDFs become asymmetrical. A progressive stretching of the bell curve towards the left is evident, indicating a positive asymmetry. This implies that the denser fluid component mixes more rapidly than the less dense component. This behavior is also observed in turbulence simulations driven by buoyancy effects conducted by Livescu and Ristorcelli [1].

It is also noteworthy that the differences between the various density ratios $s$, are minimal, Especially concerning the asymmetry, which does not intensify with the increasing in $s$. This asymmetry is substantiated by the skewness
<figure><img src="image_2.png" /><figcaption>FIG. 2: xz-Cross-section plot of the density field at different times during simulation run 256-4. Re = 1000, s = 4, Sc = 1, N = 256.</figcaption></figure>

coefficient of the density fluctuations

$$S_{\rho} = \frac{\langle (\rho - \bar{\rho})^3 \rangle_{\Omega}}{\langle (\rho - \bar{\rho})^2 \rangle_{\Omega}^{3/2}}, \qquad (19)$$

which remains positive throughout the simulation as depicted in fig. 5. The evolution of $S_{\rho}$ reveals that the asymmetry diminishes as mixing progresses. However, towards the end of the simulation, it appears to reach a plateau, with a skewness coefficient of approximately 1, which remains far from a Gaussian behavior, characterized by a zero coefficient.

The evolution of the flatness coefficient of the density fluctuations,

$$F_{\rho} = \frac{\langle (\rho - \bar{\rho})^4 \rangle_{\Omega}}{\langle (\rho - \bar{\rho})^2 \rangle_{\Omega}^2}, \qquad (20)$$

also presented in fig. 5, demonstrate a gradual convergence towards Gaussian behavior ($F = 3$) as mixing proceeds.

The mixing can also be characterized by the PDFs of the gradient of the density field. Figures 6 and 7 illustrate the PDFs of the longitudinal $\partial\rho/\partial x$ and transversal $\partial\rho/\partial y$ derivatives of the density field. Initially, the inhomogeneous density field comprises regions of pure fluids where the directional density derivatives are zero, except for the transition zone between the two regions. This results in a PDF where values close to zero have the highest frequency, while the high-amplitude values occur with low frequency.

During the initial moments of the simulation, the PDF broadens due to the emergence of multiple smaller-scale mixing zones. This widening suggests that variable-density effects are rapidly transported to smaller scales during mixing. Subsequently, over time, the PDFs narrow as gradients gradually dissipate due to diffusion.
<figure><img src="image_3.png" /><figcaption>FIG. 2: xz-Cross-section plot of the density field at different times during simulation run 256-4. Re = 1000, s = 4, Sc = 1, N = 256.</figcaption></figure>

coefficient of the density fluctuations

$$S_{\rho} = \frac{\langle (\rho - \bar{\rho})^3 \rangle_{\Omega}}{\langle (\rho - \bar{\rho})^2 \rangle_{\Omega}^{3/2}}, \qquad (19)$$

which remains positive throughout the simulation as depicted in fig. 5. The evolution of $S_{\rho}$ reveals that the asymmetry diminishes as mixing progresses. However, towards the end of the simulation, it appears to reach a plateau, with a skewness coefficient of approximately 1, which remains far from a Gaussian behavior, characterized by a zero coefficient.

The evolution of the flatness coefficient of the density fluctuations,

$$F_{\rho} = \frac{\langle (\rho - \bar{\rho})^4 \rangle_{\Omega}}{\langle (\rho - \bar{\rho})^2 \rangle_{\Omega}^2}, \qquad (20)$$

also presented in fig. 5, demonstrate a gradual convergence towards Gaussian behavior ($F = 3$) as mixing proceeds.

The mixing can also be characterized by the PDFs of the gradient of the density field. Figures 6 and 7 illustrate the PDFs of the longitudinal $\partial\rho/\partial x$ and transversal $\partial\rho/\partial y$ derivatives of the density field. Initially, the inhomogeneous density field comprises regions of pure fluids where the directional density derivatives are zero, except for the transition zone between the two regions. This results in a PDF where values close to zero have the highest frequency, while the high-amplitude values occur with low frequency.

During the initial moments of the simulation, the PDF broadens due to the emergence of multiple smaller-scale mixing zones. This widening suggests that variable-density effects are rapidly transported to smaller scales during mixing. Subsequently, over time, the PDFs narrow as gradients gradually dissipate due to diffusion.
<figure><img src="image_3.png" /><figcaption>FIG. 3: Evolution of the volumetric total kinetic energy, isotropy degree [33] (calculation detailed in appendix B), microscale Reynolds number $Re_λ$, density field standard deviation in space $σ(ρ)$, velocity longitudinal-derivative skewness $S$, and eddy turnover time $T$ during simulation run 256-4. Vertical lines delimit a time interval $Γ$ of 2 eddy turnover times, where turbulent parameters do not fluctuate much, and where mass diffusion effects are still significant. $Re = 1000$, $s = 4$, $Sc = 1$, $N = 256$.</figcaption></figure>
<figure><img src="image_4.png" /><figcaption>FIG. 3: Evolution of the volumetric total kinetic energy, isotropy degree [33] (calculation detailed in appendix B), microscale Reynolds number $Re_λ$, density field standard deviation in space $σ(ρ)$, velocity longitudinal-derivative skewness $S$, and eddy turnover time $T$ during simulation run 256-4. Vertical lines delimit a time interval $Γ$ of 2 eddy turnover times, where turbulent parameters do not fluctuate much, and where mass diffusion effects are still significant. $Re = 1000$, $s = 4$, $Sc = 1$, $N = 256$.</figcaption></figure>
<figure><img src="image_4.png" /><figcaption>FIG. 4: Density field PDF evolution during DNSs of variable-density HIT. Simulation parameters for variable-density HIT DNSs are listed in table I page 5.</figcaption></figure>
<figure><img src="image_5.png" /><figcaption>FIG. 4: Density field PDF evolution during DNSs of variable-density HIT. Simulation parameters for variable-density HIT DNSs are listed in table I page 5.</figcaption></figure>
<figure><img src="image_5.png" /><figcaption>FIG. 5: Density field moments evolution during DNSs of variable-density HIT: (left) skewness defined in eq. (19) ; (right) flatness defined in eq. (20). Simulation parameters for variable-density HIT DNSs are listed in table I page 5.</figcaption></figure>

It is important to note that the distribution of gradients is symmetrical. This fully isotropic dynamics contrast with that observed by Livescu and Ristorcelli [1], where the PDFs of derivatives in the buoyancy direction exhibit asymmetry. It is also important to highlight that the density ratio *s* does not influence the distribution of gradients, except for the production of more intense gradients.

C. Statistical analysis

Small-scale turbulence intermittency requires that the PDFs exhibit flared-out tails [14]. In particular, Jiménez et al. [10], Gotoh et al. [34] and Ishihara et al. [13] demonstrate that in constant-density flows, the PDFs of the velocity gradients $\partial u/\partial x$, $\partial v/\partial x$ and of the vorticity components $\omega_x$ and $\omega_y$ become increasingly non-Gaussian, and feature wide tails as the turbulent Reynolds number $Re_\lambda$ increases. In the variable-density flows, the inhomogeneity of the density field has a direct influence onto the turbulence scales and dynamics as reported in fig. 8, where the temporal evolutions of $Re_\lambda$ are displayed for various initial density ratios. The figure shows that for all flow cases, $Re_\lambda$ reaches a nearly constant value for $2 < t < 6$, which corresponds to a time interval of 2 to 3 eddy turnover times. We also observe that the microscale Reynolds number is decreasing with $s$.

Figure 9 shows PDFs of the longitudinal velocity-derivative $\partial u/\partial x$, the transversal velocity-derivative $\partial v/\partial x$, vorticity components $\omega_x$ and $\omega_y$, and Lagrangian acceleration component $A_x \equiv (\partial/\partial t + (u \cdot \nabla))u$ for the different values of $s$ that are considered in this study. Following the discussion in section III A, quantities are averaged over the 500 time-steps corresponding to the time interval $\Gamma$. In fig. 9(a) and 9(b) the PDFs of $\partial u/\partial x$ are depicted with and without normalization by the standard deviation, respectively. In panel (b), we observe that the PDFs for variable-density flows have wide tails, that become wider for larger $s$-values. This suggests, that as $s$ increases, energy is transferred to smaller-scale eddies with high velocity gradients. In panel (a), the normalized PDFs do not exhibit increasingly flare-out tails for $s \le 4$. In this case, the observed differences for the smallest scales, i.e. high velocity gradients, should be mainly attributed to a Reynolds number effect. For $s = 8$, the PDF becomes symmetrical and presents a flare-out tail. Hence, it seems to indicate that small-case intermittency remains unaffected by mass diffusion effects until a threshold value of $s$ between 4 and 8, given our control parameters. Beyond this threshold, the widening of the PDFs describes an increase in $\partial u/\partial x$ small-case intermittency. Intermittency of the transverse velocity-derivative and vorticity components remains unchanged regardless of the density ratio as shown in fig. 9(c) to 9(e).

As mentioned by Ishihara et al. [13], turbulence intermittency occurs not only in space but also in time. Therefore, it is observed not only in velocity gradients but also in the time-derivatives of turbulent velocities. In fig. 9(f), a similar behavior to that of $\partial u/\partial x$ is seen in the PDF of the $A_x$ component of the Lagrangian acceleration: for $s \le 4$, the PDF do not exhibit increasing flare-out tails, and for $s = 8$, the PDF becomes wider and more symmetrical. Thus, both the space- and time-intermittency remains unaffected by variable-density effects until a specific threshold value
<figure><img src="image_6.png" /><figcaption>FIG. 5: Density field moments evolution during DNSs of variable-density HIT: (left) skewness defined in eq. (19) ; (right) flatness defined in eq. (20). Simulation parameters for variable-density HIT DNSs are listed in table I page 5.</figcaption></figure>

It is important to note that the distribution of gradients is symmetrical. This fully isotropic dynamics contrast with that observed by Livescu and Ristorcelli [1], where the PDFs of derivatives in the buoyancy direction exhibit asymmetry. It is also important to highlight that the density ratio *s* does not influence the distribution of gradients, except for the production of more intense gradients.

C. Statistical analysis

Small-scale turbulence intermittency requires that the PDFs exhibit flared-out tails [14]. In particular, Jiménez et al. [10], Gotoh et al. [34] and Ishihara et al. [13] demonstrate that in constant-density flows, the PDFs of the velocity gradients $\partial u/\partial x$, $\partial v/\partial x$ and of the vorticity components $\omega_x$ and $\omega_y$ become increasingly non-Gaussian, and feature wide tails as the turbulent Reynolds number $Re_\lambda$ increases. In the variable-density flows, the inhomogeneity of the density field has a direct influence onto the turbulence scales and dynamics as reported in fig. 8, where the temporal evolutions of $Re_\lambda$ are displayed for various initial density ratios. The figure shows that for all flow cases, $Re_\lambda$ reaches a nearly constant value for $2 < t < 6$, which corresponds to a time interval of 2 to 3 eddy turnover times. We also observe that the microscale Reynolds number is decreasing with $s$.

Figure 9 shows PDFs of the longitudinal velocity-derivative $\partial u/\partial x$, the transversal velocity-derivative $\partial v/\partial x$, vorticity components $\omega_x$ and $\omega_y$, and Lagrangian acceleration component $A_x \equiv (\partial/\partial t + (u \cdot \nabla))u$ for the different values of $s$ that are considered in this study. Following the discussion in section III A, quantities are averaged over the 500 time-steps corresponding to the time interval $\Gamma$. In fig. 9(a) and 9(b) the PDFs of $\partial u/\partial x$ are depicted with and without normalization by the standard deviation, respectively. In panel (b), we observe that the PDFs for variable-density flows have wide tails, that become wider for larger $s$-values. This suggests, that as $s$ increases, energy is transferred to smaller-scale eddies with high velocity gradients. In panel (a), the normalized PDFs do not exhibit increasingly flare-out tails for $s \le 4$. In this case, the observed differences for the smallest scales, i.e. high velocity gradients, should be mainly attributed to a Reynolds number effect. For $s = 8$, the PDF becomes symmetrical and presents a flare-out tail. Hence, it seems to indicate that small-case intermittency remains unaffected by mass diffusion effects until a threshold value of $s$ between 4 and 8, given our control parameters. Beyond this threshold, the widening of the PDFs describes an increase in $\partial u/\partial x$ small-case intermittency. Intermittency of the transverse velocity-derivative and vorticity components remains unchanged regardless of the density ratio as shown in fig. 9(c) to 9(e).

As mentioned by Ishihara et al. [13], turbulence intermittency occurs not only in space but also in time. Therefore, it is observed not only in velocity gradients but also in the time-derivatives of turbulent velocities. In fig. 9(f), a similar behavior to that of $\partial u/\partial x$ is seen in the PDF of the $A_x$ component of the Lagrangian acceleration: for $s \le 4$, the PDF do not exhibit increasing flare-out tails, and for $s = 8$, the PDF becomes wider and more symmetrical. Thus, both the space- and time-intermittency remains unaffected by variable-density effects until a specific threshold value
<figure><img src="image_6.png" /><figcaption>FIG. 6: $\partial\rho/\partial x$ PDF evolution during DNSs of variable-density HIT. Simulation parameters for variable-density HIT DNSs are listed in table I page 5.</figcaption></figure>
<figure><img src="image_7.png" /><figcaption>FIG. 6: $\partial\rho/\partial x$ PDF evolution during DNSs of variable-density HIT. Simulation parameters for variable-density HIT DNSs are listed in table I page 5.</figcaption></figure>
<figure><img src="image_7.png" /><figcaption>FIG. 7: $\partial\rho/\partial y$ PDF evolution during DNSs of variable-density HIT. Simulation parameters for variable-density HIT DNSs are listed in table I page 5.</figcaption></figure>
<figure><img src="image_8.png" /><figcaption>FIG. 7: $\partial\rho/\partial y$ PDF evolution during DNSs of variable-density HIT. Simulation parameters for variable-density HIT DNSs are listed in table I page 5.</figcaption></figure>
<figure><img src="image_8.png" /><figcaption>FIG. 8: Time evolution of the turbulent Reynolds number $Re_λ = (λ/L) Re$ during variable-density simulations. DNS parameters are listed in table I. For comparison, dashed line depicts the evolution throughout constant-density simulation 256-1 for $t \in [40, 50]$. Vertical lines delimit the study time interval $Γ$ which is considered for the statistical computations.</figcaption></figure>

of the density ratio. Beyond this limit, an increase in intermittency is observed.

The non-Gaussian statistics of the velocity increments may also be characterized quantitatively by the third- and fourth-order moments of the longitudinal velocity-derivative, the skewness $S$ and flatness $F$ factors, respectively defined as

$$S = \frac{\langle (\partial u / \partial x)^3 \rangle}{\langle (\partial u / \partial x)^2 \rangle^{3/2}}, \qquad (21)$$

and

$$F = \frac{\langle (\partial u / \partial x)^4 \rangle}{\langle (\partial u / \partial x)^2 \rangle^2}. \qquad (22)$$

For constant-density flows, Ishihara et al. [13] and Gylfason et al. [35] show that $S$ remains constant for $Re_λ < 200$, and decreases slowly for larger $Re_λ$. The authors propose the scaling law

$$-S \approx (0.32 \mp 0.02) Re_{\lambda}^{0.11 \pm 0.01}. \qquad (23)$$

Table I presents the skewness factor averaged over the statistic study interval $Γ$ for different $s$-values. Notably, the skewness factor converges toward different values depending on $s$, even if the microscale Reynolds number remains below 200 for all cases examined. For instance, for $s = 2$, an averaged skewness factor of $S = -0.43$ is obtained while for $s = 4$, the skewness factor has increased up to $S = -0.30$. Thus, the assumption of a constant skewness factor for $Re_λ < 200$ does not hold for variable-density flows. Moreover, we note that an increase in $s$ tends also to increase the skewness factor. This suggests that PDFs would become progressively more symmetrical as the density ratio is increasing. This characteristic is even more obvious in fig. 9(a), where the PDF of $∂u/∂x$ appears more symmetrical for $s = 8$ compared to smaller values of $s$. Finally, the scaling law provided by eq. (23) tends to overestimate the absolute value of the skewness factor in the case of variable-density flows, as can be observed in fig. 10.

Prior investigations on constant-density flows have demonstrated that the flatness factor exhibits an increasing trend as the microscale Reynolds number increases for $100 < Re_λ < 1000$. Ishihara et al. [13] propose the scaling law

$$F \approx (1.14 \pm 0.19) Re_{\lambda}^{0.34 \pm 0.03}. \qquad (24)$$

While our findings align with an increasing trend of $F$ with respect to $Re_λ$, (see fig. 10) our dataset is not large enough to draw definitive conclusions.

## D. Spectral analysis

From previous statistical investigations, the modification of energy distributions over the different scales has been clearly illustrated as variable-density effects increase. In an effort to explain this behavior, this section is devoted to
<figure><img src="image_9.png" /><figcaption>FIG. 8: Time evolution of the turbulent Reynolds number $Re_λ = (λ/L) Re$ during variable-density simulations. DNS parameters are listed in table I. For comparison, dashed line depicts the evolution throughout constant-density simulation 256-1 for $t \in [40, 50]$. Vertical lines delimit the study time interval $Γ$ which is considered for the statistical computations.</figcaption></figure>

of the density ratio. Beyond this limit, an increase in intermittency is observed.

The non-Gaussian statistics of the velocity increments may also be characterized quantitatively by the third- and fourth-order moments of the longitudinal velocity-derivative, the skewness $S$ and flatness $F$ factors, respectively defined as

$$S = \frac{\langle (\partial u / \partial x)^3 \rangle}{\langle (\partial u / \partial x)^2 \rangle^{3/2}}, \qquad (21)$$

and

$$F = \frac{\langle (\partial u / \partial x)^4 \rangle}{\langle (\partial u / \partial x)^2 \rangle^2}. \qquad (22)$$

For constant-density flows, Ishihara et al. [13] and Gylfason et al. [35] show that $S$ remains constant for $Re_λ < 200$, and decreases slowly for larger $Re_λ$. The authors propose the scaling law

$$-S \approx (0.32 \mp 0.02) Re_{\lambda}^{0.11 \pm 0.01}. \qquad (23)$$

Table I presents the skewness factor averaged over the statistic study interval $Γ$ for different $s$-values. Notably, the skewness factor converges toward different values depending on $s$, even if the microscale Reynolds number remains below 200 for all cases examined. For instance, for $s = 2$, an averaged skewness factor of $S = -0.43$ is obtained while for $s = 4$, the skewness factor has increased up to $S = -0.30$. Thus, the assumption of a constant skewness factor for $Re_λ < 200$ does not hold for variable-density flows. Moreover, we note that an increase in $s$ tends also to increase the skewness factor. This suggests that PDFs would become progressively more symmetrical as the density ratio is increasing. This characteristic is even more obvious in fig. 9(a), where the PDF of $∂u/∂x$ appears more symmetrical for $s = 8$ compared to smaller values of $s$. Finally, the scaling law provided by eq. (23) tends to overestimate the absolute value of the skewness factor in the case of variable-density flows, as can be observed in fig. 10.

Prior investigations on constant-density flows have demonstrated that the flatness factor exhibits an increasing trend as the microscale Reynolds number increases for $100 < Re_λ < 1000$. Ishihara et al. [13] propose the scaling law

$$F \approx (1.14 \pm 0.19) Re_{\lambda}^{0.34 \pm 0.03}. \qquad (24)$$

While our findings align with an increasing trend of $F$ with respect to $Re_λ$, (see fig. 10) our dataset is not large enough to draw definitive conclusions.

## D. Spectral analysis

From previous statistical investigations, the modification of energy distributions over the different scales has been clearly illustrated as variable-density effects increase. In an effort to explain this behavior, this section is devoted to
<figure><img src="image_9.png" /><figcaption>FIG. 9: $\zeta$-PDFs for (a), (b) $\zeta = \partial u/\partial x$, (c) $\zeta = \partial v/\partial x$, (d) $\zeta = \omega_x$, (e) $\zeta = \omega_y$, (f) $\zeta = A_x$ (Lagrangian acceleration), averaged over the study time interval $\Gamma$ of variable-density simulations, and compared with the constant-density initial condition (dashed). DNS parameters are listed in table I. $\sigma$ is the standard deviation of $\zeta$ in each PDF. $\langle Re_\lambda \rangle$ represents the averaged turbulent Reynolds number over time interval $\Gamma$. In (b), the PDF has not been normalized by $\sigma$.</figcaption></figure>
<figure><img src="image_10.png" /><figcaption>FIG. 9: $\zeta$-PDFs for (a), (b) $\zeta = \partial u/\partial x$, (c) $\zeta = \partial v/\partial x$, (d) $\zeta = \omega_x$, (e) $\zeta = \omega_y$, (f) $\zeta = A_x$ (Lagrangian acceleration), averaged over the study time interval $\Gamma$ of variable-density simulations, and compared with the constant-density initial condition (dashed). DNS parameters are listed in table I. $\sigma$ is the standard deviation of $\zeta$ in each PDF. $\langle Re_\lambda \rangle$ represents the averaged turbulent Reynolds number over time interval $\Gamma$. In (b), the PDF has not been normalized by $\sigma$.</figcaption></figure>
<figure><img src="image_10.png" /><figcaption>FIG. 10: Skewness (S) and flatness (F) coefficients from variable-density HIT DNSs compared with constant-density HIT scaling laws (dashed), eqs. (23) and (24), proposed by Ishihara et al. [13].</figcaption></figure>

give some physical insights from a structural point of view using Fourier spectral analysis and flow-fields visualization.

The first difference is highlighted in fig. 11 where cross-sections of vorticity norm isosurfaces for both constant- and variable-density flows (referenced 256-1 and 256-4, respectively) are depicted. It shows that the variable-density flow exhibits smaller coherent structures than constant-density flows.

Let us now further examine the mechanisms behind the generation of these motions. For that purpose, the spectral density of energy is represented in fig. 12 for the variable-density simulations and its constant-density counterpart. From the figure, it is obvious that the Kolmogorov slope remains unchanged in the larger scales; however, there is a discernible alternation in energy levels within the range of the smallest scales. An analogous pattern can also be observed when calculating the spectrum based on the kinetic energy per unit of volume (i.e. $\rho u$). Indeed, fig. 13 presents comparison between spectra calculated traditionally using the energy per unit of mass and spectra calculated using energy per unit of volume (dotted lines). Both spectra exhibit the Kolmogorov cascade, and a noticeable increase in energy can be observed towards the end of the inertial range. It then suggests an additional mechanism associated with the emergence of these smallest scales.

In summary, the spectra of the variable-density HIT exhibit an energy bump in the smaller scales of the inertial range that is absent in the constant density spectra. This modification confirms the presence of additional dynamics introduced by mass effects, which become dominant beyond a certain structure size. This observation completes the analysis of the velocity PDFs, which reveals a greater distribution of energetic structures in the variable-density flow compared to the constant-density flow. The spectral analysis indicates that this increase in energy distribution is associated with an energy gain within the smallest structures of the flow. It appears that in a variable-density HIT flow, the final scale of the inertial range, which in the constant-density case is related to viscous diffusion, exhibits an energy transfer dynamic induced by mass variations.

As a result, viscous dissipation emerges as the predominant process at smaller scales compared to the constant-density case. This outcome is also illustrated by the reduction of the Taylor and Kolmogorov scales as a function of the density ratio $s$, as shown in table I and fig. 14. In this figure, the vertical line representing the Taylor scales of the variable-density simulation gradually moves away from the constant-density lines as $s$ increases. For $s < 1$, the Taylor scale increases, and the viscous cutoff appears earlier in the spectrum; in constant-density DNS of HIT, the reduction in scales and the shifting of the viscous cutoff to smaller scales are the consequence of an increase in the Reynolds number. Therefore, increasing $s$ leads to a phenomenon similar to increasing the Reynolds number $Re$. This observation aligns with the results of the PDF analysis, where an increase in $s$ results in broadening, which is generally associated with an increase in $Re$.

This behavior can be deduced from the momentum equation eq. (1a), where the viscous dissipation term is transformed, upon introducing the variable-density field, from $Re^{-1} \nabla^2 u$ to

$$ (\rho Re)^{-1} \nabla^2 u \approx \left( \frac{s-1}{2} Re \right)^{-1} \nabla^2 u \equiv Re'^{-1} \nabla^2 u, \qquad (25) $$

where $\frac{1}{2}(s-1)$ represents the mean value of the density field. Thus, for $s > 1$, the introduction of the variable-density field leads to an increase in the equivalent Reynolds number $Re'$ and, consequently, a decrease in the viscous length
<figure><img src="image_11.png" /><figcaption>FIG. 10: Skewness (S) and flatness (F) coefficients from variable-density HIT DNSs compared with constant-density HIT scaling laws (dashed), eqs. (23) and (24), proposed by Ishihara et al. [13].</figcaption></figure>

give some physical insights from a structural point of view using Fourier spectral analysis and flow-fields visualization.

The first difference is highlighted in fig. 11 where cross-sections of vorticity norm isosurfaces for both constant- and variable-density flows (referenced 256-1 and 256-4, respectively) are depicted. It shows that the variable-density flow exhibits smaller coherent structures than constant-density flows.

Let us now further examine the mechanisms behind the generation of these motions. For that purpose, the spectral density of energy is represented in fig. 12 for the variable-density simulations and its constant-density counterpart. From the figure, it is obvious that the Kolmogorov slope remains unchanged in the larger scales; however, there is a discernible alternation in energy levels within the range of the smallest scales. An analogous pattern can also be observed when calculating the spectrum based on the kinetic energy per unit of volume (i.e. $\rho u$). Indeed, fig. 13 presents comparison between spectra calculated traditionally using the energy per unit of mass and spectra calculated using energy per unit of volume (dotted lines). Both spectra exhibit the Kolmogorov cascade, and a noticeable increase in energy can be observed towards the end of the inertial range. It then suggests an additional mechanism associated with the emergence of these smallest scales.

In summary, the spectra of the variable-density HIT exhibit an energy bump in the smaller scales of the inertial range that is absent in the constant density spectra. This modification confirms the presence of additional dynamics introduced by mass effects, which become dominant beyond a certain structure size. This observation completes the analysis of the velocity PDFs, which reveals a greater distribution of energetic structures in the variable-density flow compared to the constant-density flow. The spectral analysis indicates that this increase in energy distribution is associated with an energy gain within the smallest structures of the flow. It appears that in a variable-density HIT flow, the final scale of the inertial range, which in the constant-density case is related to viscous diffusion, exhibits an energy transfer dynamic induced by mass variations.

As a result, viscous dissipation emerges as the predominant process at smaller scales compared to the constant-density case. This outcome is also illustrated by the reduction of the Taylor and Kolmogorov scales as a function of the density ratio $s$, as shown in table I and fig. 14. In this figure, the vertical line representing the Taylor scales of the variable-density simulation gradually moves away from the constant-density lines as $s$ increases. For $s < 1$, the Taylor scale increases, and the viscous cutoff appears earlier in the spectrum; in constant-density DNS of HIT, the reduction in scales and the shifting of the viscous cutoff to smaller scales are the consequence of an increase in the Reynolds number. Therefore, increasing $s$ leads to a phenomenon similar to increasing the Reynolds number $Re$. This observation aligns with the results of the PDF analysis, where an increase in $s$ results in broadening, which is generally associated with an increase in $Re$.

This behavior can be deducted from the momentum equation eq. (1a), where the viscous dissipation term is transformed, upon introducing the variable-density field, from $Re^{-1} \nabla^2 u$ to

$$ (\rho Re)^{-1} \nabla^2 u \approx \left( \frac{s-1}{2} Re \right)^{-1} \nabla^2 u \equiv Re'^{-1} \nabla^2 u, \quad (25) $$

where $\frac{1}{2}(s-1)$ represents the mean value of the density field. Thus, for $s > 1$, the introduction of the variable-density field leads to an increase in the equivalent Reynolds number $Re'$ and, consequently, a decrease in the viscous length
<figure><img src="image_11.png" /><figcaption>FIG. 11: xz-Cross-section plot of the vorticity norm for both (a), (b) the constant-density initial condition and (c), (d) the variable-density simulation 256-4 (at time t = 4). Re = 1000, s = 4, Sc = 1, N = 256.</figcaption></figure>

scale. Conversely, for $s < 1$, $Re'$ decreases, resulting in an increase in the viscous length scale. This increase in the viscous length scale explains the observed behavior in the simulation with $s < 1$.

It is important to recall that the introduction of a density field where the lighter fluid is confined within the sphere for the initial condition has little influence on turbulence dynamics. This is reflected in velocity field PDFs that are less spread out than those of the constant-density simulation, as well as spectra exhibiting a viscous cutoff at larger scales. Consequently, in this configuration, the viscous cutoff occurs before the range of scale where variable-density introduce the energy peak. Therefore, it has a limited impact on turbulent dynamics. It is essential to note that the flow does not exhibit $s \to s^{-1}$ symmetry since the mean values of $\rho$ are not equal, resulting in different equivalent $Re$ numbers and, consequently, different dynamics.

It is widely recognized that density variations can give rise to various phenomena, with the Rayleigh-Taylor instability (RTI) being one of the most prominent examples [36–39]. Now, we will investigate whether the emergence of higher energy structures within the lower region of the spectrum is a result of vortical structures developing under the influence of the Rayleigh-Taylor instability. Recently, Jacques et al. [40] prove that $R < 0$, with

$$R = \text{local acceleration} \cdot \text{local mass variation} = ((u \cdot \nabla)u) \cdot \nabla \rho, \qquad (26)$$

is a sufficient condition for RTI in rotating flows. Figure 15 displays a cross-section of the RTI criterion, presenting both unfiltered and filtered flow fields. When comparing it with cross-sectional view of the vorticity depicted in fig. 11(d), it becomes evident that the RTI criterion is satisfied for similar flow structures. Finally, fig. 16 illustrates a comparison between the spectra and the RTI criterion (26) spectral density of energy

$$\text{ESD } R = |\mathcal{F}R|^2. \qquad (27)$$
<figure><img src="image_12.png" /><figcaption>FIG. 11: xz-Cross-section plot of the vorticity norm for both (a), (b) the constant-density initial condition and (c), (d) the variable-density simulation 256-4 (at time t = 4). Re = 1000, s = 4, Sc = 1, N = 256.</figcaption></figure>

scale. Conversely, for $s < 1$, $Re'$ decreases, resulting in an increase in the viscous length scale. This increase in the viscous length scale explains the observed behavior in the simulation with $s < 1$.

It is important to recall that the introduction of a density field where the lighter fluid is confined within the sphere for the initial condition has little influence on turbulence dynamics. This is reflected in velocity field PDFs that are less spread out than those of the constant-density simulation, as well as spectra exhibiting a viscous cutoff at larger scales. Consequently, in this configuration, the viscous cutoff occurs before the range of scale where variable-density introduce the energy peak. Therefore, it has a limited impact on turbulent dynamics. It is essential to note that the flow does not exhibit $s \to s^{-1}$ symmetry since the mean values of $\rho$ are not equal, resulting in different equivalent $Re$ numbers and, consequently, different dynamics.

It is widely recognized that density variations can give rise to various phenomena, with the Rayleigh-Taylor instability (RTI) being one of the most prominent examples [36-39]. Now, we will investigate whether the emergence of higher energy structures within the lower region of the spectrum is a result of vortical structures developing under the influence of the Rayleigh-Taylor instability. Recently, Jacques et al. [40] prove that $R < 0$, with

$$R = \text{local acceleration} \cdot \text{local mass variation} = ((\mathbf{u} \cdot \nabla) \mathbf{u}) \cdot \nabla \rho, \qquad (26)$$

is a sufficient condition for RTI in rotating flows. Figure 15 displays a cross-section of the RTI criterion, presenting both unfiltered and filtered flow fields. When comparing it with cross-sectional view of the vorticity depicted in fig. 11(d), it becomes evident that the RTI criterion is satisfied for similar flow structures. Finally, fig. 16 illustrates a comparison between the spectra and the RTI criterion (26) spectral density of energy

$$\text{ESD } R = |\mathcal{F}R|^2. \qquad (27)$$
<figure><img src="image_12.png" /><figcaption>FIG. 12: Spectra averaged over the study time interval $Γ$ of variable-density simulations, and compared with the constant-density initial condition (dashed). Simulation parameters for variable-density HIT DNSs are listed in table I page 5.</figcaption></figure>
<figure><img src="image_13.png" /><figcaption>FIG. 12: Spectra averaged over the study time interval $Γ$ of variable-density simulations, and compared with the constant-density initial condition (dashed). Simulation parameters for variable-density HIT DNSs are listed in table I page 5.</figcaption></figure>
<figure><img src="image_13.png" /><figcaption>FIG. 13: Spectra averaged over the study time interval $Γ$ of variable-density simulations, and compared with the constant-density initial condition (dashed). Dotted lines correspond to spectra $E_ρ$ computed with kinetic energy per unit of volume $ρu$. Simulation parameters for variable-density HIT DNSs are listed in table I page 5.</figcaption></figure>

<figure><img src="image_14.png" /><figcaption>FIG. 14: Spectra averaged over the study time interval $Γ$ of variable-density simulations, and compared with the constant-density initial condition (dashed). The vertical lines represent the integral and Taylor scales associated with the equations defined in eqs. (16) and (17). Simulation parameters for variable-density HIT DNSs are listed in table I page 5.</figcaption></figure>
<figure><img src="image_14.png" /><figcaption>FIG. 13: Spectra averaged over the study time interval $Γ$ of variable-density simulations, and compared with the constant-density initial condition (dashed). Dotted lines correspond to spectra $E_ρ$ computed with kinetic energy per unit of volume $ρu$. Simulation parameters for variable-density HIT DNSs are listed in table I page 5.</figcaption></figure>

<figure><img src="image_15.png" /><figcaption>FIG. 14: Spectra averaged over the study time interval $Γ$ of variable-density simulations, and compared with the constant-density initial condition (dashed). The vertical lines represent the integral and Taylor scales associated with the equations defined in eqs. (16) and (17). Simulation parameters for variable-density HIT DNSs are listed in table I page 5.</figcaption></figure>
<figure><img src="image_15.png" /><figcaption>FIG. 15: xz-Cross-section plot of the Rayleigh-Taylor instability criterion (R) at time t = 4 of simulation 256-4. Re = 1000, s = 4, Sc = 1, N = 256.</figcaption></figure>

<figure><img src="image_16.png" /><figcaption>FIG. 16: Spectra averaged over the study time interval $Γ$ of variable-density simulations, and compared with the spectral density of the Rayleigh-Taylor instability criterion R (dotted). DNS parameters are listed in table I.</figcaption></figure>

Notably, the region where the criterion exhibits the highest energy aligns with the range of modes where we observe the alterations in energy levels.

This study supports that differences observed in PDFs could be attributed to the generation of smallest scales structures under the effect of Rayleigh-Taylor-like mechanism.

## IV. CONCLUSIONS

In this paper, we investigate the effects of density variations on statistically stationary forced isotropic turbulence. Simulations are carried out by introducing a density field in the shape of a sphere into a homogeneous isotropic turbulent state obtained from a preliminary constant-density simulation. The preservation of the turbulent kinetic energy is achieved with the linear forcing method proposed by Lundgren, extended here to the incompressible variable-
<figure><img src="image_16.png" /><figcaption>FIG. 15: $xz$-Cross-section plot of the Rayleigh-Taylor instability criterion ($R$) at time $t = 4$ of simulation 256-4. $Re = 1000$, $s = 4$, $Sc = 1$, $N = 256$.</figcaption></figure>

<figure><img src="image_17.png" /><figcaption>FIG. 16: Spectra averaged over the study time interval $Γ$ of variable-density simulations, and compared with the spectral density of the Rayleigh-Taylor instability criterion $R$ (dotted). DNS parameters are listed in table I.</figcaption></figure>

Notably, the region where the criterion exhibits the highest energy aligns with the range of modes where we observe the alterations in energy levels.

This study supports that differences observed in PDFs could be attributed to the generation of smallest scales structures under the effect of Rayleigh-Taylor-like mechanism.

## IV. CONCLUSIONS

In this paper, we investigate the effects of density variations on statistically stationary forced isotropic turbulence. Simulations are carried out by introducing a density field in the shape of a sphere into a homogeneous isotropic turbulent state obtained from a preliminary constant-density simulation. The preservation of the turbulent kinetic energy is achieved with the linear forcing method proposed by Lundgren, extended here to the incompressible variable-
density equations. We explore density ratios ranging from $s = 0.5$ to $s = 8$, while maintaining a Reynolds number of 1000 and a Schmidt number of 1. Initially, it is noted that the microscale Reynolds number $Re_λ$ reaches different values depending on $s$. For a better understanding of variable-density effects, statistical analyses are conducted within a time interval of approximately 2 eddy turnover times. In this interval, $Re_λ$ and the skewness factor remain nearly constant, while the effects of mass diffusion stay significant. In order to mitigate the occurrence of transient effects resulting from the abrupt introduction of the inhomogeneous density field into the turbulent initial condition, a filter is applied during the first hundred iterations of the simulation. Additionally, the statistical analysis interval commences two eddy turnover times later than the start of the simulation.

Special attention is paid to the PDFs of velocity gradients and Lagrangian acceleration. DNS data show that PDFs for variable-density flows exhibit wide tails, which become even wider as $s$ increases. As a result, a noticeable increase in space- and time-intermittency is observed once a threshold $s$ value is exceeded. As quantitative measurements of the deviation from Gaussian distributions, we study the skewness $S$ and flatness factors $F$ of the longitudinal velocity-derivative for various $s$ values. Our findings indicate that the skewness factor increases with the density ratio, even for $Re_λ < 200$. Consequently, the PDFs gradually become more symmetrical as the value of $s$ increases. Additionally, the flatness factor shows an increasing trend with respect to $Re_λ$. These results are consistent with trends observed in constant-density scenarios. However, the $Re_λ$-power-laws found in the literature for homogeneous density fields tend to overestimate both $S$ and $F$ values.

Flow visualization reveals that variable-density flows exhibit smaller coherent structures compared to constant-density flows. Fourier analyses show there is a significant modification in energy levels within the smallest scales of the inertial range. In particular, the region where the energy levels are altered matches area where the Rayleigh-Taylor instability criterion peaks. It suggests that the introduction of variable-density effects into isotropic turbulence results in the generation of smaller eddy structures through Rayleigh-Taylor-like instability.

As a perspective, the logical continuation of this study should be devoted to the analysis of variable-density jets that undergo a transition from laminar to turbulent flow. In particular, density variations are known to shift the jet dynamics though subtle mechanisms. In that respect, the work of Di Pierro and Abid [41], Di Pierro and Abid [37], Ravier et al. [42], and Jacques et al. [40] dealing with inhomogeneous jets can be mentioned. Hence, it would be interesting to evaluate if statistics obtained within the framework of homogeneous isotropic turbulence can be compared to those computed from variable-density turbulent jet simulations. Specifically, one may wonder whether the alteration of energy levels in the small scales observed in the present study — which could be associated to RTI — is observed.

ACKNOWLEDGMENTS

This work was granted access to the HPC resources of PMCS2I (Pôle de Modélisation et de Calcul en Sciences de l'Ingénieur de l'Information) of École Centrale de Lyon, and P2CHPD (Pôle de Calcul Hautes Performances Dédiés) of Université Lyon1.

Appendix A: Relation between mean-mass and mean-volume velocity fields

The incompressible formulation serves as an efficient and self-contained tool for describing the physics of low-Mach number flows, whether they exhibit constant or variable-density. This formulation is based on the assumption that velocity field is solenoidal. However, in a mixing flow, there are multiple possible representations of the velocity field. Each species $k$ has a velocity $v_k$, defined as the Eulerian-averaged velocity of the $k$ species molecule. In a DNS approach aimed at resolving the entire flow within the study domain $\Omega$, the goal is to work with an ensemble velocity that encompasses the behavior of all these species. For modeling the ensemble velocity, two approaches are conceivable:

* The mean-mass velocity $\bar{u}$, averaged over the total mass of the mixture and weighted by the partial mass $\rho_k V_k$ of each species

$$ \bar{u}(\mathbf{r}, t) = \frac{\sum_{k=1}^{N} \rho_k(t) V_k(t) v_k(\mathbf{r}, t)}{\int_{\Omega} \rho(\mathbf{r}, t) dV} \quad (A1) $$

where $\rho_k$ and $V_k$ represent the mass density and partial volume of species $k$, respectively. In the case of a binary A-B mixture, this velocity is expressed as follows:

$$ \bar{u}(\mathbf{r}, t) = \frac{(\rho_A(t)V_A(t)) v_A(\mathbf{r}, t) + (\rho_B(t)V_B(t)) v_B(\mathbf{r}, t)}{\int_{\Omega} \rho(\mathbf{r}, t) dV} \quad (A2) $$
density equations. We explore density ratios ranging from $s = 0.5$ to $s = 8$, while maintaining a Reynolds number of 1000 and a Schmidt number of 1. Initially, it is noted that the microscale Reynolds number $Re_λ$ reaches different values depending on $s$. For a better understanding of variable-density effects, statistical analyses are conducted within a time interval of approximately 2 eddy turnover times. In this interval, $Re_λ$ and the skewness factor remain nearly constant, while the effects of mass diffusion stay significant. In order to mitigate the occurrence of transient effects resulting from the abrupt introduction of the inhomogeneous density field into the turbulent initial condition, a filter is applied during the first hundred iterations of the simulation. Additionally, the statistical analysis interval commences two eddy turnover times later than the start of the simulation.

Special attention is paid to the PDFs of velocity gradients and Lagrangian acceleration. DNS data show that PDFs for variable-density flows exhibit wide tails, which become even wider as $s$ increases. As a result, a noticeable increase in space- and time-intermittency is observed once a threshold $s$ value is exceeded. As quantitative measurements of the deviation from Gaussian distributions, we study the skewness $S$ and flatness factors $F$ of the longitudinal velocity-derivative for various $s$ values. Our findings indicate that the skewness factor increases with the density ratio, even for $Re_λ < 200$. Consequently, the PDFs gradually become more symmetrical as the value of $s$ increases. Additionally, the flatness factor shows an increasing trend with respect to $Re_λ$. These results are consistent with trends observed in constant-density scenarios. However, the $Re_λ$-power-laws found in the literature for homogeneous density fields tend to overestimate both $S$ and $F$ values.

Flow visualization reveals that variable-density flows exhibit smaller coherent structures compared to constant-density flows. Fourier analyses show there is a significant modification in energy levels within the smallest scales of the inertial range. In particular, the region where the energy levels are altered matches area where the Rayleigh-Taylor instability criterion peaks. It suggests that the introduction of variable-density effects into isotropic turbulence results in the generation of smaller eddy structures through Rayleigh-Taylor-like instability.

As a perspective, the logical continuation of this study should be devoted to the analysis of variable-density jets that undergo a transition from laminar to turbulent flow. In particular, density variations are known to shift the jet dynamics though subtle mechanisms. In that respect, the work of Di Pierro and Abid [41], Di Pierro and Abid [37], Ravier et al. [42], and Jacques et al. [40] dealing with inhomogeneous jets can be mentioned. Hence, it would be interesting to evaluate if statistics obtained within the framework of homogeneous isotropic turbulence can be compared to those computed from variable-density turbulent jet simulations. Specifically, one may wonder whether the alteration of energy levels in the small scales observed in the present study — which could be associated to RTI — is observed.

## ACKNOWLEDGMENTS

This work was granted access to the HPC resources of PMCS2I (Pôle de Modélisation et de Calcul en Sciences de l'Ingénieur de l'Information) of École Centrale de Lyon, and P2CHPD (Pôle de Calcul Hautes Performances Dédiés) of Université Lyon1.

## Appendix A: Relation between mean-mass and mean-volume velocity fields

The incompressible formulation serves as an efficient and self-contained tool for describing the physics of low-Mach number flows, whether they exhibit constant or variable-density. This formulation is based on the assumption that velocity field is solenoidal. However, in a mixing flow, there are multiple possible representations of the velocity field. Each species $k$ has a velocity $v_k$, defined as the Eulerian-averaged velocity of the $k$ species molecule. In a DNS approach aimed at resolving the entire flow within the study domain $\Omega$, the goal is to work with an ensemble velocity that encompasses the behavior of all these species. For modeling the ensemble velocity, two approaches are conceivable:

* The mean-mass velocity $\bar{u}$, averaged over the total mass of the mixture and weighted by the partial mass $\rho_k V_k$ of each species

$$ \bar{u}(\mathbf{r}, t) = \frac{\sum_{k=1}^{N} \rho_k(t) V_k(t) v_k(\mathbf{r}, t)}{\int_{\Omega} \rho(\mathbf{r}, t) dV} \quad (A1) $$

where $\rho_k$ and $V_k$ represent the mass density and partial volume of species $k$, respectively. In the case of a binary A-B mixture, this velocity is expressed as follows:

$$ \bar{u}(\mathbf{r}, t) = \frac{(\rho_A(t)V_A(t)) v_A(\mathbf{r}, t) + (\rho_B(t)V_B(t)) v_B(\mathbf{r}, t)}{\int_{\Omega} \rho(\mathbf{r}, t) dV} \quad (A2) $$
In simpler terms, this ensemble velocity field can be referred to as the “barycentric velocity”. This baricentric representation is employed for the momentum balances and consequently appears in the Navier–Stokes equations.

• The mean-volume velocity $\mathbf{u}$, averaged over the total volume of the mixture $V$ and weighted by the volume fraction of each species

$$ \mathbf{u}(\mathbf{r}, t) = \frac{\sum_{k=1}^{N} V_k(t) \mathbf{v}_k(\mathbf{r}, t)}{V} = \sum_{k=1}^{N} \phi_k(t) \mathbf{v}_k(\mathbf{r}, t), \qquad (A3) $$

where $\phi_k$ represent the volumetric fraction of each species. In the case of a binary A–B mixture, this velocity is expressed as follows:

$$ \mathbf{u}(\mathbf{r}, t) = \frac{V_A(t)\mathbf{v}_A(\mathbf{r}, t) + V_B(t)\mathbf{v}_B(\mathbf{r}, t)}{V}. \qquad (A4) $$

In more common terms, this ensemble velocity field can be referred to as “kinematic velocity”.

By definition, a flow is considered incompressible if the elemental volume of a given fluid particle remains constant over time. Mathematically, this is expressed as follows:

$$ \lim_{V \to 0} \frac{1}{V} \frac{dV}{dt} = 0. \qquad (A5) $$

In a variable-density context, the change in the elemental volume $V$ is equal to the sum of the changes in the partial volumes of each species. Using a transport equation for $V_k$ with $\mathbf{v}_k$ within the elemental domain $\Delta$ of the particle, it can be expressed as follows:

$$ \frac{dV}{dt} = \sum_{k=1}^{N} \frac{dV_k}{dt} = \sum_{k=1}^{N} \int_{\Delta} \left( \frac{\partial \phi_k}{\partial t} + \nabla \cdot (\phi_k \mathbf{v}_k) \right) dV = \int_{\Delta} \nabla \cdot \mathbf{u} \, dV. \qquad (A6) $$

The constraint of incompressibility in the case of a mixing flow is thus expressed as follows:

$$ \lim_{V \to 0} \frac{1}{V} \frac{dV}{dt} = \nabla \cdot \mathbf{u} = 0 \qquad (A7) $$

which states that the mean-volume velocity $\mathbf{u}$ must have zero divergence. It is important to note that mean-mass velocity $\bar{\mathbf{u}}$ is not necessarily solenoidal in the case of a variable-density flow. This particularity arises because this representation of the velocity field depends on the mass variation within the elemental volumes associated with each particle. However, these two representations are not independent and can be connected through the mass diffusive flux $\mathbf{J}_k$ here expressed by Fick’s law:

$$ \bar{\mathbf{u}} = \mathbf{u} - \sum_{k=1}^{N} \phi_k (\mathbf{v}_k - \bar{\mathbf{u}}) = \mathbf{u} - \sum_{k=1}^{N} \frac{\phi_k}{\rho_k} \mathbf{J}_k = \mathbf{u} + \rho \sum_{k=1}^{N} \sum_{i \neq k}^{N} D_{ki} \frac{\phi_k}{\rho_k} \nabla w_k \qquad (A8) $$

where $w_k(\mathbf{r}, t) = m_k(\mathbf{r}, t)/M$ represents the local mass fraction of species $k$ with respect to the total mixture mass $M$ at a given position $\mathbf{r}$ and time $t$. In the case of a binary mixture with equal diffusion coefficient, the relation simplifies to

$$ \bar{\mathbf{u}} = \mathbf{u} + \rho D \nabla \left( \frac{\phi_A}{\rho_A} w_A + \frac{\phi_B}{\rho_B} w_B \right) = \mathbf{u} + \rho D \nabla \left( \frac{1}{\rho_A} + \frac{1}{\rho_B} \right) = \mathbf{u} - \frac{D}{\rho} \nabla \rho, \qquad (A9) $$

where $D$ represent the common diffusion coefficient for both species.

The incompressible formulation of the equations of motion of the mixing is derived by combining the Navier–Stokes equations with the incompressibility constraint developed above (eq. (A7)). The Navier–Stokes equations are obtained by applying Newont’s second law to fluid motion and thus use the baricentric velocity $\bar{\mathbf{u}}$. By employing Fick’s law from eq. (A9) to relate the two representations of the velocity field, the equations are rewritten using only the mean-mass velocity. This approach, notably detailed by Guillén-González et al. [28], leads to the system (1) used in this study.
In simpler terms, this ensemble velocity field can be referred to as the “barycentric velocity”. This baricentric representation is employed for the momentum balances and consequently appears in the Navier–Stokes equations.

• The mean-volume velocity $\mathbf{u}$, averaged over the total volume of the mixture $V$ and weighted by the volume fraction of each species

$$ \mathbf{u}(\mathbf{r}, t) = \frac{\sum_{k=1}^{N} V_k(t) \mathbf{v}_k(\mathbf{r}, t)}{V} = \sum_{k=1}^{N} \phi_k(t) \mathbf{v}_k(\mathbf{r}, t), \qquad (A3) $$

where $\phi_k$ represent the volumetric fraction of each species. In the case of a binary A–B mixture, this velocity is expressed as follows:

$$ \mathbf{u}(\mathbf{r}, t) = \frac{V_A(t)\mathbf{v}_A(\mathbf{r}, t) + V_B(t)\mathbf{v}_B(\mathbf{r}, t)}{V}. \qquad (A4) $$

In more common terms, this ensemble velocity field can be referred to as “kinematic velocity”.

By definition, a flow is considered incompressible if the elemental volume of a given fluid particle remains constant over time. Mathematically, this is expressed as follows:

$$ \lim_{V \to 0} \frac{1}{V} \frac{dV}{dt} = 0. \qquad (A5) $$

In a variable-density context, the change in the elemental volume $V$ is equal to the sum of the changes in the partial volumes of each species. Using a transport equation for $V_k$ with $\mathbf{v}_k$ within the elemental domain $\Delta$ of the particle, it can be expressed as follows:

$$ \frac{dV}{dt} = \sum_{k=1}^{N} \frac{dV_k}{dt} = \sum_{k=1}^{N} \int_{\Delta} \left( \frac{\partial \phi_k}{\partial t} + \nabla \cdot (\phi_k \mathbf{v}_k) \right) dV = \int_{\Delta} \nabla \cdot \mathbf{u} \, dV. \qquad (A6) $$

The constraint of incompressibility in the case of a mixing flow is thus expressed as follows:

$$ \lim_{V \to 0} \frac{1}{V} \frac{dV}{dt} = \nabla \cdot \mathbf{u} = 0 \qquad (A7) $$

which states that the mean-volume velocity $\mathbf{u}$ must have zero divergence. It is important to note that mean-mass velocity $\bar{\mathbf{u}}$ is not necessarily solenoidal in the case of a variable-density flow. This particularity arises because this representation of the velocity field depends on the mass variation within the elemental volumes associated with each particle. However, these two representations are not independent and can be connected through the mass diffusive flux $\mathbf{J}_k$ here expressed by Fick’s law:

$$ \bar{\mathbf{u}} = \mathbf{u} - \sum_{k=1}^{N} \phi_k (\mathbf{v}_k - \bar{\mathbf{u}}) = \mathbf{u} - \sum_{k=1}^{N} \frac{\phi_k}{\rho_k} \mathbf{J}_k = \mathbf{u} + \rho \sum_{k=1}^{N} \sum_{i \neq k}^{N} D_{ki} \frac{\phi_k}{\rho_k} \nabla w_k \qquad (A8) $$

where $w_k(\mathbf{r}, t) = m_k(\mathbf{r}, t)/M$ represents the local mass fraction of species $k$ with respect to the total mixture mass $M$ at a given position $\mathbf{r}$ and time $t$. In the case of a binary mixture with equal diffusion coefficient, the relation simplifies to

$$ \bar{\mathbf{u}} = \mathbf{u} + \rho D \nabla \left( \frac{\phi_A}{\rho_A} w_A + \frac{\phi_B}{\rho_B} w_B \right) = \mathbf{u} + \rho D \nabla \left( \frac{1}{\rho_A} + \frac{1}{\rho_B} \right) = \mathbf{u} - \frac{D}{\rho} \nabla \rho, \qquad (A9) $$

where $D$ represent the common diffusion coefficient for both species.

The incompressible formulation of the equations of motion of the mixing is derived by combining the Navier–Stokes equations with the incompressibility constraint developed above (eq. (A7)). The Navier–Stokes equations are obtained by applying Newont’s second law to fluid motion and thus use the baricentric velocity $\bar{\mathbf{u}}$. By employing Fick’s law from eq. (A9) to relate the two representations of the velocity field, the equations are rewritten using only the mean-mass velocity. This approach, notably detailed by Guillén-González et al. [28], leads to the system (1) used in this study.
## Appendix B: Isotropy degree calculation

The HIT modeling requires a condition of statistical isotropy within the considered flow. This implies that the statistical properties of turbulent fluctuations exhibit uniformity in all spatial directions. Departure from isotropy can be quantified through the trace-deviator splitting of the real part of the spectral tensor $\hat{R}_{ij}(k, t)$

$$ \mathrm{Re} \left( \hat{R}_{ij}(k, t) \right) = \underbrace{\frac{E(k, t)}{4\pi k^2} \mathbf{P}_{ij}(k)}_{\text{isotropic part}} + \underbrace{\left( \mathcal{E}(k, t) - \frac{E(k, t)}{4\pi k^2} \right) \mathbf{P}_{ij}(k)}_{\text{directional anisotropy}} + \underbrace{\mathrm{Re} \left( \hat{R}_{ij}^{(\mathrm{pol})}(k, t) \right)}_{\text{polarization anisotropy}}, \quad \mathbf{P}_{ij}(k) = \delta_{ij} - \frac{k_i k_j}{k^2} $$

in terms of isotropic, directional and polarization parts [43–45], where $E(k, t)$ is the energy spectra that depends only on the wave-number modulus $k$ and no longer on the orientation of the vector $k$.

### 1. Polarization anisotropy

Polarization anisotropy is measured by computing the degree of isotropy of the velocity field after each temporal integration using a measurement inspired from the methodology introduced by Curry et al. [33]. The measurement is notably depicted in fig. 3 for simulation 256-4.

This approach involves constructing a vector basis composed of a unit vector $e_z$ of the Cartesian coordinates in physical space, and two spectral unit vectors

$$ e_1(k) = \frac{e_z \times k}{\|e_z \times k\|}, \quad e_2(k) = \frac{e_1 \times k}{\|e_1 \times k\|} $$

In the absence of carrier effects along the $e_z$ direction, nonlinear effects along with mixing lead to an equal distribution of energy of $\hat{u}(k)$ along the $e_1(k)$ and $e_2(k)$ axes [33]. A measure of isotropy is thus provided by comparing the contributions $\psi_1$ and $\psi_2$ of the kinetic energy per unit mass $\frac{1}{2}u$ in the spectral space:

$$ \psi_1 = \langle \|e_1 \cdot \hat{u}(k)\|^2 \rangle, \quad \psi_2 = \langle \|e_2 \cdot \hat{u}(k)\|^2 \rangle $$

where the angular brackets $\langle \cdot \rangle$ represent an average over the entire periodic volume, or in other words, a summation over all available wave numbers. For a fully isotropic flow, $\psi_1 = \psi_2$; thus, an approximate measure of deviation from isotropy is given by:

$$ I^2 = \psi_1/\psi_2. $$

A flow is considered close to a fully isotropic state when $0.95 \le I \le 1.05$.

It is crucial to note that when a nonisotropic forcing method is implemented to maintain the HIT, the modes associated with energy injection should not be considered in the isotropy measurement. In this study, simulations are initialized similarly to the method of [8], which excites the first modes in each direction (see eq. (12)). This initial condition is not isotropic. To eliminate a signature in the isotropy measurement, the first and last modes of each direction are not counted. The exclusion of these modes does not pose an issue concerning the HIT modeling constraint since it requires isotropy of fluctuations. The first modes are associated with the large flow structures, while the fluctuations at other scales are independent of the behavior of the large vortices, in accordance with Kolmogorov’s theory [7].

### 2. Directional anisotropy

Directional anisotropy can be measured by computing the angular distribution of $\mathcal{E} - E/(4\pi k^2)$. By using the spherical polar coordinate system $\{e_1, e_2, e_z\}$ in $k$-space defined earlier, the spectral Reynolds tensor and $\mathcal{E}(k, t)$ simplifies [43] as

$$ [\hat{R}_{ij}] = \begin{bmatrix} \phi_{11} & \phi_{12} & 0 \\ \phi_{12}^* & \phi_{22} & 0 \\ 0 & 0 & 0 \end{bmatrix}, \quad \mathcal{E}(k, t) = \frac{1}{2} \hat{R}_{ii} = \frac{1}{2} (\phi_{11} + \phi_{22}). $$
Appendix B: Isotropy degree calculation

The HIT modeling requires a condition of statistical isotropy within the considered flow. This implies that the statistical properties of turbulent fluctuations exhibit uniformity in all spatial directions. Departure from isotropy can be quantified through the trace-deviator splitting of the real part of the spectral tensor $\hat{R}_{ij}(k, t)$

$$ \mathrm{Re} \left( \hat{R}_{ij}(k, t) \right) = \underbrace{\frac{E(k, t)}{4\pi k^2} \mathbf{P}_{ij}(k)}_{\text{isotropic part}} + \underbrace{\left( \mathcal{E}(k, t) - \frac{E(k, t)}{4\pi k^2} \right) \mathbf{P}_{ij}(k)}_{\text{directional anisotropy}} + \underbrace{\mathrm{Re} \left( \hat{R}_{ij}^{(\mathrm{pol})}(k, t) \right)}_{\text{polarization anisotropy}}, \quad \mathbf{P}_{ij}(k) = \delta_{ij} - \frac{k_i k_j}{k^2} $$

in terms of isotropic, directional and polarization parts [43–45], where $E(k, t)$ is the energy spectra that depends only on the wave-number modulus $k$ and no longer on the orientation of the vector $k$.

1. Polarization anisotropy

Polarization anisotropy is measured by computing the degree of isotropy of the velocity field after each temporal integration using a measurement inspired from the methodology introduced by Curry et al. [33]. The measurement is notably depicted in fig. 3 for simulation 256-4.

This approach involves constructing a vector basis composed of a unit vector $e_z$ of the Cartesian coordinates in physical space, and two spectral unit vectors

$$ e_1(k) = \frac{e_z \times k}{\|e_z \times k\|}, \quad e_2(k) = \frac{e_1 \times k}{\|e_1 \times k\|} $$

In the absence of carrier effects along the $e_z$ direction, nonlinear effects along with mixing lead to an equal distribution of energy of $\hat{u}(k)$ along the $e_1(k)$ and $e_2(k)$ axes [33]. A measure of isotropy is thus provided by comparing the contributions $\psi_1$ and $\psi_2$ of the kinetic energy per unit mass $\frac{1}{2}u$ in the spectral space:

$$ \psi_1 = \langle \|e_1 \cdot \hat{u}(k)\|^2 \rangle, \quad \psi_2 = \langle \|e_2 \cdot \hat{u}(k)\|^2 \rangle $$

where the angular brackets $\langle \cdot \rangle$ represent an average over the entire periodic volume, or in other words, a summation over all available wave numbers. For a fully isotropic flow, $\psi_1 = \psi_2$; thus, an approximate measure of deviation from isotropy is given by:

$$ I^2 = \psi_1/\psi_2. $$

A flow is considered close to a fully isotropic state when $0.95 \le I \le 1.05$.

It is crucial to note that when a nonisotropic forcing method is implemented to maintain the HIT, the modes associated with energy injection should not be considered in the isotropy measurement. In this study, simulations are initialized similarly to the method of [8], which excites the first modes in each direction (see eq. (12)). This initial condition is not isotropic. To eliminate a signature in the isotropy measurement, the first and last modes of each direction are not counted. The exclusion of these modes does not pose an issue concerning the HIT modeling constraint since it requires isotropy of fluctuations. The first modes are associated with the large flow structures, while the fluctuations at other scales are independent of the behavior of the large vortices, in accordance with Kolmogorov’s theory [7].

2. Directional anisotropy

Directional anisotropy can be measured by computing the angular distribution of $\mathcal{E} - E/(4\pi k^2)$. By using the spherical polar coordinate system $\{e_1, e_2, e_z\}$ in $k$-space defined earlier, the spectral Reynolds tensor and $\mathcal{E}(k, t)$ simplifies [43] as

$$ [\hat{R}_{ij}] = \begin{bmatrix} \phi_{11} & \phi_{12} & 0 \\ \phi_{12}^* & \phi_{22} & 0 \\ 0 & 0 & 0 \end{bmatrix}, \quad \mathcal{E}(k, t) = \frac{1}{2} \hat{R}_{ii} = \frac{1}{2} (\phi_{11} + \phi_{22}). $$
<figure><img src="image_17.png" /><figcaption>FIG. 17: Evolution of the relative difference between $\mathcal{E}(k)$ and $E(k, t)/4\pi k^2$ during simulation run 256-8. The difference between these two terms measures directional anisotropy. Here, a difference of a few percents is observed, indicating that the flow is isotropic. A similar behavior is obtained in all simulation runs.</figcaption></figure>

From this, one can compute the angular distribution $\mathcal{E}-E/(4\pi k^2)$ in order to verify that this directional anisotropy term is indeed nearly zero. Figure 17 shows the evolution of the maximum relative difference between $\mathcal{E}$ and $E/(4\pi k^2)$ during simulation run 256-8. Here too, the modes associated with energy injection from forcing are excluded. It is observed that the relative difference ranges between 0.1 and 5% which confirms that the departure from isotropy due to directional anisotropy is negligible. Similar results are observed for all other simulations runs, with a maximum relative difference ranging from 5-6%.

## Appendix C: Evolution of turbulence characteristics

Figures 18 to 21 present the time evolution of various turbulence characteristic during the variable density HIT DNSs. Simulation parameters for variable-density HIT DNSs are listed in table I page 5. Time evolutions for run 256-4 are presented in the main text in fig. 3 page 8.

[1] D. Livescu and J. R. Ristorcelli, Buoyancy-driven variable-density turbulence, Journal of Fluid Mechanics **591**, 43 (2007).

[2] D. Livescu, J. R. Ristorcelli, M. R. Petersen, and R. A. Gore, New phenomena in variable-density Rayleigh-Taylor turbulence, Physica Scripta **T142**, 014015 (2010).

[3] P. Rao, C.-C. P. Caulfield, and J. D. Gibbon, Nonlinear effects in buoyancy-driven variable-density turbulence, Journal of Fluid Mechanics **810**, 362 (2017).

[4] K. K. Nomura and S. E. Elghobashi, The structure of inhomogeneous turbulence in variable density nonpremixed flames, Theoretical and Computational Fluid Dynamics **5-5**, 153 (1993).

[5] M. Mohaghar, J. Carter, B. Musci, D. Reilly, J. McFarland, and D. Ranjan, Evaluation of turbulent mixing transition in a shock-driven variable-density flow, Journal of Fluid Mechanics **831**, 779 (2017).

[6] Y. Tian, F. A. Jaberi, and D. Livescu, Density effects on post-shock turbulence structure and dynamics, Journal of Fluid Mechanics **880**, 935 (2019).

[7] A. N. Kolmogorov, The local structure of turbulence in incompressible viscous fluid for very large Reynolds numbers, Proceedings of the Royal Society of London. Series A: Mathematical and Physical Sciences **434**, 9 (1991).

[8] R. M. Kerr, Higher-order derivative correlations and the alignment of small-scale structures in isotropic numerical turbulence, Journal of Fluid Mechanics **153**, 31 (1985).

[9] A. Vincent and M. Meneguzzi, The spatial structure and statistical properties of homogeneous turbulence, Journal of Fluid Mechanics **225**, 1 (1991).

[10] J. Jiménez, A. A. Wray, P. G. Saffman, and R. S. Rogallo, The structure of intense vorticity in isotropic turbulence, Journal of Fluid Mechanics **255**, 65 (1993).

[11] P. Moin and K. Mahesh, Direct Numerical Simulation: A tool in turbulence research, Annual Review of Fluid Mechanics **30**, 539 (1998).
<figure><img src="image_18.png" /><figcaption>FIG. 17: Evolution of the relative difference between $\mathcal{E}(k)$ and $E(k, t)/4\pi k^2$ during simulation run 256-8. The difference between these two terms measures directional anisotropy. Here, a difference of a few percents is observed, indicating that the flow is isotropic. A similar behavior is obtained in all simulation runs.</figcaption></figure>

From this, one can compute the angular distribution $\mathcal{E}-E/(4\pi k^2)$ in order to verify that this directional anisotropy term is indeed nearly zero. Figure 17 shows the evolution of the maximum relative difference between $\mathcal{E}$ and $E/(4\pi k^2)$ during simulation run 256-8. Here too, the modes associated with energy injection from forcing are excluded. It is observed that the relative difference ranges between 0.1 and 5% which confirms that the departure from isotropy due to directional anisotropy is negligible. Similar results are observed for all other simulations runs, with a maximum relative difference ranging from 5-6%.

## Appendix C: Evolution of turbulence characteristics

Figures 18 to 21 present the time evolution of various turbulence characteristic during the variable density HIT DNSs. Simulation parameters for variable-density HIT DNSs are listed in table I page 5. Time evolutions for run 256-4 are presented in the main text in fig. 3 page 8.

[1] D. Livescu and J. R. Ristorcelli, Buoyancy-driven variable-density turbulence, Journal of Fluid Mechanics **591**, 43 (2007).

[2] D. Livescu, J. R. Ristorcelli, M. R. Petersen, and R. A. Gore, New phenomena in variable-density Rayleigh-Taylor turbulence, Physica Scripta **T142**, 014015 (2010).

[3] P. Rao, C.-C. P. Caulfield, and J. D. Gibbon, Nonlinear effects in buoyancy-driven variable-density turbulence, Journal of Fluid Mechanics **810**, 362 (2017).

[4] K. K. Nomura and S. E. Elghobashi, The structure of inhomogeneous turbulence in variable density nonpremixed flames, Theoretical and Computational Fluid Dynamics **5-5**, 153 (1993).

[5] M. Mohaghar, J. Carter, B. Musci, D. Reilly, J. McFarland, and D. Ranjan, Evaluation of turbulent mixing transition in a shock-driven variable-density flow, Journal of Fluid Mechanics **831**, 779 (2017).

[6] Y. Tian, F. A. Jaberi, and D. Livescu, Density effects on post-shock turbulence structure and dynamics, Journal of Fluid Mechanics **880**, 935 (2019).

[7] A. N. Kolmogorov, The local structure of turbulence in incompressible viscous fluid for very large Reynolds numbers, Proceedings of the Royal Society of London. Series A: Mathematical and Physical Sciences **434**, 9 (1991).

[8] R. M. Kerr, Higher-order derivative correlations and the alignment of small-scale structures in isotropic numerical turbulence, Journal of Fluid Mechanics **153**, 31 (1985).

[9] A. Vincent and M. Meneguzzi, The spatial structure and statistical properties of homogeneous turbulence, Journal of Fluid Mechanics **225**, 1 (1991).

[10] J. Jiménez, A. A. Wray, P. G. Saffman, and R. S. Rogallo, The structure of intense vorticity in isotropic turbulence, Journal of Fluid Mechanics **255**, 65 (1993).

[11] P. Moin and K. Mahesh, Direct Numerical Simulation: A tool in turbulence research, Annual Review of Fluid Mechanics **30**, 539 (1998).
<figure><img src="image_18.png" /><figcaption>FIG. 18: Evolution of the volumetric total kinetic energy, isotropy degree [33] (calculation detailed in appendix B), microscale Reynolds number $Re_λ$, density field standard deviation in space $σ(ρ)$, velocity longitudinal-derivative skewness $S$, and eddy turnover time $T$ during simulation run 256-2. Vertical lines delimit the statistical analysis time interval $Γ$. Simulation parameters for variable-density HIT DNSs are listed in table I page 5.</figcaption></figure>
<figure><img src="image_19.png" /><figcaption>FIG. 18: Evolution of the volumetric total kinetic energy, isotropy degree [33] (calculation detailed in appendix B), microscale Reynolds number $Re_λ$, density field standard deviation in space $σ(ρ)$, velocity longitudinal-derivative skewness $S$, and eddy turnover time $T$ during simulation run 256-2. Vertical lines delimit the statistical analysis time interval $Γ$. Simulation parameters for variable-density HIT DNSs are listed in table I page 5.</figcaption></figure>
<figure><img src="image_19.png" /><figcaption>FIG. 19: Evolution of the volumetric total kinetic energy, isotropy degree [33] (calculation detailed in appendix B), microscale Reynolds number $Re_λ$, density field standard deviation in space $σ(ρ)$, velocity longitudinal-derivative skewness $S$, and eddy turnover time $T$ during simulation run 256-6. Vertical lines delimit the statistical analysis time interval $Γ$. Simulation parameters for variable-density HIT DNSs are listed in table I page 5.</figcaption></figure>
<figure><img src="image_20.png" /><figcaption>FIG. 19: Evolution of the volumetric total kinetic energy, isotropy degree [33] (calculation detailed in appendix B), microscale Reynolds number $Re_λ$, density field standard deviation in space $σ(ρ)$, velocity longitudinal-derivative skewness $S$, and eddy turnover time $T$ during simulation run 256-6. Vertical lines delimit the statistical analysis time interval $Γ$. Simulation parameters for variable-density HIT DNSs are listed in table I page 5.</figcaption></figure>
<figure><img src="image_20.png" /><figcaption>FIG. 20: Evolution of the volumetric total kinetic energy, isotropy degree [33] (calculation detailed in appendix B), microscale Reynolds number $Re_λ$, density field standard deviation in space $σ(ρ)$, velocity longitudinal-derivative skewness $S$, and eddy turnover time $T$ during simulation run 256-8. Vertical lines delimit the statistical analysis time interval $Γ$. Simulation parameters for variable-density HIT DNSs are listed in table I page 5.</figcaption></figure>
<figure><img src="image_21.png" /><figcaption>FIG. 20: Evolution of the volumetric total kinetic energy, isotropy degree [33] (calculation detailed in appendix B), microscale Reynolds number $Re_λ$, density field standard deviation in space $σ(ρ)$, velocity longitudinal-derivative skewness $S$, and eddy turnover time $T$ during simulation run 256-8. Vertical lines delimit the statistical analysis time interval $Γ$. Simulation parameters for variable-density HIT DNSs are listed in table I page 5.</figcaption></figure>
<figure><img src="image_21.png" /><figcaption>FIG. 21: Evolution of the volumetric total kinetic energy, isotropy degree [33] (calculation detailed in appendix B), microscale Reynolds number $Re_λ$, density field standard deviation in space $σ(ρ)$, velocity longitudinal-derivative skewness $S$, and eddy turnover time $T$ during simulation run 256-0.5. Vertical lines delimit the statistical analysis time interval $Γ$. Simulation parameters for variable-density HIT DNSs are listed in table I page 5.</figcaption></figure>
<figure><img src="image_22.png" /><figcaption>FIG. 21: Evolution of the volumetric total kinetic energy, isotropy degree [33] (calculation detailed in appendix B), microscale Reynolds number $Re_λ$, density field standard deviation in space $σ(ρ)$, velocity longitudinal-derivative skewness $S$, and eddy turnover time $T$ during simulation run 256-0.5. Vertical lines delimit the statistical analysis time interval $Γ$. Simulation parameters for variable-density HIT DNSs are listed in table I page 5.</figcaption></figure>
[12] Z. Qin, L. Fang, and J. Fang, How isotropic are turbulent flows generated by using periodic conditions in a cube?, Physics Letters A **380**, 1310 (2016).

[13] T. Ishihara, Y. Kaneda, M. Yokokawa, K. Itakura, and A. Uno, Small-scale statistics in high-resolution direct numerical simulation of turbulence: Reynolds number dependence of one-point velocity gradient statistics, Journal of Fluid Mechanics **592**, 335 (2007).

[14] K. Sreenivasan and R. A. Antonia, The phenomenology of small-scale turbulence., Annu. Rev. Fluid Mech. **29**, 435 (1997).

[15] B. Dubrulle, Beyond Kolmogorov cascades, Journal of Fluid Mechanics **867**, 1 (2019).

[16] E. D. Siggia and G. S. Patterson, Intermittency effects in a numerical simulation of stationary three-dimensional turbulence, Journal of Fluid Mechanics **86**, 567 (1978).

[17] M. R. Overholt and S. B. Pope, A deterministic forcing scheme for direct numerical simulations of turbulence, Computers & Fluids **27**, 11 (1998).

[18] V. Eswaran and S. B. Pope, An examination of forcing in direct numerical simulations of turbulence, Computers & Fluids **16**, 257 (1988).

[19] T. S. Lundgren, Linearly forced isotropic turbulence, in *Annual Research Briefs* (Center for Turbulence Research, Stanford, 2003) pp. 461–473.

[20] P. L. Carroll and G. Blanquart, A proposed modification to Lundgren’s physical space velocity forcing method for isotropic turbulence, Physics of Fluids **25**, 105114 (2013).

[21] C. Rosales and C. Meneveau, Linear forcing in numerical simulations of isotropic turbulence: Physical space implementations and convergence properties, Physics of Fluids **17**, 095106 (2005).

[22] J. Janin, F. Duval, C. Friess, and P. Sagaut, A new linear forcing method for isotropic turbulence with controlled integral length scale, Physics of Fluids **33**, 045127 (2021).

[23] O. Desjardins, G. Blanquart, G. Balarac, and H. Pitsch, High order conservative finite difference scheme for variable density low mach number turbulent flows, Journal of Computational Physics **227**, 7125 (2008).

[24] L. Reynier, B. Di Pierro, F. Alizard, A. Cadiou, L. Le Penven, and M. Buffat, A preconditioning for the spectral solution of incompressible variable-density flows, Computers & Fluids **266**, 106024 (2023).

[25] D. A. Frank-Kamenetskii, *Diffusion and Heat Exchange in Chemical Kinetics*, Vol. 2171 (Princeton University Press, 1955).

[26] A. V. Kazhikhov and S. S. Smagulov, The correctness of boundary-value problems in a diffusion model of an inhomogeneous liquid, Dokl. Akad. Nauk SSSR **234**, 330 (1977).

[27] S. Antontsev, A. Kazhiktov, and V. Monakhov, *Boundary Value Problems in Mechanics of Nonhomogeneous Fluids* (Elsevier, North Holland Amsterdam, 1990).

[28] F. Guillén-González, P. D. Damázio, and M. A. Rojas-Medar, Approximation by an iterative method for regular solutions for incompressible fluids with mass diffusion, Journal of Mathematical Analysis and Applications **326**, 468 (2007).

[29] R. Peyret, *Spectral Methods for Incompressible Viscous Flow*, Applied mathematical sciences No. 148 (Springer New York, New York, 2002).

[30] J. B. Bell and D. L. Marcus, A second-order projection method for variable-density flows, Journal of Computational Physics **101**, 334 (1992).

[31] C. Tadjeran, Stability analysis of the Crank–Nicholson method for variable coefficient diffusion equation, Communications in Numerical Methods in Engineering **23**, 29 (2006).

[32] B. Di Pierro and M. Abid, A projection method for the spectral solution of non-homogeneous and incompressible Navier–Stokes equations, International Journal for Numerical Methods in Fluids **71**, 1029 (2012).

[33] J. H. Curry, J. R. Herring, J. Loncaric, and S. A. Orszag, Order and disorder in two- and three-dimensional Bénard convection, Journal of Fluid Mechanics **147**, 1 (1984).

[34] T. Gotoh, D. Fukayama, and T. Nakano, Velocity field statistics in homogeneous steady turbulence obtained using a high-resolution direct numerical simulation, Physics of Fluids **14**, 1065 (2002).

[35] A. Gylfason, S. Ayyalasomayajula, and Z. Warhaft, Intermittency, pressure and acceleration statistics from hot-wire measurements in wind-tunnel turbulence, Journal of Fluid Mechanics **501**, 213 (2004).

[36] M. M. Scase and R. J. A. Hill, Centrifugally forced Rayleigh–Taylor instability, Journal of Fluid Mechanics **852**, 543 (2018).

[37] B. Di Pierro and M. Abid, Rayleigh-Taylor instability in variable density swirling flows, The European Physical Journal B **85**, 69 (2012).

[38] C. Jacques, B. Di Pierro, and M. Buffat, Stability of variable density jets (2023), submitted.

[39] R. H. Zeng, J. J. Tao, and Y. B. Sun, Three-dimensional viscous Rayleigh-Taylor instability at the cylindrical interface, Physical Review E **102**, 023112 (2020).

[40] C. Jacques, B. Di Pierro, F. Alizard, M. Buffat, A. Cadiou, and L. Le Penven, Stability of variable density rotating flows: Inviscid case and viscous effects in the limit of large Reynolds numbers, Physical Review Fluids **8**, 033901 (2023).

[41] B. Di Pierro and M. Abid, Instabilities of variable-density swirling flows, Physical Review E **82**, 046312 (2010).

[42] S. Ravier, M. Abid, M. Amielh, and F. Anselmet, Direct numerical simulations of variable-density plane jets, Journal of Fluid Mechanics **546**, 153 (2005).

[43] P. Sagaut and C. Cambon, *Homogeneous turbulence dynamics*, 2nd ed. (Springer, 2008).

[44] V. Mons, C. Cambon, and P. Sagaut, A spectral model for homogeneous shear-driven anisotropic turbulence in terms of spherically averaged descriptors, Journal of Fluid Mechanics **788**, 147 (2015).

[45] A. Briard, M. Iyer, and T. Gomez, Anisotropic spectral modeling for unstably stratified homogeneous turbulence, Physical Review Fluids **2**, 10.1103/physrevfluids.2.044604 (2017).
[12] Z. Qin, L. Fang, and J. Fang, How isotropic are turbulent flows generated by using periodic conditions in a cube?, Physics Letters A **380**, 1310 (2016).

[13] T. Ishihara, Y. Kaneda, M. Yokokawa, K. Itakura, and A. Uno, Small-scale statistics in high-resolution direct numerical simulation of turbulence: Reynolds number dependence of one-point velocity gradient statistics, Journal of Fluid Mechanics **592**, 335 (2007).

[14] K. Sreenivasan and R. A. Antonia, The phenomenology of small-scale turbulence., Annu. Rev. Fluid Mech. **29**, 435 (1997).

[15] B. Dubrulle, Beyond Kolmogorov cascades, Journal of Fluid Mechanics **867**, 1 (2019).

[16] E. D. Siggia and G. S. Patterson, Intermittency effects in a numerical simulation of stationary three-dimensional turbulence, Journal of Fluid Mechanics **86**, 567 (1978).

[17] M. R. Overholt and S. B. Pope, A deterministic forcing scheme for direct numerical simulations of turbulence, Computers & Fluids **27**, 11 (1998).

[18] V. Eswaran and S. B. Pope, An examination of forcing in direct numerical simulations of turbulence, Computers & Fluids **16**, 257 (1988).

[19] T. S. Lundgren, Linearly forced isotropic turbulence, in *Annual Research Briefs* (Center for Turbulence Research, Stanford, 2003) pp. 461–473.

[20] P. L. Carroll and G. Blanquart, A proposed modification to Lundgren’s physical space velocity forcing method for isotropic turbulence, Physics of Fluids **25**, 105114 (2013).

[21] C. Rosales and C. Meneveau, Linear forcing in numerical simulations of isotropic turbulence: Physical space implementations and convergence properties, Physics of Fluids **17**, 095106 (2005).

[22] J. Janin, F. Duval, C. Friess, and P. Sagaut, A new linear forcing method for isotropic turbulence with controlled integral length scale, Physics of Fluids **33**, 045127 (2021).

[23] O. Desjardins, G. Blanquart, G. Balarac, and H. Pitsch, High order conservative finite difference scheme for variable density low mach number turbulent flows, Journal of Computational Physics **227**, 7125 (2008).

[24] L. Reynier, B. Di Pierro, F. Alizard, A. Cadiou, L. Le Penven, and M. Buffat, A preconditioning for the spectral solution of incompressible variable-density flows, Computers & Fluids **266**, 106024 (2023).

[25] D. A. Frank-Kamenetskii, *Diffusion and Heat Exchange in Chemical Kinetics*, Vol. 2171 (Princeton University Press, 1955).

[26] A. V. Kazhikhov and S. S. Smagulov, The correctness of boundary-value problems in a diffusion model of an inhomogeneous liquid, Dokl. Akad. Nauk SSSR **234**, 330 (1977).

[27] S. Antontsev, A. Kazhiktov, and V. Monakhov, *Boundary Value Problems in Mechanics of Nonhomogeneous Fluids* (Elsevier, North Holland Amsterdam, 1990).

[28] F. Guillén-González, P. D. Damázio, and M. A. Rojas-Medar, Approximation by an iterative method for regular solutions for incompressible fluids with mass diffusion, Journal of Mathematical Analysis and Applications **326**, 468 (2007).

[29] R. Peyret, *Spectral Methods for Incompressible Viscous Flow*, Applied mathematical sciences No. 148 (Springer New York, New York, 2002).

[30] J. B. Bell and D. L. Marcus, A second-order projection method for variable-density flows, Journal of Computational Physics **101**, 334 (1992).

[31] C. Tadjeran, Stability analysis of the Crank–Nicholson method for variable coefficient diffusion equation, Communications in Numerical Methods in Engineering **23**, 29 (2006).

[32] B. Di Pierro and M. Abid, A projection method for the spectral solution of non-homogeneous and incompressible Navier–Stokes equations, International Journal for Numerical Methods in Fluids **71**, 1029 (2012).

[33] J. H. Curry, J. R. Herring, J. Loncaric, and S. A. Orszag, Order and disorder in two- and three-dimensional Bénard convection, Journal of Fluid Mechanics **147**, 1 (1984).

[34] T. Gotoh, D. Fukayama, and T. Nakano, Velocity field statistics in homogeneous steady turbulence obtained using a high-resolution direct numerical simulation, Physics of Fluids **14**, 1065 (2002).

[35] A. Gylfason, S. Ayyalasomayajula, and Z. Warhaft, Intermittency, pressure and acceleration statistics from hot-wire measurements in wind-tunnel turbulence, Journal of Fluid Mechanics **501**, 213 (2004).

[36] M. M. Scase and R. J. A. Hill, Centrifugally forced Rayleigh–Taylor instability, Journal of Fluid Mechanics **852**, 543 (2018).

[37] B. Di Pierro and M. Abid, Rayleigh-Taylor instability in variable density swirling flows, The European Physical Journal B **85**, 69 (2012).

[38] C. Jacques, B. Di Pierro, and M. Buffat, Stability of variable density jets (2023), submitted.

[39] R. H. Zeng, J. J. Tao, and Y. B. Sun, Three-dimensional viscous Rayleigh-Taylor instability at the cylindrical interface, Physical Review E **102**, 023112 (2020).

[40] C. Jacques, B. Di Pierro, F. Alizard, M. Buffat, A. Cadiou, and L. Le Penven, Stability of variable density rotating flows: Inviscid case and viscous effects in the limit of large Reynolds numbers, Physical Review Fluids **8**, 033901 (2023).

[41] B. Di Pierro and M. Abid, Instabilities of variable-density swirling flows, Physical Review E **82**, 046312 (2010).

[42] S. Ravier, M. Abid, M. Amielh, and F. Anselmet, Direct numerical simulations of variable-density plane jets, Journal of Fluid Mechanics **546**, 153 (2005).

[43] P. Sagaut and C. Cambon, *Homogeneous turbulence dynamics*, 2nd ed. (Springer, 2008).

[44] V. Mons, C. Cambon, and P. Sagaut, A spectral model for homogeneous shear-driven anisotropic turbulence in terms of spherically averaged descriptors, Journal of Fluid Mechanics **788**, 147 (2015).

[45] A. Briard, M. Iyer, and T. Gomez, Anisotropic spectral modeling for unstably stratified homogeneous turbulence, Physical Review Fluids **2**, 10.1103/physrevfluids.2.044604 (2017).