Variable-density effects in incompressible non-buoyant shear-driven turbulent mixing layers

The asymmetries that arise when a mixing layer involves two miscible fluids of differing densities are investigated using incompressible (low-speed) direct numerical simulations. The simulations are performed in the temporal configuration with very large domain sizes, to allow the mixing layers to reach prolonged states of fully-turbulent self-similar growth. Imposing a mean density variation breaks the mean symmetry relative to the classical single-fluid temporal mixing layer problem. Unlike prior variable-density mixing layer simulations in which the streams are composed of the same fluids with dissimilar thermodynamic properties, the density variations are presently due to compositional differences between the fluid streams, leading to different mixing dynamics. Variable-density (non-Boussinesq) effects introduce strong asymmetries in the flow statistics that can be explained by the strongest turbulence increasingly migrating to the lighter fluid side as free stream density difference increases. Interface thickness growth rates also reduce, with some thickness definitions particularly sensitive to the corresponding changes in alignment between density and streamwise velocity profiles. Additional asymmetries in the sense of statistical distributions of densities at a given position within the mixing layer reveal that fine scales of turbulence are preferentially sustained in lighter fluid, which also is where fastest mixing occurs. These effects influence statistics involving density fluctuations, which have important implications for mixing and more complicated phenomena that are sensitive to the mixing dynamics, such as combustion.


Introduction
A wide range of applications include the fundamental phenomenon of turbulence sustained by shear between streams of fluids. Frequently, the streams may have different densities because they consist of different fluids. Such flows can involve miscible or immiscible fluids; we are here concerned only with the miscible case. Miscible applications exist in combustion, industrial chemical mixing, and geophysical flows. The relevance of mixing layer simulations to combustion is reviewed in Givi (1989), and other complex applications of sheared variable-density flows are summarized in Akula et al. (2013).
In many cases, the density differences can be large, producing significant changes to the flow evolution. Dimotakis (2005), in a review of turbulent mixing, classified mixing into three categories based on the complexity (physics coupling) of the mixing phenomena and the importance of correctly capturing the mixing dynamics to the overall predictions. In the simplest (Level-1) cases, capturing the turbulence but not the mixing itself is sufficient to predict the flow dynamics. Level-2 indicates that mixing alters the flow dynamics. Inertial effects of the large density variations of the mixing layers investigated herein place the flow in Level-2 with increased complexity that cannot be captured by extending single-density mixing layer results with passive mixing.
In combustion, very large density variations can exist due to differing fluid compositions and thermodynamic variations. Combustion is among the most complex mixing flows (classified as Level-3) because the mixing strongly affects reactions that produce changes in the fluids (including heat release) which then couple back to the mixing dynamics. Capturing the inertial effects associated with compositional variations during the mixing of reactants and reaction products can be a significant component of predicting combustion. Bilger (1976) noted the importance of density differences in turbulent jet diffusion flames. In configurations such as a jet of hydrogen fuel released into air, the density differences can be very large simply due to the different molar masses of the fluids.
Several recent incompressible studies have revealed interesting effects on turbulent mixing when density differences are large solely due to differing compositions. The Atwood number A characterizes the difference in densities between streams of fluids: where ρ 1 and ρ 2 are the densities of each pure fluid. Pure helium mixing with air (or nitrogen) corresponds to an Atwood number of 0.75, while pure hydrogen mixing with air corresponds to A = 0.85. Studying Rayleigh-Taylor (RT) instability in the classical configuration and a triply-periodic version (i.e. homogeneous buoyancy-driven turbulence), Livescu & Ristorcelli (2008, 2009) found significant changes in behavior when Atwood number was increased to high values. Atwood numbers of A 0.05 are typically considered to be the limit of the Boussinesq approximation . Flows of sufficiently high Atwood number to vary significantly from the Boussinesq approximation have been termed variable-density. Livescu et al. (2010) showed that changes in alignment between density gradient and local strain is a variable-density effect associated with reduced mixing in the heavy fluid regions. Much of the simulation studies of density effects on mixing have occurred in buoyancy-driven turbulence, such as the small density variation study of Batchelor et al. (1992) that was later extended to non-Boussinesq flow by Sandoval (1995). Sandoval (1995) also considered decaying isotropic turbulence without buoyancy, which was further studied by Jang & de Bruyn Kops (2007). Movahed & Johnsen (2015) studied variable-density mixing in two fluids with decaying isotropic turbulence initially separated by a planar interface. Notable classical RT studies include Cabot & Cook (2006) and . Shear-driven mixing layers have historically received a great deal of attention, but mainly for single-fluid configurations. Rogers & Moser (1994) simulated an incompressible mixing layer in the temporal (streamwise-periodic) configuration to self-similar fullyturbulent growth. A similar configuration was simulated by Balaras et al. (2001) to study the effects of initial conditions. More powerful computational resources have recently enabled performing spatially-developing simulations, which more closely approximate mixing layer experiments. These require much longer streamwise domains to attain a desired mixing layer thickness since they thicken with downstream distance rather than in time as is the case for temporal simulations. (However, meaningful temporal simulations implicitly require sufficiently large domains to not interfere with the growth of turbulent structures.) Wang et al. (2007) designed a spatially-developing mixing layer simulation to be comparable to the temporal mixing layer of Tanahashi et al. (2001) and observed similar energy dissipation rates but increased turbulent kinetic energy. The DNS of Attili & Bisetti (2012) advanced spatially-developing mixing layer simulations to a very long domain that enabled attaining a relatively large Reynolds number. During self-similar growth, they found remarkable agreement between their self-similar dissipation values and that of the Rogers & Moser (1994) temporal simulation, as well as close agreement for most other statistics. Relevant low-speed experimental studies include those of Spencer & Jones (1971), Bell & Mehta (1990) and Loucks & Wallace (2012). Experiments addressing detailed turbulent structure include those of Olsen & Dutton (2003) (which also contained a weak density difference) and Li et al. (2010). In several studies, mixing properties have been investigated with shear-driven mixing layers, but in the absence of density differences between the participating fluids (e.g., Sharan et al. 2019).
High-speed compressible mixing layers have also received a great deal of attention, particularly due to the strong reduction in mixing layer growth rate that occurs with increasing Mach number. Though density effects associated with compressibility were once thought to affect growth rate (as discussed in Brown & Roshko 1974), DNS simulations have clarified how compressibility effects reduce the growth due to decreased turbulent kinetic energy production as compressibility decorrelates the strain and pressure fluctuations (Vreman et al. 1996;Sarkar 1996;Freund et al. 2000;Pantano & Sarkar 2002;Livescu & Madnia 2004). Research has continued on this mechanism in compressible mixing layer experiments (e.g., Barre & Bonnet 2015). Recent simulations have further investigated the mixing characteristics of compressible mixing layers (e.g., Jahanbakhshi et al. 2015).
Non-buoyant mixing layers with significant density variations (i.e. density ratios larger than 2) have begun to receive attention. 2D and 3D simulations demonstrated that differing free-stream densities significantly changed the early-time growth and Kelvin-Helmholtz (KH) flow structures (Joly et al. 2001;Joly 2002). The pioneering 3D temporal simulations of Pantano & Sarkar (2002) included an investigation of different freestream densities within a broader study of compressible mixing layers. The differing densities were established by varying the temperature for a single fluid. They found that increasing Atwood number decreased the temporal thickness growth rate, though the extent depended on how thickness was defined. During self-similar growth, the Reynolds shear stress changed little in magnitude but shifted to the light fluid side with increasing Atwood number. They also developed a model characterizing the shift of the mean velocity profile to the light fluid side and the associated decrease in momentum thickness growth rate. Mild compressibility effects were likely present because the convective Mach number was M c = 0.7. More recently, Almagro et al. (2017) performed DNS using a lowspeed approximation for the flow of Pantano & Sarkar (2002). Two streams of a single fluid with different temperatures again create the density difference, but compressibility effects are considered negligible at low speeds. They also developed a semi-empirical model for the reduction in momentum thickness growth rate with density ratio.
Details of mixing layers with variable density due to differing fluid compositions are much less understood. Detailed studies of mixing layers involving two different miscible fluids have been rare, particularly when not complicated by other effects such as buoyancy or compressibility, despite earlier attention. The historic low-speed experiments of Brown & Roshko (1974) using two gases with different densities found reductions in the growth rates as large as 50% for density ratios up to 7. These measurements were limited to mean density and streamwise velocity profiles and no details of the changes to turbulence and mixing properties are available. Our present investigation focuses on this flow but in a temporal configuration. The governing equations for this incompressible flow differ from those for a single fluid with thermal-induced density variations, as used by Pantano & Sarkar (2002) and, in a low speed limit, by Almagro et al. (2017). The relationship between the equations governing these flows has been reviewed in detail by Livescu (2020). Baltzer & Livescu (2020) focused this analysis on applications to mixing layer simulations and found that mean statistical profiles showed little difference when the density difference between free streams was compositionally-induced or thermallyinduced. However, these cases had significant differences in their mixing and density probability density function behaviors.
The present temporal simulations are relevant to understanding variable-density effects on growth in the spatially-developing configuration. 2D simulations of early-time spatially-developing mixing layers show strong differences in entrainment depending on whether the low or high speed stream has lower or higher density (Reinaud 2000;Joly 2002); we are unaware of any spatial simulations of fully turbulent growth. Based on experiments, Brown (1974) studied the thickness growth rate of variable-density spatiallydeveloping fully-turbulent mixing layers. He assumed that the temporal growth rate (i.e., from a frame of reference moving with the mixing layer convection velocity) is independent of the density difference between the streams, which is contrary to the reductions observed by Pantano & Sarkar (2002) and Almagro et al. (2017). As discussed in Pantano & Sarkar (2002), Brown (1974) combined this with the observation that the convection velocity is closer to the velocity of the high-density stream to propose a formula for growth rate reduction with Atwood number. Dimotakis (1984) refined the formula to account for asymmetric entrainment that is present only in spatiallydeveloping mixing layers. Ashurst & Kerstein (2005) studied variable density effects in temporal and spatial mixing layers using the one-dimensional turbulence stochastic simulation method; they captured many of the effects observed in Pantano & Sarkar (2002).
Other studies have addressed variable-density shear-driven mixing layers with buoyancy or other complicating physics playing a significant role. Olson et al. (2011) simulated mixing layers with mixed RT (buoyant) and KH (shear) instability and Atwood numbers ranging up to 0.71 using the same governing equations as for our present study. They focused on early times when complicated interactions between the instabilities produce complex effects on the growth rate. The linear stability study of Zhang et al. (2005) also considered a similar configuration. Barros & Choi (2011) performed linear stability analysis in a similar configuration representative of some environmental flows and highlighted the importance of the variable-density inertial terms beyond a Boussinesq approximation. Experimentally, Akula et al. (2013) studied mixed RT and KH instability with air and air/helium mixture streams shearing past each other, following a number of water-based experiments (also reviewed therein); buoyancy was the principal density effect and the Atwood numbers were low (< 0.04). Gat et al. (2017) simulated the mixing of vertical columns of fluid with different densities and perturbed interfaces. Gravity accelerates the perturbed heavy column downward within the triply-periodic domain to induce KH instability. Their configuration contains some of the same physics (shear aligned with buoyancy) as the more complex configuration of a buoyant jet, which was recently studied experimentally by Charonko & Prestridge (2017) and received more detailed analysis of the cascade of energy between scales by Lai et al. (2018). Additional multi-composition variable-density shear studies in the presence of other complicating physics include simulations of hydrogen and air streams to address supersonic turbulent combustion by O'Brien et al. (2014), reacting mixing layer simulations by Miller et al. (2001), and hybrid motor simulations with oxidizer and gasified fuel by Haapanen (2008).
Our present investigation seeks to elucidate the fundamental changes to the self-similar growth in a free shear flow produced by differences in the density of each stream with differing compositions. We perform direct numerical simulations in the simple incompressible temporally-developing configuration with two miscible fluids. In particular, we seek to quantify the asymmetries that appear in the flow statistics due to variable-density effects (whereas the analogous single-density incompressible temporal mixing layer configuration is statistically symmetrical) and explain their effect on growth characteristics. The paper is structured as follows: Section 2 describes the simulation approach and governing equations, followed by a description of the initial conditions in Section 3. Section 4 discusses flow properties that can be adduced from the governing equations and introduces definitions of flow measurements. Section 5 presents an overview of mean and fluctuation statistics from the simulations and relates growth rates to statistical profiles. Section 6 briefly addresses the local effects of density on velocity-related statistics, leading to the conclusions of Section 7. This is followed by appendices addressing (a) the relationship between density profiles and mean crossstream velocity and (b) contrasts between the present variable-composition flow and variable-thermodynamic-property flow.

Simulation Approach
The simulations are performed in the canonical temporal configuration, with two velocity streams of equal magnitudes flowing in opposite directions. The temporal configuration can be regarded as the limit of mean convection velocity of a spatial mixing layer approaching zero. In this case, the mixing layer develops with time instead of with spatial position as the flow convects downstream for the latter configuration. By using periodic boundary conditions in the streamwise (and spanwise) directions, the temporal configuration avoids the need for choosing inflow and outflow conditions and focuses on the variable-density effects on mixing in the simplest configuration possible. To our knowledge, this is the first study focusing on variable-density effects due to composition variation without additional effects such as compressibility, reactions, etc. Following the typical set-up (e.g., Rogers & Moser 1994), the coordinates are oriented such that 1 (x) denotes the streamwise direction aligned with the mean velocities, 2 (y) denotes the cross-stream (transverse) direction normal to the fluid interface, and 3 (z) denotes the spanwise direction (figure 1).

Governing Equations
To study incompressible mixing layers involving two fluid streams with strongly differing densities, the governing equations are formed by considering the full compressible flow equations for a miscible binary fluid mixture and then obtaining the infinite speed of sound incompressible limit (Livescu 2013). Gravity is not included here, but otherwise the governing equations are identical to those describing variable-density (non-Boussinesq) RT flow, as simulated by Cook & Dimotakis (2001), Livescu & Ristorcelli (2007) and Wei & Livescu (2012). To our knowledge, the present study is the first application of these equations to purely shear-driven variable-density fully-turbulent mixing layers.
The equations for the instantaneous variables (with partial derivatives denoted following the comma in the subscript, namely t representing the time variable t and an index i representing the relevant spatial direction x i ) are Figure 1. Variable-density mixing layer simulation set-up and coordinate system.
where the viscous stress, assumed to be Newtonian, is The governing equations are supplemented by slip boundary conditions in the y direction and periodic boundary conditions in x and z directions. Equation (2.3) represents the nonzero divergence of velocity that occurs due to the change in volume during mixing (while the flow is incompressible). The Fickian form with diffusion coefficient D represents the infinite sound speed limit of the full multicomponent diffusion operator (Livescu 2013). Equation (2.3) can be derived from the mixture rule ρ = 1/ (Y 1 /ρ 1 + Y 2 /ρ 2 ) (where Y 1 and Y 2 are species mass fractions of pure fluids with constant densities ρ 1 and ρ 2 , respectively) and species mass fraction transport equations for each species (ρY m ) ,t + (ρY m u j ) ,j = D (ρY m,j ) ,j (Sandoval 1995;Cook & Dimotakis 2001;Livescu & Ristorcelli 2007). The mixture rule can also be connected to the infinite speed of sound limit of the ideal gas mixture equation of state. Alternately, the same divergence relation can be derived as the infinite sound speed limit of the energy transport equation, which demonstrates the consistency of the VD governing equations (Livescu 2013). The dynamic viscosity of mixed fluid obeys a relation analogous to the density: µ = 1/ (Y 1 /µ 1 + Y 2 /µ 2 ), where µ 1 is the viscosity of the pure fluid with density ρ 1 and µ 2 is the viscosity of the pure fluid with density ρ 2 , which ensures a uniform Schmidt number, Sc = µ/(ρD), throughout the mixture.

Notations
Many of the statistics are based on averages, which are indicated by the symbol . Generically, the Reynolds mean of a quantity q is denoted by q and Reynolds fluctuation is q = q − q . For simple expressions, the Reynolds mean will also be indicated by an overbar, i.e.q, which is equal to q . As is typical for compressible flows, Favre averaging is employed for the mean governing equations to account for density variations. The Favre mean of a velocity component, u i , is denoted byŨ i = ρu i / ρ and the Favre fluctuation is u i = u i −Ũ i , in contrast to the Reynolds mean,Ū i = u i , and fluctuation, Numerical quantities presented in sections below are obtained from averages computed based on homogeneities present within the flows. Since the flow is periodic and homogeneous in the streamwise and spanwise coordinates x and z, area averages are computed across y-normal planes. Self-similar statistics will also be considered in which profiles should not change with time (except for noise due to lack of statistical convergence) when the y coordinate is scaled by an appropriate length scale. For these statistics, time averaging is also performed over the self-similar growth duration to improve statistical convergence ( §5.2). The averages computed to obtain Reynolds (Ū i ) and Favre (Ũ i ) averages are x − z area averages only when the statistic is a function of time or not in self-similar coordinates, but time averages are taken of the area averages when selfsimilar statistics are presented and the same set of notations is used for the averaged quantities.

Numerical Approach
The governing equations (2.1-2.4) are solved numerically using a pseudo-spectral scheme for spatial discretization in the periodic (streamwise and spanwise) directions and a compact difference scheme for the inhomogeneous (cross-stream) direction of the flow. The algorithm and code are slightly modified from those employed and described by Wei & Livescu (2012); Livescu et al. (2010Livescu et al. ( , 2011 for variable-density RT simulations; the equations solved are the same except non-zero mean streamwise velocity is present in the mixing layer. The cross-stream (normal) velocities at the lower and upper slip wall boundaries are maintained at zero, and this is consistent with the governing equations for this temporal mixing layer. Averaging the divergence equation (2.3) with diffusivity D assumed constant, and then omitting the terms of the summed indices that vanish due to the homogeneities present in the flow results in u 2 ,2 = −D ln ρ ,22 . (2.5) Integrating across the y domain, this expression becomes u 2 (y max ) − u 2 (y min ) = −D {[ln ρ] ,2 (y max ) − [ln ρ] ,2 (y min )}. Since density remains constant at the free streams existing at the upper and lower walls, it follows that u 2 (y max ) − u 2 (y min ) = 0. Thus, the variable-density equations are consistent with the boundary conditions u 2 (y min ) = u 2 (y max ) = 0. This argument also holds for thermally-induced single-fluid variable-density mixing layers (for which the governing equations are summarized and contrasted with the present equations in Livescu 2020; Baltzer & Livescu 2020) if the heat conduction coefficient is constant. More complicated cases such as heat release with chemical reaction necessitates nonzero normal velocity at the boundaries, e.g., Higuera & Moser (1994). Spatially-developing mixing layers also include streamwise gradients in the streamwise velocity, leading to another term remaining in the left-hand side of the divergence equation, which leads to cross-stream velocities at the upper and lower domain velocities associated with entrainment in even the single-density case. The third-order accurate variable time stepping Adams-Bashforth-Moulton method is used for time integration, coupled with the usual fractional step method. This is adapted for the pressure equation with variable coefficients due to non-zero velocity divergence associated with the variable-density equations. Fourier representations in the periodic coordinate directions allow the variable coefficient Poisson equation for pressure to reduce to an ordinary differential equation in the inhomogeneous direction. Taking advantage of the structure of the compact derivative, direct solvers can be employed for constant coefficient Poisson equations. The algorithm was initially devised for triplyperiodic buoyant turbulence simulations by Livescu & Ristorcelli (2007) to provide an exact divergence of momentum and thus avoid degrading the overall order of accuracy. This was an advancement from the algorithm used by Sandoval (1995) that required an extrapolation of velocity in time in order to determine the divergence of momentum but could degrade the overall temporal order of accuracy from second-order.
The variable coefficient Poisson equation for pressure is decomposed into the form ∇p/ρ (n+1) = ∇q + ∇ × A + L , which results in a constant coefficient equation corresponding to the dilatational (curl-free) component, ∇q, and implicit equations for the curl (divergence free), ∇× A, and mean components. The implicit equations are solved iteratively, using the direct Poisson solvers at each step. Due to the periodic boundary conditions, the mean term L is non-zero only in y direction. Differences compared to the RT algorithm appear in the mean term for the mixing layer because of the mean flow in the streamwise direction. For the RT case, the mean velocity is zero in both (periodic) horizontal directions, while for the mixing layer case, it is zero only in the (periodic) spanwise direction.
This algorithm avoids introducing additional errors that could affect mass conservation or degrade the accuracy from the time stepping method. The dilatational component of ∇p/ρ is related to mass conservation, which is enforced to machine precision due to the direct solvers involved. The curl component, ∇× A, is related to the baroclinic production of vorticity. The iterative procedure is performed until the maximum x-z planar average squared change in ∇ × A relative to the previous iteration value reduces to 0.01 times the squared value of ∇ × A averaged within the plane, for each component α: where n denotes the iteration number. This tolerance ensures small differences compared to convergence to machine precision. Note that each step of the iterative procedure is based on a direct Poisson solver.
No filtering was used in the simulations, so that the small scales are not affected by numerical artifacts. The spatial resolutions were determined by the requirement that the Kolmogorov scale is well resolved and a series of lower resolution, early time mesh convergence studies. The higher Atwood number cases have more stringent spatial resolution requirements, but for consistency, the same resolution was used for all simulations with Atwood number of 0.75 or below. Therefore, the lowest Atwood number simulations are over-resolved but should yield very high-quality vorticity and velocity gradient statistics. As described below in the discussion of self-similarity, at late times the peak local dissipation decays linearly with time, so the simulations require the finest resolution during the initial growth stage. Moin & Mahesh (1998) note that the Kolmogorov length scale is often cited as the smallest scale that needs to be resolved, but suggest that this requirement is more stringent than necessary for reliable first-and second-order statistics. For spectral methods, resolution is often expressed as k max η, where η is the average Kolmogorov length scale (ν 3 / ) 1/4 and k max = a(2π/L) for a spectral representation of N grid points in a domain of length L. The leading coefficient of the k max definition depends on the dealiasing employed, up to a maximum of N/2 if no truncation is used. The present simulations calculate the advective terms in skew-symmetric form to reduce the aliasing errors for cubic terms (Blaisdell et al. 1991). In DNS intended to maximize Reynolds number, typical values are 1 k max η 2 (Gotoh & Yeung 2013), with k max η ≈ 1.5 typical for adequately-resolved DNS of isotropic turbulence Pope 2000). Greater resolution may be required when special attention is focused on certain features, such as fine scale structure associated with stretched spiral vortices in isotropic turbulence that requires k max η 4 (Horiuti & Fujisawa 2008) or the alignment of strain rate and vorticity (Hamlington et al. 2008).
In the present mixing layer at negligible Atwood number, k max η for the Fourier spectral representation of each homogeneous direction reaches a minimum of ≈ 1.7 at early times at the centerline (where turbulence is most developed) and continuously increases thereafter. Pantano & Sarkar (2002) report final values of k max η ≈ 1.0, and they rely on spatial filtering that was shown to produce a relatively small amount of nonphysical dissipation to improve stability in their simulations. Resolution can also be be quantified in terms of grid spacing relative to the average Kolmogorov scale. Almagro et al. (2017) reported horizontal grid spacing finer than 1.8η during the self-similar growth, whereas the corresponding values in Pantano & Sarkar (2002) are 3-4η. In the present low A simulation, the horizontal grid spacing (∆x and ∆z) peaks at 1.8η during the earlytime transition and reduces to 1.0η during self-similar growth. Since the mixing layer is inhomogeneous and the Kolmogorov microscales shown above calculated from the dissipation at the peak y position does not account for inhomogeneities in the flow scales, these values merely represent a guideline.
For the present high Atwood number simulations, resolutions can be similarly estimated using the isotropic turbulence formula for η that does not address how scales may vary with local density variations. For the present A = 0.75 simulation, which has the same grid spacing as the A = 0.001 simulation, k max η attains a minimum value of 1.8 at early times and is 3.2 to 3.7 during the self-similar growth (which is similar to the values attained in the A = 0.001 case). For A = 0.75, the horizontal grid spacing corresponds to a maximum of 1.8η at early time and decreases to 1.0η by the end of self-similar growth. For A = 0.87, the simulation requires a greater number of grid points for the same physical domain size to maintain numerical stability. The calculated k max η reaches a minimum value of 2.7 at early times but remains between 4.4 and 5.3 during the identified self-similar growth interval. The horizontal grid spacing corresponds to a maximum of 1.2η at early time and decreases to 0.6η by the end of self-similar growth for A = 0.87. Nonetheless, these values based on isotropic turbulence η are not sensitive to localized steep velocity and density gradients at increased Atwood number that are hypothesized to necessitate greater resolution for numerical stability.
The compact finite difference scheme used for the cross-stream (y) direction is 6th order accurate for both the momentum and pressure equations. The uniform grid spacing is finer (reduced to a factor of 0.8: ∆ y = 0.8∆ x = 0.8∆ z ) in the inhomogeneous direction, in order to compensate for the lower accuracy relative to the Fourier directions. Modified wavenumber analysis for 6th order compact difference equations indicates errors in differentiating modes become larger at higher wavenumbers . Since differentiation with the Fourier method is exact up to its highest resolved wavenumber, the Fourier method has no error until the Nyquist frequency. This corresponds to a grid spacing of 2η if k max η=1.5. Requiring the compact difference method to produce less than 25% error in differentiating a mode with this same wavelength dictates that the grid spacing must be refined relative to that of the spectral method by a factor of 0.8. Note that the vast majority of the energy in the flow is at longer wavelengths that have negligible error, according to the modified wavenumber analysis: the lowest 3/4 of the wavenumbers have errors of less than 3.5%.
The pressure determined by the fractional step method restores the velocity field divergence to be consistent with (2.3); however, it represents the average pressure over the time step. To recover the instantaneous pressure for calculating budgets and other statistics, the Poisson equation resulting from obtaining the divergence of (2.2) is computed as a post-processing step after the flow has been advanced in time by the fractional step method. The numerical algorithm has been verified to accurately satisfy the governing equations by comparing the time derivatives calculated for various quantity budgets that appear throughout this paper with the appropriate budget right-hand sides.

Domain Size
The domain lengths in the homogeneous streamwise and spanwise directions L x and L z are directly related to the convergence of statistical quantities obtained by planar averaging. In addition, these dimensions potentially affect the sizes of structures that grow within the domain. Convergence can be improved either by enlarging the domain size or by using an ensemble of smaller domain simulations. However, a sufficiently large domain is necessary to achieve correct structure growth and interactions.
Several domain sizes were tested and the final dimensions used were found to have minimal evidence of structure growth restriction compared to smaller sizes. From the perspective of initial KH rollup structures with an assumed streamwise wavelength of the most unstable linear instability mode λ ls , the present mixing layer domain accommodates 64λ ls in the streamwise direction. This corresponds to 6 successive mergers; Vreman et al. (1997) found that lengths of 8λ ls (i.e., three successive mergers) were required to reach reasonable self-similarity. In shear flows, the longest scales are oriented along the streamwise direction. The domain therefore has a L x /L z ratio of 4, which was adopted by a number of previous temporal mixing layer simulations (e.g., Rogers & Moser 1994;O'Brien et al. 2014).
The cross-stream domain size, L y , must also be sufficiently large that the mixing layer evolves freely without the slip walls at the y domain boundaries influencing the growth. A series of simulations with different thicknesses has been performed to ensure the statistics are not influenced by the walls for the self-similar time of interest. The initial interface is positioned so that it is nearer the heavy-fluid wall than the light fluid wall in proportion to the Atwood number, since the mean velocity neutral point (interface center) and the most intense turbulence drift to the light fluid side as the flow develops ( §5.3). The interface is centered within the domain for the A = 0.001 case (as this effect is negligible at low density ratios). The domain sizes are summarized in Table 1. Although initial momentum thickness δ m,0 (defined below) is somewhat ill-defined for making comparisons, comparing L x /δ m,0 suggests that the domain lengths are approximately 10 times those of Pantano & Sarkar (2002)

Initial Conditions
Mixing layer simulations are typically designed either to approximate a physical mixing layer experiment or to be in a generic configuration commencing from a simple disturbance. The latter approach is here adopted for generality and to promote quickly reaching self similarity without artifacts from the initial condition. Nonetheless, parameters are broadly within the range of those found in experiments.

Mean Velocity and Density Profiles
The initial mean velocity profile that approaches the free-stream velocities of ±∆U/2 at the y boundaries is specified as where the momentum thickness δ m,0 specifies the initial thickness of the interface. The hyperbolic tangent profile is commonly used in a wide range of mixing layer simulations, such as Riley et al. An initial density profile is prescribed to specify the differing compositions (and thus densities) of the fluid streams. The simulations focus on the simplest case of two separate streams of different velocities and densities meeting at a thin interface, so the initial density profiles are aligned with and of the same thickness as the velocity profiles. Thus, the initial density profile is with density profile thickness δ ρ,0 chosen to equal δ m,0 . This specification of aligned tanh profiles of density and velocity is similar to the approach of Pantano & Sarkar (2002) and Almagro et al. (2017), though their density variations were attained by varying the thermodynamic properties for a single fluid. In either approach, the mean density of the lower and upper streams of fluid ρ 0 = (ρ 1 + ρ 2 ) /2 is matched between all of the simulations within the set. The desired Atwood numbers A are then attained by specifying free-stream densities ρ 1 = ρ 0 − ∆ρ/2 and ρ 2 = ρ 0 + ∆ρ/2, where ∆ρ = ρ 2 − ρ 1 = 2Aρ 0 . Symmetries present in the temporal mixing layer (but not the spatially-developing case) result in the flow behaviors being equivalent whether the negative mean streamwise velocity is associated with the light fluid and the positive velocity is associated with the heavy fluid or vice versa, as also noted by Pantano & Sarkar (2002). Thus, results from a different profile convention can be compared by selecting coordinates to match density profiles and then changing the sign of the mean streamwise velocity to also match.

Initial Disturbance
Only the velocity field is perturbed relative to the mean profile given above to induce the transition to turbulence. This is appropriate because the velocity field drives the instability and turbulence, as observed in the single-density case; this approach also allows the disturbance to be consistent between Atwood numbers. Different velocity disturbances can produce significantly different growth rates at early times in mixing layers (Fathali et al. 2008), but the present goal is to quickly establish self-similar growth and minimize long-lived large-scale structures that are uniquely associated with initial disturbances. To roughly resemble physical experiments, the velocity perturbation is confined to a thin (in y) region centered at the mean velocity profile interface.
In the present simulations, this is accomplished by generating a random field (filling the full domain) that is divergence-free and has a 3D energy spectrum obeying a Gaussian behavior at high wavenumbers with k 4 behavior at low wavenumbers as E(k) = (k/k 0 ) 4 e −2(k/k0) 2 . Here, k = k 2 1 + k 2 2 + k 2 3 is wavenumber and k 0 is the prescribed peak wavenumber. k 0 is selected to be λ ls /4, where λ ls is the streamwise wavelength of the least stable mode calculated from temporal linear stability analysis for the base velocity profile (λ ls = 28δ m,0 for the present set-up). This places much of the energy at small scales to quickly establish turbulent motions. The disturbance spectrum is that used by Pantano & Sarkar (2002) and the positions of the peak wavelength (relative to the least stable wavelength) are similar. The field is then tapered to a thin interface region by multiplying by the Gaussian profile in y to obey u i u i (y) = Ae − 1 2 (y/σ) 2 , where σ is the intensity profile thickness chosen to be 2δ m . This is nearly equivalent to the thickness used in Riley et al. (1986) simulations based on measurements of the intensity profile in a mixing layer experiment and to the thickness used by Pantano & Sarkar (2002). The peak amplitude A is specified for peak intensity u i u i of 0.03∆U 2 by prescribing a 0.1∆U RMS fluctuation for each velocity component. This relatively strong disturbance reduces the time to reach self-similar growth. The self-similar value of the streamwise turbulent velocity fluctuation intensity reaches approximately 2.5 times this initial value.
This initial velocity disturbance is similar to those used by Riley et al. (1986) (further described in Riley & Metcalfe 1979) and Pantano & Sarkar (2002) (further described in Pantano-Rubino 2000), but details of the implementations differ. The present approach of multiplying the field by the y-intensity profile produces divergence, which is corrected by applying the pressure step of the projection method to the velocity field. This step slightly weakens the intensity of the u 2 velocity component. Alternatives exist (e.g., applying the profile to a vorticity field, thereby producing a divergence-free velocity field as in Pantano-Rubino 2000), but the present method produces an initial velocity field divergence fully consistent with the variable-density incompressible divergence condition (2.3). A small mean u 2 velocity is also produced by this step, which is consistent with the divergence condition (as further explained in Appendix A). This mean velocity is concentrated at the interface and decays toward the y boundaries; the magnitude is also very small (< 1% of ∆U in all simulations shown).

Viscosity and Diffusivity
Momentum thickness Reynolds number, Re m = ∆U δ m /ν, can be maximized during the self-similar stage by either growing to a large final thickness δ m or having a small viscosity ν. The initial configuration is chosen to maximize the thickness growth so that the fully-turbulent state is less affected by the initial disturbance. This is achieved by selecting a relatively small initial momentum thickness and appropriate viscosity such that all scales are well resolved and the initial growth is not overly damped. The fundamental velocity scale ∆U to initialize the simulation is arbitrary and can be scaled out. In consistent units, ∆U = 1 is prescribed with initial momentum thickness of 0.5 and viscosity of 0.00625. This initialization results in a Reynolds number Re m of 80; however, this value has limited meaning before mixing layer evolution sustains the scales of motion.
The Schmidt number Sc = ν/D is chosen to maintain a constant value of 1 everywhere as the fluids mix. This is imposed by selecting the same values of kinematic viscosity ν = µ/ρ for each of the participating fluids (i.e., ν 1 = ν 2 ) with constant diffusivity D. The choice of constant kinematic viscosity to maintain constant Schmidt number of 1 is frequently used in other multi-fluid mixing studies (e.g., Sandoval 1995;Cook & Dimotakis 2001;Livescu & Ristorcelli 2007), though maintaining Sc = 0.7 (which is typical for gases) is also common (e.g., Olson et al. 2011). Note that the choice of constant ν implies that µ ∼ ρ, whereas with real fluids there is typically a weaker dependence on density such as µ ∼ √ ρ ).

Basic Definitions and Theoretical Flow Properties
While detailed simulations are necessary to obtain many quantities describing the flow, several characteristics of the flow can be deduced from the governing equations and flow configuration. The Favre mean equations obtained from (2.1-2.2) arē where the Favre Reynolds stresses arẽ These equations apply to incompressible variable-density flows as well as fully compressible flows. When the equations are applied to the geometry and flow conditions of the temporallydeveloping mixing layer, many of the terms vanish due to homogeneity and symmetries of the flow. The expanded equations after these simplifications arē ,2 + (ρR 12 ) ,2 =τ 12,2 (4.5) The slip wall boundary condition in the y direction requires thatŪ 2 =Ũ 2 = 0,R 12 = 0, andτ 12,2 = 0 at the boundary. These conditions are consistent with the variations outside the mixing layer, whereρ andŨ 1 are constant. As shown in Appendix A, for the incompressible flow considered here, the mean cross velocity can be expressed solely in terms of density moment statistics and their derivatives; the cross-stream velocity is necessarily zero if the flow contains no density variations.

Conservation Properties
Integrating the mean density conservation equation (4.4) over the y domain indicates that y2 y1ρ dy is constant with respect to time (total mass within the domain is conserved). The mean momentum equations (4.5)-(4.6), when similarly integrated over the y domain, show that y2 y1ρŨ i dy are also constant with respect to time (total momentum within the domain is conserved), when the remaining terms vanish at the boundaries. This is approximately satisfied for (4.5) and (4.6) throughout the duration of the simulation, since the velocity fluctuations remain at low values near the slip walls, and therefore the advective term and Reynolds stress are negligible at the y domain boundaries, while the mean pressure gradients and viscous stresses have relatively little effect.

Self-Similarity
Another property expected of mixing layers is attaining states of self-similar growth. For the temporal configuration, the statistics are functions only of time and the inhomogeneous y position. Assuming self-similarity and that both mean density and velocity profiles are initially centered at y = 0 (so that no other length scale is introduced in the problem), the time-and y-dependencies are eliminated by introducing a new variable η = y/h, where for the present purpose, h generically represents a length scale that characterizes the y-thickness of the mixing layer and grows with time. Specific choices for defining this thickness are discussed below. The scaled coordinate η defined here is separate from the Kolmogorov length scale η of §2.3.
As described in Appendix B, the mean mass conservation equation (4.1) and Favre mean streamwise momentum equation (4.2) are satisfied for self-similar growth when the growth rate dh/dt is constant and the mean variables are non-dimensionalized as Analyzing the resulting self-similar mass conservation and streamwise momentum equations (Appendix B) reveals relations between the scaled y positions at which features in the statistical profiles occur. Let η 2 be defined as the η point where the Favre crossstream velocity inflection point occurs [dÛ 2 /dη(η 2 ) = 0] and η 12 as the point where Favre shear stress has its inflection [dR 12 /dη(η 12 ) = 0]. Then the self-similar analysis proves that η 12 < η 2 < 0. That is, the Reynolds stress peak is located further in the light fluid than the peak of mean cross-stream velocity. This analysis does not determine the position η 1 of the zero-crossing of Favre streamwise velocity [Û 1 (η 1 ) = 0], but this can be empirically investigated in the simulations.
The above analysis and arguments reach similar conclusions to those presented by Pantano & Sarkar (2002) after developing the self-similar analysis framework while analyzing their variable-density flow. It should be noted that these self-similar equations and results pertain to any variable-density mixing layer that obeys the compressible mass conservation and streamwise momentum equations. Specifying particular cases of the flow (in this case, incompressible binary mixing of species, as opposed to thermodynamic variations or high-speed flow) influences the specific forms of the self-similar quantities (assuming states of self-similar growth are reached).

Thickness Definitions
The thicknesses of the density and streamwise velocity profiles are among the most basic global quantities characterizing mixing layers growth. Though the density and velocity mean profiles initially coincide, they need not grow identically as the flow evolves, so various thickness measurements are defined based on both profiles.
Thickness of a mixing layer is traditionally quantified based on the mean streamwise momentum profile, which has a clear connection to the momentum equation (4.2). Momentum thickness is defined as: As the first form emphasizes, this corresponds to the integral of the product representing deficits relative to free streams, which have streamwise velocities of U − = −∆U/2 and U + = ∆U/2. An analogous thickness could also be defined on a per-mass basis to depend only upon the mean velocity profile: This definition uses Reynolds-averaged streamwise velocity rather than Favre-averaged to avoid any explicit dependence on the density field. For single-density mixing layers, (4.12) is commonly given as the definition of the momentum thickness because (4.11) reduces to this when density is constant, though (4.11) is the most formal definition. Several other quantities also are commonly used to characterize mixing layer thickness based on the mean velocity profile, but these are generally less smooth (i.e., more sensitive to lack of statistical convergence) than the integral thicknesses defined above. These other measurements include lengths based on gradients of profiles. Vorticity thickness is obtained from gradients of the Reynolds mean velocity profiles as as the vorticity magnitude reduces to |dŪ /dy| in the absence of a mean streamwise gradient in cross-stream velocity. This measure based on only a small portion of the mixing layer (where the mean gradient is steepest) has the potential to produce a misleading representation of the thickness of the layer when significant asymmetries are present. The distance between positions at which the mean velocity reaches specific percents (e.g., 10%) of the difference ∆U between its free-steam values U − and U + (which are associated with fluids having densities ρ 1 and ρ 2 , respectively): While momentum thickness and vorticity thickness have been the most commonly used thickness measurements in the historic mixing layer literature, Pope (2000) adopts h 0.1 in treating planar mixing layers, and it has also been recently used by Schwarzkopf et al. (2016), for example. For brevity, h will be used herein to indicate h 0.1 . This choice of velocity percent produces measurements that are smoother and less sensitive to statistical fluctuations than selecting a smaller fraction (e.g., h 0.01 ) that would yield thicknesses based on the flow far out in the intermittently turbulent / non-turbulent interface. Favreaveraged velocity is used for h, though it could alternatively be based on Reynoldsaveraged velocity, as could any of the other thickness quantities. For even the highest Atwood numbers, the effect of averaging type on the calculated thickness is negligible: using Reynolds averages instead of Favre averages for A = 0.87 produces about 1% larger values for h and 5% larger values for δ m . Favre averaging is used for all of the velocity-based thicknesses shown except for δ m,pm and δ ω .
For variable-density mixing layers, similar thicknesses may be defined based instead on the density profiles as Note that δ ρ is equivalent to the width measurement introduced by Youngs (1991Youngs ( , 2009 and also used by Livescu et al. (2010); it is typically written as W = ∞ −∞ F 1 F 2 dy, defined based on the mean volume fractions of each species F 1 = (ρ − ρ 1 ) / (ρ 2 − ρ 1 ) and F 2 = (ρ 2 −ρ) / (ρ 2 − ρ 1 ). Typically, a scaling constant β is used with W to approximate bubble height in RT flows as h * b = βW ; β depends on Atwood number in order to represent asymmetries that develop in RT flow structure as Atwood number increases (Youngs 2013;Livescu et al. 2010). For brevity, h ρ will be used herein to indicate h ρ,0.1 . As shown below, the mean density profiles develop significant asymmetries at high Atwood numbers, which implies that (4.17) can only accurately represent the layer thickness at very low density ratios.
Additional width quantities commonly used for variable-density flows (particularly RT instabilities, e.g., Zhou & Cabot 2019) are also relevant. One such quantity, used by Cook & Dimotakis (2001) where X P represents the amount of product in a hypothetical fast reaction between the two species (4.18) X P (ρ) corresponds to the mole fraction of fluid fully mixed to the mean density. Physically, h Xρ is the thickness of mixed fluid that would result if the two fluids were perfectly homogenized within the mixing layer.

Time Evolution of Mean Profiles and Thickness Growth
Area-averaged mean profiles of streamwise velocity and density illustrate the basic properties of the mixing layers' evolution with respect to time. These profiles are shown for two representative Atwood numbers (almost single-density and strongly variable density) in Figure 2. These profiles form the basis for the thickness scales defined in §4.3. Figure 3 displays the time evolution of thickness by several definitions involving the above profiles. All measurements indicate that simulations for each Atwood number approach linear thickness growth with respect to time at late times. Regardless of the specific thickness definition, thickness growth is retarded with increasing Atwood number. The momentum thickness (a) indicates a strong reduction in growth with Atwood number, whereas the momentum thickness per mass (b) and h (c) quantities both indicate weaker growth reduction, as does vorticity thickness (not shown). The thickness evolutions also highlight that the mixing layers grow to many times their initial thicknesses, as desired to reach self-similar growth.
It should be noted that δ m,0 , used for nondimensionalization, is based on initially aligned profiles at t = 0, before the shifts of mean streamwise velocity relative to mean  density have developed. Thus, δ m and δ m,pm are initially essentially equal but evolve differently as the profile shifts develop. The correspondence between initial δ m and other initial length scales is h 0 = 4.39δ m,0 and δ ω,0 = 4δ m,0 ; similar relations apply to the analogous initial density thicknesses δ ρ , h ρ , and δ dρ/dy as well. However, as the profile shapes evolve in transition and turbulent flow, these relations no longer apply.
To evaluate whether constant values for self-similar temporal growth rates are reached, the time derivatives of thicknesses are shown as functions of time for each Atwood number in Figure 4. Thicknesses based on integral measures produce relatively smooth growth rates that in each case asymptote to constant values at late time. Growth rates based on h contain more noise than the rates based on the integral quantities, but applying a Hann filter to smooth the thickness vs. time functions produces the result shown in Figure 4c. These results are also consistent with asymptoting growth rate (though statistical fluctuations are present). For vorticity and density gradient thicknesses, calculating the gradient of a mean profile and then extracting its y-maximum makes these measurements more sensitive to noise associated with lack of statistical convergence. The sensitivity of the gradients to small-scale noise dictates that a small amount of spatial smoothing (via a Hann filter) first be applied to the instantaneous mean profiles to remove the finest scales of noise before calculating peak gradients.

Determining the Time Interval of Self-Similar Growth
In addition to constant growth rate, another consequence of self-similar growth is the statistical profiles collapsing when appropriately scaled. For example, the mean streamwise velocity and density profiles would collapse to single curves for all times during self-similar growth when y is scaled by thickness (e.g., δ m or h). As observed by Rogers & Moser (1994), mean velocity profiles are relatively insensitive to deviations from self-similar growth. However, fluctuation intensity profiles generally continue to converge after the mean velocities reach their self-similar profiles. Statistical profiles for many quantities are expected to have constant peak values and thus linearly increasing integral values as thickness grows linearly with time. Directly evaluating the time histories of statistics' peak values comprises a more stringent test of self-similarity, but evaluating their corresponding integral quantities instead is less sensitive to noise.
One statistic that is meaningful for evaluating self-similar growth is integral of crossstream velocity fluctuation intensity: (5.1) In earlier simulations emphasizing roll-ups of KH vortex structures and their subsequent mergers, Moser & Rogers (1993) showed that large values of V are associated with these features. Conversely, when Rogers & Moser (1994) began a mixing layer simulation from a fully turbulent field, no large values were attained but instead V slowly increased and then asymptoted to the self-similar value. Attili & Bisetti (2013) examined V for their spatially-developing mixing layer beginning from a thin disturbance (similar to that for the present simulations). It overshot the self-similar growth value when the vortices played an important role at early time, but decreased and asymptoted thereafter as the mixing layer reached a self-similar growth regime. This behavior is compared to that of the present simulation with negligible Atwood number in Figure 5a. The present simulation produces a much weaker peak in V than the Attili & Bisetti (2013) simulation. Despite the weaker peak, the present simulation follows similar behavior of approaching selfsimilarity after the peak. This behavior contrasts with the asymptoting from below that appears to occur for the fully-turbulent initial condition of Rogers & Moser (1994). All of the simulations shown in Figure 5 display V values remaining approximately constant throughout their respective self-similar growth periods, and these values are in good agreement between the simulations. In the present simulations, similar behavior also occurs at increased Atwood numbers. An important indication of self-similarity employed by Rogers & Moser (1994) is total dissipation of turbulent kinetic energy (TKE), which is planar-averaged dissipation ε = − τ ij u i,j (from the TKE budget equation) integrated across the entire mixing layer: The rate at which TKE ultimately is dissipated is set by the large-scale motions that scale (in magnitude) with the velocity difference between streams ∆U . Since E has units of velocity cubed, it can be argued on dimensional grounds that E scales with ∆U only and therefore is constant with respect to time during self-similar growth (Rogers & Moser 1994). Unlike the velocity fluctuation intensities, the dissipation peak value does not remain constant with respect to time but instead decays in magnitude proportionally with the mixing layer thickness. Thus, its integral over the increasing width as the mixing layer thickens remains constant.
For the essentially single-density case, the dissipation evolution is compared with those of other mixing layer simulations in Figure 5b. The self-similar growth durations are marked as identified in each corresponding reference. Depending on the route of transition, the peak dissipation may also correspond to an overshoot in dissipation prior to self-similar growth or to part of the self-similar growth regime. The former scenario applies to the simulation of Attili & Bisetti (2013) that begin from a thin disturbance. The latter applies to the simulation of Rogers & Moser (1994) that begins from a field containing fully-turbulent fluid and slowly approaches the self-similar state from below (in terms of dissipation). Attili & Bisetti (2013) discuss these differences and the role of KH structures in the transition in further detail. The present flow corresponds to the former scenario, beginning from a thin disturbance leading to structures that cause dissipation to overshoot, though this is weaker than in Attili & Bisetti (2013) likely due to the form of the disturbance and the temporally-developing nature of the flow.
Compared to the close agreement of self-similar V value with the other simulations in the literature, there is significantly more variation among the self-similar integrated dissipation values. However, the Attili & Bisetti (2013) mixing layer appears to be asymptoting to a value near that observed in the present A = 0.001 simulation. The self-similar time interval shown for this present simulation (for which E is one of the determining considerations) maintains E to a nearly constant value.
The dimensional argument described above for constant E in self-similar growth holds for the variable-density mixing layers as well. For variable-density mixing layers, the TKE budget equation terms are often defined to include density (e.g., ), unlike the typical budgets written for single-fluid incompressible mixing layers (e.g., Rogers & Moser 1994). Therefore, the integrated dissipation must be divided by density to have the units of (∆U ) 3 . One option is to nondimensionalize by ρ 0 , the average of the two streams. However, the most typical treatment is to divide ε by the mean densitȳ ρ, in analogy to Favre averaging other quantities: Figure 6 demonstrates that E andẼ become constant in self-similar growth for each Atwood number. The values for E scaled by ρ 0 and ∆U 3 decrease strongly with increasing Atwood number, whileẼ scaled by ∆U 3 displays a much weaker dependence. While linear growth of thickness and constant integrated dissipation are key indicators of self-similar growth (which have been long been employed, e.g., Rogers & Moser 1994), comparing additional flow statistics profiles produces further useful indications. This was recognized by Vreman et al. (1997), who determined mixing layer growth to be selfsimilar when "the development of the shear layer thickness is linear in time and profiles of normalized statistical quantities at different times coincide." The time evolutions of profiles can be evaluated by monitoring the peak values of these statistics or examining their integrals in y divided by the thickness (as with V). This latter approach is less sensitive to statistical variability than the peaks. A number of profile quantities are considered in determining the self-similar growth time interval; integral velocity variances and Reynolds stresses are shown in Figure 7(a-c), while additional profiles (e.g., crosscorrelations between velocity and density) are considered but not shown for brevity. For each Atwood number, the integral turbulence intensities match very closely with the corresponding integral Favre-averaged Reynolds stresses and are nearly identical for A = 0.5 and below. Comparing between Atwood numbers, there is a consistent trend to lower intensities with increasing A during transition (when the values peak); during self-similar growth, the trend is weak and easily obscured by statistical variability. The y-integrated values shown may conceal some of the complexity in weakly changing profile shapes. For the cross-stream component (Figure 7b), the intensity increasing at late time is hypothesized to be associated the turbulent fluctuations reaching and accumulating near the slip walls to affect the interior of the mixing layer. This is expected to occur soonest for the lowest Atwood numbers because they experience the fastest growth. The self-similar time interval is determined to end before this phenomenon affects the flow.
Variable-density mixing layers introduce additional quantities to be considered for selfsimilarity, most importantly the density fluctuation intensity ρ ρ . The integral values of this planar-mean quantity are shown for each Atwood number in Figure 7d. ρ ρ can remain within a tolerance of a constant value later than other statistics and thus determine when the self-similar interval begins. These profiles are related to the mixing of the two streams, which is dependent on how fluid is transported into the cores of the mixing layers. Despite the complex mixing behavior, the simulations indicate that the density fluctuation intensity profiles for each Atwood number approach a unique self-similar scaled profile that remains approximately constant with respect to time.
The integral ρ ρ for A = 0.75 (blue curve) is suggestive of reaching self-similar growth at particularly late time, with a leveling occurring at earlier time before it again increases and levels off. It appears that the flow configuration changes during the second period of rapid increase. This behavior is responsible for the late starting time of the selfsimilar period. This increase in maximum density fluctuation intensity also appears to be associated with a smaller increase in integral cross-stream component velocity fluctuation intensity, as shown in Figure 7b. In summary, the self-similar periods are determined by seeking constant thickness growth rates, constant values of integrated dissipation, and statistical profiles that remain constant when the cross-stream coordinate is self-similarly scaled. In addition to the velocity intensity profiles, density fluctuation intensity profiles must also be considered for variable-density mixing layers. To identify self-similar growth periods in a consistent manner for all Atwood numbers, these conditions are approximated by requiring that thickness growth rates as well as integrals across the cross-stream domain of dissipation, velocity fluctuation intensity u i u j , and density fluctuation intensity ρ 2 be constant to within a specified threshold. The integrals of fluctuation intensity profiles are scaled by thickness (h) to attain constant values (or equivalently are integrated with respect to y/h), since the integrals would grow proportional to thickness if the self-similar scaled profiles remain constant. Mean profile convergence is accomplished by ensuring the more sensitive fluctuation intensity profiles are converged. This algorithm is consistently applied by determining the longest time interval that each of the quantities specified above remains within 10% of any value and then retaining the intersection of these time intervals as the self-similar time interval. The very large simulations produce satisfactory adherence to a relatively stringent set of criteria that must be simultaneously satisfied, as indicated by the self-similar periods marked in Figure 6. The self-similar periods for other simulations compared in Figure 5 are taken from their respective publications. Due to the effects of differing initial momentum thicknesses (and how they relate to the disturbances), the scaled times t∆U/δ m,0 (or scaled downstream position x/δ m,0 for the spatial-developing case) in this comparison cannot be meaningfully related between simulations. The significantly smaller domains that were feasible for many previous studies could contribute to the difficulties reported in reaching self-similarity (e.g., Vreman et al. 1996Vreman et al. , 1997Pantano & Sarkar 2002). In general, questions remain about the universality of the self-similar state (e.g., Dimotakis & Brown 1976;Rogers & Moser 1994;Vreman et al. 1997). However, the thin and broadband disturbance is intended to reduce idiosyncratic large-scale vortices that persist after transition as a result of the initial condition so the present simulations reach generic self-similar states.
Another consideration relevant to the self-similar growth regime is flow Reynolds number. For the flow statistics to be representative of the fully turbulent mixing in practical applications, the Reynolds numbers must be sufficiently large throughout the averaging time duration. In general, significant changes in mixing behavior have been observed to occur at a Reynolds number threshold (i.e., the mixing transition, Dimotakis 2000). Relevant Reynolds numbers are typically defined using the mixing layer thickness or the Taylor microscale. Both scales continuously grow as the mixing layers thicken with time. According to Dimotakis (2000), general necessary conditions for passing the mixing transition for turbulent flows are that the outer-flow Reynolds number exceeds Re ≈ 1-2 × 10 4 and that Taylor Reynolds number exceeds Re λ ≈ 100-140. Dimotakis defines the former Reynolds number using a visual thickness scale δ sh that is used in experiments; it has been estimated as δ sh ≈ 2δ ω for numerical simulations (e.g., Rogers  1994). This criterion corresponds to attaining Re ω ≈ 0.5-1 × 10 4 . Table 2 confirms that this condition is satisfied for the self-similar growth statistical averaging periods. The decrease of Re m values with Atwood number is a consequence of δ m decreasing as the velocity profiles shift into lighter density fluid. This complicates interpreting Re m in variable-density mixing layers. Though Taylor microscale is anisotropic in its most fundamental definition, it is estimated using a relation that strictly only applies to homogeneous isotropic turbulence, λ g = 10k ν ε . (Averaging the homogeneous-coordinate components of the fundamental Taylor microscale shows good agreement with this estimate for the present mixing layers.) The velocity scale is also taken as u rms ≈ ( u 1 u 1 + u 2 u 2 + u 3 u 3 ) /3 = 2k/3. Using the turbulent kinetic energy and dissipation at the y position of most intense turbulence, the estimate of Taylor microscale Reynolds number is Re λ =k 20 3 1 εν . Using u rms produces consistency with the velocity scale used in the λ g definition as well as consistency between the turbulent kinetic energy and dissipation included in turbulent kinetic energy budget (in analogy to isotropic turbulence). Though similar definitions are also used for other relevant flows (e.g., Sekimoto et al. 2016), mixing layer literature often uses √ 2k as the velocity scale (rather than u rms = 2k/3) to form Re λ = (2k)λ g /ν = k 20/(εν) (e.g., Pantano & Sarkar 2002;O'Brien et al. 2014;Almagro et al. 2017). Renormalized to the present convention, the Re λ range during self-similar growth for the single-density mixing layer of Rogers & Moser (1994) is 84-99 and for Almagro et al. (2017) is 81-87, for example. The present simulations generally satisfy the Re λ ≈ 100 (with Re λ is defined in this way) mixing transition guideline given by Dimotakis (2000) before their self-similar growth periods end. The consistency of the statistics within the self-similar growth periods suggests the turbulence is well-developed throughout. The initial condition that produces rapid transition is expected to lead to this state more quickly than the large-scale features that persist through other mixing layers' transitions. Figure 8. The times included in the plots correspond to the self-similar growth regimes, for which the determination is explained below ( §5.2). Figure 8 demonstrates that the time series of mean streamwise velocity and density profiles collapse to single curves when the cross-stream coordinate is scaled by the thickness measurement h. Similar collapse is also observed when the cross-stream coordinate is instead scaled by δ m , δ m,pm , or δ ω . While δ m was used as the thickness length scale in the discussions above to allow comparison with other studies, scaling statistics in terms of the h scale offers interpretive advantages in variable-density flow. For consistency, h will be used as the thickness scale henceforth, except for when making certain comparisons with other studies. The collapse of mean profiles is one indication that self-similar growth is achieved. During self-similar  For Favre mean streamwise velocityŨ1 (a-b, e-f) and scaled density (c-d, g-h), scaling the cross-stream coordinate by the initial thickness h0 shows only the self-similar time interval curves from Figure 2 and demonstrates the growth of the mixing layer (a,c,e,g), whereas scaling instead by each curve's thickness h causes this same time series to collapse to a single curve for each simulation (b,d,f,h). growth, it is thus appropriate to time-average the scaled profiles to improve statistical convergence. This averaging is also applied to all of the other scaled statistics presented below.

Time-Averaged Self-Similar Statistical Profiles
Comparing the self-similar scaled profiles among Atwood numbers (Figure 9) illustrates several basic changes that occur as the density difference between streams increases. For A = 0.001, the mean streamwise velocity and mean density profiles are essentially centered at y = 0 and symmetric about that point. A shift in the mean streamwise velocity profiles to the light fluid side (i.e., η 1 < 0) that increases in magnitude with increasing Atwood number is apparent. With increasing Atwood number, the shapes of these velocity profiles remain generally similar as they shift to the light fluid side, but the asymmetry in their gradients (Figure 9b) reveals an additional steepening on the light fluid side and shallowing on the heavy fluid side. Conversely, the neutral points of the density profiles (whereρ(y) = ρ 0 ) remain relatively stationary while the density profiles steepen on the heavy fluid side but become shallower on the light fluid side with increasing Atwood number. Figure 10 displays the corresponding profiles for the cross-stream mean velocity component. The magnitudes are much smaller than those of the streamwise velocity. However, as the self-similar analysis indicates, the cross-stream velocity has an important relationship with mass conservation and mixing layer growth in variable-density mixing layers. In Figure 10(b-c), these velocity profiles are shown with the scaling suggested by the selfsimilar analysis, using h based on the Favre mean streamwise velocity for the thickness scaling. The Reynolds averaged cross stream velocity is much smaller in magnitude than the corresponding Favre average quantity. In addition, V can be shown to strongly depend on the mean density gradient (Appendix A) and therefore not reach a time-constant magnitude during self-similar growth; the averages in Figure 10(c) should be understood to pertain only to their particular averaging time periods. It is shown below ( §5.6) that V is dominated by the turbulent mass flux, which does approach a constant value during self-similar growth.
The positions of the neutral points (i.e., ρ = ρ 0 andŨ 1 = 0) and positions of extrema for various statistical quantities (e.g., min(Ũ 2 )) are important in characterizing the shape of the mixing layer during the self-similar regime. The mixing layer growth and its asymmetry can be summarized by tracking the points at which the mean streamwise velocity is equal to 10 and 90 percent of the free-stream difference ∆U : y [Ũ1=U−+0.1 * ∆U] and y [Ũ1=U+−0.1 * ∆U] . These are the points whose separation define h in (4.14). In Figure 11a, the linear growth of these positions (scaled by initial thickness) with respect to time, approximately extending from y = 0 at t = 0, is consistent with the positions collapsing to fixed self-similar scaled (e.g., y/h) values. TheŨ 1 -based positions also evolve linearly and likewise collapse to fixed y/h values. Plotting the scaled positions of these points as a function of Atwood number (Figure 11b) highlights the prominent features observed in Figure 9: an increasing drift of the mean streamwise velocity profile to the light fluid side with increasing Atwood number, while the density profile remains approximately centered at the initial interface. In addition, Figure 11b indicates that the mean cross-stream velocityŨ 2 peak similarly drifts to the light fluid side, as well as the peak Reynolds stressR 12 ( §5.4). The relative magnitudes of the drifts confirm the predictions of the self-similar analysis ( §4.2) and are consistent with previous simulations of other variable-density mixing layers (e.g., Pantano & Sarkar 2002;Almagro et al. 2017). For the range of Atwood numbers simulated, η 12 < η 1 < η 2 < 0: the Reynolds stress peak is located further in the light fluid than the neutral point of mean streamwise velocity, which itself is further than the peak of mean cross-stream velocity.

Velocity Fluctuation Intensity Profiles
Statistical profiles for velocity fluctuations are similarly obtained using self-similar scaling applied to the y coordinate. It has also been verified that these profiles collapse over the self-similar growth time period (apart from a small amount of statistical variability) when scaled in this manner. These time-averaged profiles are compared among Atwood numbers in Figure 12.
Overall, the behaviors of the velocity variances for the low Atwood number case agree well with other published single-density mixing layer simulations. However, there can be significant differences in the magnitudes. The peak variance magnitudes of Rogers & Moser (1994) are 23% larger than those of the present A = 0.001 simulation. The peak magnitudes of the density ratio 1 simulation of Almagro et al. (2017) are on average 52% larger than those of the present simulation. The magnitudes for the Reynolds stressR 12 peak likewise differ between the simulations by similar amounts. The spatially developing mixing layer simulations of Attili & Bisetti (2012) that reach relatively high Reynolds numbers have peak magnitudes on average 19% greater than the present results.
One factor likely contributing to the differences of intensity magnitude is the determination of self-similar averaging time. With the present initial disturbance, an overshoot in the turbulence intensities occurs, and after a significant period of time the overshoot decays and asymptotes to the final self-similar growth state as the mixing layer thickens. Other simulations approach self-similar growth differently, and the self-similar period may be determined differently. Despite the difference of the spatial vs. temporally developing configuration, the Attili & Bisetti (2012) intensity profiles appear to agree most closely with the present simulation. Their simulation attains higher Reynolds number and greater thickness growth than the other temporal simulations cited. Differences in simulation domain sizes could potentially alter the turbulence dynamics by restricting structure growth and thereby affect fluctuation intensities. An additional factor may be persisting  Figure 12. Self-similar profiles for all Atwood numbers (colored as in Figure 3) showing velocity variances and Reynolds stresses. In (a-d), the solid line represents the velocity variance u i u j while the dashed line represents Favre-averaged Reynolds stressRij = ρu i u j /ρ. (e)R12 as in (d) but scaled by ∆U (dh/dt) as suggested by the self-similar analysis. (f) compares the Favre shear stress of (d) not scaled by local average density but Rij = ρu i u j scaled by the average of the free-stream densities ρ0.
effects of the differing initial disturbances. Among experiments, there is significant scatter in the intensity magnitudes, e.g., the differences between Bell & Mehta (1990) and Spencer & Jones (1971) as shown in Almagro et al. (2017). Rogers & Moser (1994) summarized the wide range of magnitudes for streamwise velocity variances measured in experiments (as well as the mixing layer growth rates, which are closely related to u 1 u 2 ). They also noted the perspective of Dimotakis & Brown (1976) that persisting influence of the initial conditions may be responsible. When Atwood number is increased, the behavior of the intensities and Reynolds stresses remain similar to the A = 0.001 case, except they shift to the light fluid side and generally decay slightly in magnitude. As shear moves to the light fluid side with increasing Atwood number, the turbulence intensity peak moves to the light fluid side as well. (The close relation between mean shear and the production of turbulent kinetic energyk =R ii /2 is apparent from the shear production term −ρR 12Ũ1,2 that dominates the budget forρk.) Velocity variances u i u j and Reynolds stressesR ij [ Figure 12(a-d)] both increasingly shift to the light fluid side with increasing Atwood number; this applies to the on-diagonal (i = j) elements as well as the streamwise-cross-stream (i = 1, j = 2). Figure 12(a-c) suggests that there is only a weak reduction in peak turbulent kinetic energy with increasing Atwood number. The reduction in peakR 12 Reynolds stress (or u 1 u 2 ) is as strong as that experienced by any of the on-diagonal turbulent kinetic energy contributions, yet it is reduced by no more than about 30% from A = 0.001 to A = 0.87. WhenR 12 Reynolds stress is scaled by ∆U and dh/dt as suggested by the self-similar analysis, rather than by ∆U 2 as is typically reported, the peak magnitudes weakly increase with increasing Atwood number (Figure 12e).
If Reynolds stress is scaled using the average density of the two free streams (ρ 0 ) rather than the local mean density, the reduction in peak value with Atwood number is enormous (Figure 12f). This is further confirmation that the intense turbulent motions move to (and are sustained in) light density fluid. u i u j andR ij = ρu i u j /ρ agree very closely for even the highest Atwood number throughout self-similar growth (while there are significant differences during transition with high A). This agreement is remarkable because these quantities do not agree well with R ij = ρR ij due to the shift of strong fluctuations to fluid on average lighter than ρ 0 . In other words, at elevated A,ρ is much smaller than ρ 0 at the position of peak turbulence intensity, but ρu i u j is also commensurately smaller so their ratio is nearly the same as for low A. Details of the local density distributions and how they correlate with velocity-based fluctuations will be further considered ( §6).

Analysis of Thickness Growth Rate During Self-Similar Growth
The statistical profiles discussed above can be related to growth rate attained by each mixing layer during its self-similar growth regime. The average growth rates calculated over the self-similar growth time intervals obtained above are first summarized as a function of Atwood number. At very low Atwood number (A = 0.001), the momentum thickness growth rate dδ m /dt/∆U = 0.012 agrees well with the value of 0.014 reported by Rogers & Moser (1994). Almagro et al. (2017) reports a somewhat higher growth rate of 0.017 when density is constant. In terms of thickness measured by h, the present simulation's growth rate dh/dt/∆U of 0.069 is consistent with the 0.062 value for the simulation of Rogers & Moser (1994), though both of these growth rates are toward the lower end of the 0.06-0.11 values typically observed in experiments (Pope 2000;Dimotakis 1991).
Assessing the growth rate reductions as a function of Atwood number across the present simulations, Figure 13(a) shows that the momentum thickness growth rate for A = 0.87 is reduced by 77% from the rate for the single-density case, while the momentum thickness per mass growth rate (δ m,pm ) and the analogous integral growth rate for the density profile (δ ρ ) experience lesser but nonetheless significant reductions. The stronger reduction for δ m can largely be explained by a misalignment that develops between density and velocity profiles ( §5.3). The reductions in growth rate based on h and h ρ (Figure 13b) are similar to those of the δ m,pm and δ ρ integral quantities ( Figure 13a); the reductions for h Xp are also similar (Figure 13d). The thicknesses derived from mean profile y-gradients (δ ω and δ dρ/dy ), however, display less smooth growth rate reduction behavior (Figure 13c). This could be a consequence of greater sensitivity to noise associated with a lack of statistical convergence when the maximum gradient is calculated on the smoothed profile, but could also appear due to the greater sensitivity of gradient-based measurements to the details of the profile asymmetries that appear with increased Atwood number. Though h Xp is an integral measurement and therefore lacks the extreme sensitivity to local noise in the mean profile of gradient-based measurements, it appears to display less smooth reductions in growth rate than the other integral thickness quantities. This suggests that some measurements may be particularly sensitive to specific features of the profile shapes. The close correspondence between most of the growth rates (particularly for δ m , δ m,pm , and h) confirms that any of the corresponding thickness measurements would be acceptable for scaling the flow statistics profiles during self-similar growth. Generally, the growths of density thickness quantities (due to mixing of fluids) also behave similarly to the growths of velocity thickness quantities.
To compare the growth rate effects of Atwood numbers for other variable-density mixing layers, Figure 14 Pantano & Sarkar (2002). The growth rates are scaled by the corresponding growth rate with A ≈ 0 for each data set. ible 0.7 convective Mach number). In contrast to the present species mixing governed by simplified INBM (incompressible non-Boussinesq mixing) equations, the latter cases with thermal variations are governed by the LMNOB (low-Mach number non-Oberbeck-Boussinesq) and fully compressible non-Oberbeck-Boussinesq (NOB) equations (Livescu 2020). A detailed comparison between simplified INBM equations and LMNOB equations has been made at A = 0.75 by Baltzer & Livescu (2020). Atwood number for the wide range shown in Figure 14 affects both momentum and vorticity thickness growth rates relatively similarly between density variation mechanisms (Figure 14), particularly given the differences between simulations in addition to species vs. thermal transport, e.g. domain sizes, initial disturbance, determination of self-similar growth period, etc. The growth rates are normalized by the A ≈ 0 growth rate for each configuration (simulation set), and doing so conceals important physical differences and their effects between cases: the growth rate of the A = 0 Pantano & Sarkar (2002) mixing layer had reduced by 40% relative to their low-speed simulation solely due to compressibility effects. Their simulations investigating the range of Atwood numbers are only available at a Mach number of 0.7. Since vorticity thickness is based on the maximum gradient of the mean streamwise velocity profile, more scatter in these growth rates is to be expected as they are sensitive both to the shape of the profile and noise in this quantity (smoothing was applied when calculating for the present simulations).
The influence of Atwood number on growth rate can be further explained based on the behavior of the statistical profiles. One useful property of the momentum thickness definitions (4.11-4.12) is that they straightforwardly lead to relations between growth rates and flow statistical quantities. This was explored by Vreman et al. (1996), who showed that an informative relation can be formed based on the time derivative of (4.11) that yields as other terms vanish using the y-integrated averaged continuity and momentum equations. Multiplying the Favre mean momentum equation (4.2) for i = 1 byŨ 1 produces an equation relating the time derivative of mean kinetic energy to Reynolds stress as Many terms are cast in conservative forms, so they vanish when integrated over the y domain in (5.4). Then, (5.4) becomes after neglecting the mean viscous term, which is consistent with the self-similar analysis ( §4.2) and has been shown to be small by scaling arguments of Pantano & Sarkar (2002). The growth rates calculated from this equation using the self-similar averaged profiles match those directly measured in the flow to within several percent. The above relation shows that the momentum thickness growth rate depends on the density, mean streamwise velocity, and Reynolds shear stress profiles. As shown in Figure 12(d), there is relatively little change in Favre-averaged Reynolds shear stress with the density normalized by the local mean. However, when normalized by the average of the two free stream densities ρ 0 , theρR 12 quantity that appears in (5.6) reduces strongly with increasing Atwood number, as shown in Figure 12(f). The simulation suite was designed such that the average of the free-stream densities ρ 0 remains consistent across all Atwood numbers. A consequence is that, if the fluids were fully mixed in equal proportion in the core of the mixing layer, the density there would be the same for every Atwood number. Therefore, reductions as observed in the aforementioned quantity normalized by ρ 0 are indicative ofR 12 having its strongest magnitude increasingly further in the light fluid. Furthermore, this density becomes smaller relative to ρ 0 as Atwood number increases. In (5.6), the growth rate is the product ofρR 12 with the mean streamwise velocity gradient, which weakly changes in magnitude with Atwood number (Figure 9b). This density weighting reflects the dependence on density in the momentum thickness definition (4.11). Thus, the principal cause of growth rate reduction for δ m is the turbulent motions and the associated momentum deficit shifting to the light fluid side.
The dominance of the profile shifting effect is demonstrated by artificially realigning the profiles such that the peak inR 12 and inflection point inŨ 1 are returned to the point whereρ = ρ 0 (Figure 15). Eliminating the shifts in this way significantly weakens the growth rate reduction effect. The magnitudes of growth rate reductions after these artificial shifts are approximately those observed for the growth rate of momentum thickness per mass (which lacks the density weighting). For other turbulent variabledensity mixing layers, Pantano & Sarkar (2002) and Almagro et al. (2017) have developed semi-empirical formulas that estimate momentum thickness growth rate reductions as functions of density ratio (or Atwood number) largely based on the shifts that develop between the mean streamwise and mean density profiles.
Equation (5.6) also shows that the mixing layer growth rate is the integral over the entire width of the mixing layer of a term that is closely related with the production of TKE (through the shear production mechanism, −ρR 12Ũ1,2 ). It is informative to split this integral into contributions from the light fluid and heavy fluid sides of the domain: where y ρ0 is the y value at whichρ(y) = ρ 0 . Since mean density increases monotonically across the interior of the mixing layer, it follows that the first term is the light fluid contribution to the growth rate and the second term is the heavy fluid contribution (in a mean sense). The y ρ0 position remains at y = 0 for A → 0, thus splitting the mixing layer in half. As Atwood number increases, the y ρ0 position moves much more weakly compared to points on the mean velocity profiles orR 12 (and y ρ0 instead moves toward the heavy fluid side), as shown by Figure 11(b). When the density differences are very weak (A = 0.001), the contributions are essentially equally split between the light and heavy sides (Figure 15a). The heavy-fluid growth rate contribution monotonically decays with Atwood number, which is consistent with the intense turbulence and momentum deficit (relative to the free streams) drifting to the light fluid side. As Atwood number increases from A = 0.001 to 0.25, the light-fluid growth rate contribution weakly increases, but then decays for higher Atwood numbers. These growth rate changes can be interpreted in light of the time histories of mean profile positions for various Atwood numbers ( Figure 11). Beginning from the symmetric growth of the mean streamwise profile edge positions for A = 0.001 in Figure 11(a), the edge positions reduce in their penetration to the heavy fluid side but increase in penetrating the light fluid side for Atwood numbers up to A = 0.75. The penetrations into the light fluid sides stagnate as Atwood number increases beyond 0.75, and by A = 0.87 the growth rate has decreased so much that even the growth into the light fluid has reduced below that for A = 0.75 while growth into the heavy fluid side is negligible. This reduction in . The hypothetical growth rate predicted by (5.6) if the profiles' drifts were artificially removed (· · · · · · )) highlights that much of the growth rate reduction is associated with the intense turbulence and shear shifting to the light fluid. (b) Momentum thickness per mass growth rate predicted from (5.11) (-•-). The mean shear-Reynolds stress term (--•--) has the same form as the production for TKE and dominates the other growth rate terms. However, variable-density terms also appreciably reduce the growth rate at high Atwood numbers, as shown by the distance between the curves. growth into the light fluid side is manifested in the sharply reducing light fluid growth contribution in Figure 15(a).
While much of the reduction in momentum thickness growth rate with increasing Atwood number can be explained by the velocity neutral point moving into the light fluid, weaker though significant reductions in growth rate also occur in all other thickness measurements, despite their lack of explicit density profile dependencies. This differs from the assumption of constant temporal growth rates with Atwood number sometimes used to develop theories addressing growth rate for spatially-developing variable-density mixing layers, as described in the introduction. To address these subtler reductions in growth rate with Atwood number, we derive an analogous equation relating the growth rate of momentum thickness per mass to statistical profiles of the flow. A similar derivation instead beginning from (4.12) leads to Again using the homogeneities of this flow and noting that bothŪ 2 and u 1 u j vanish on the boundaries (so the conservative terms integrate to 0) simplifies the expression to It can be shown that the terms based on the Reynolds mean cross-stream velocity decay as the flow evolves (and are of small magnitude, Figure 10). In contrast, the u j,j u 1 and p ,1 /ρ terms, as well as the dominantŪ 1,2 u 1 u 2 term, can be shown to maintain constant magnitudes with the appropriate self-similar scaling. Retaining only these terms, the relation simplifies to For this particular flow in which kinematic viscosity maintains a constant value everywhere, the mean viscous term simplifies to νŪ 1,22Ū1 . It decays with time and during the self-similar growth generates a very small effect, so it is neglected. Of the variable-density terms, p ,1 /ρ Ū 1 dy is negative and dominates in magnitude over − u j,j u 1 Ū 1 dy, which is positive, for all of the Atwood numbers considered. Note that in the single-density case, withŪ 2 = 0 and u j,j everywhere zero, this equation simplifies to dδ m,pm /dt = −2/∆U 2 Ū 1,2 u 1 u 2 dy, which is consistent with the usual momentum thickness growth equation (5.6).
This growth rate equation indicates that variable-density effects can modify the δ m,pm growth rate from its single-density value both through the additional terms (related to the density variations and corresponding divergence of the velocity field) and through changes in the mean velocity gradientŪ 1,2 and/or the Reynolds stress u 1 u 2 that appear in the mean shear-Reynolds shear stress term. Applying this equation to each Atwood number (Figure 15b) demonstrates that the mean shear-Reynolds shear stress term dominates the growth rate equation for all Atwood numbers and significantly decreases in magnitude for the highest Atwood numbers. To understand the reduction in growth rate associated with the dominantŪ 1,2 u 1 u 2 term, it can be seen that there is little change in peak U 1,2 value ( Figure 9) while peak u 1 u 2 magnitude decreases moderately ( Figure 12) with increasing Atwood number. A combination of lower peak stress magnitude and subtle changes in the mean streamwise velocity and stress profile shapes account for the reduction in growth rate. These phenomena also contribute to the conventional (densityweighted) momentum thickness growth rate reductions but are weaker than the effect of shifting mean streamwise velocity and density profiles.

Profiles Involving Density Fluctuations
Density-velocity correlations are contained in the normalized mass flux quantity a i = ρ u i /ρ. The turbulent mass fluxes also quantify the relationship between Favre and Reynolds averages for the velocity and Reynolds stress quantities considered above. ). As observed in Figure 10, the Favre average cross-stream velocityŨ 2 is much larger in magnitude than the Reynolds average cross-stream velocityŪ 2 . According to the above relations,Ũ 2 =Ū 2 + a 2 is dominated by the a 2 turbulent mass flux while the Reynolds mean is relatively insignificant and is explained by the relation developed in Appendix A. In singlefluid incompressible temporal mixing layers, the mean cross-stream velocity is zero, in order to satisfy the divergence-free condition. Thus, the correlations between the crossstream velocity and density in these variable-density mixing layers dominate the Favre mean cross-stream velocity. In addition, density fluctuation properties as revealed by fluctuation intensity are relevant to the structure of the flow.
Two types of normalizations for a i and density fluctuation intensity are considered. Thē ρ denominator included in the a i definition removes the density dimensional dependency, but a consequence is that a i generally grows with Atwood number as density fluctuations become more pronounced. Likewise, as Atwood number increases, density fluctuations increase in magnitude. No-mix, denoted by nm, corresponds to the quantities' values in a hypothetical configuration in which the two fluids are distributed without any mixing, and results in the highest possible magnitudes for ρ 2 . It can be shown that ρ 2 nm = (∆ρ/2) 2 = ρ 2 0 A 2 (Livescu & Ristorcelli 2008). By analogy, A∆U is an appropriate scaling for a i that is equivalent to scaling by [∆ρ/(2ρ 0 )]∆U . Figure 16(c-d) shows the self-similar profiles of (a-b) with these scalings applied.
In contrast to the velocity fluctuations' peak magnitudes shift to the light fluid side, the strongest density fluctuations position in the heavy fluid side with increasing Atwood number (Figures 16b,d). Fluctuation profiles for these density-based quantities (intensity and a i ) are also almost completely symmetric about y = 0 for the A = 0.001 case as the intense turbulence remains at the initial interface position. At increased Atwood number, large-scale disturbances (e.g., corrugations) form on the heavy fluid side while the fine scales of motion producing faster mixing are concentrated on the light fluid side. This behavior can be adduced from the mean density and velocity fluctuation intensity profiles and by density visualizations (Figure 17). The smoother yet disturbed heavy fluid side interface suggests the dominance of large-scale structures while small scales concentrate at the light fluid side for high Atwood number. The large displacements of the largely unmixed fluids relative to the background density gradient leads to large density fluctuation magnitudes on the heavy fluid side. In contrast, the greater mixing on the light fluid side produces local densities on average nearer to the mean densityρ, thus resulting in smaller density fluctuation intensities.
In absolute terms (scaled by ρ 0 in Figure 16b), the density fluctuation intensity magnitudes increase up to A = 0.75 but then stagnate for higher Atwood number. Peak intensity positions penetrate less and less deeply into the heavy fluid side for increasing Atwood numbers. These effects are likely a consequence of the less energetic turbulence sustained on the light fluid side reducing in ability to overcome the heavy fluid side's inertia and disturb it. Scaling using the mean density ρ 0 of the two streams, which remains constant for all Atwood numbers, reveals the relative magnitudes of the fluctuation intensities. The no-mix scaled density fluctuation intensity profiles (Figure 16d) increase in magnitude up to A = 0.75 but decay for higher Atwood number. No-mix intensity is proportional to the differences in densities between streams ∆ρ, so scaling by this quantity would be expected to scale out the effect of increasing density differences for Boussinesq mixing. Thus, magnitude differences under this scaling reveal non-Boussinesq effects in the fluctuation intensities.

Conditional Statistics
Conditional statistics expose correlations with local fluid density to further reveal variable-density effects in the flow. While the unconditional statistical moments above quantify the asymmetries (in a mean sense) with respect to y position, statistics conditioned on density reveal further asymmetries with respect to local density at fixed y. The unconditional statistics have demonstrated that increasing Atwood number concentrates the most intense turbulent motions at descending y positions of mean density progressively lighter than ρ 0 . TKE provides one indication of where turbulence is concentrated, while the dissipation term of its budget is associated with intense, small-scale turbulent motions. Enstrophy is another quantity associated with intense vortical motions. Turbulent kinetic energy dissipation per unit volume can be related to fluctuation enstrophy (based on vorticity ω k = ijk u j,i where ijk is the Levi-Civita symbol) as where d = u i,i is divergence (Morinishi et al. 2004). In constant-viscosity divergence-free incompressible flow, only the first term contributes on the right-hand side. Numerical simulations have indicated that the first term dominated while the third term made a small contribution and the others were negligible, although the detailed study was performed in a wall-bounded configuration (Morinishi et al. 2004). Therefore, a reasonable approximation to the relation between enstrophy and dissipation is whereε is defined byε = ε/ρ as with (5.3) and ν is constant within the present flow ( §3.3). The enstrophy and scaled dissipation agree well for both Atwood numbers shown in Figure 18(a-b). The peaks of enstrophy and dissipation are also shown to coincide with steep mean streamwise velocity gradients and, in the case of high Atwood numbers, to be located where mean density is significantly lower than the average of the two free-stream   densities (ρ 0 ). Since the definition of enstrophy is independent of density, it is a useful quantity for investigating density effects on the kinematics of turbulence.
The y plane of peak enstrophy intensity is an informative location for which to study the relationship between local density and intense turbulence. At this position, the enstrophy conditioned on density is plotted for each Atwood number in Figure 18(cd). The mean enstrophy is related to conditional counterpart by where f is the ρ probability density function (pdf), which indicates the frequencies with which fluid of a given density exists at this location. Thus, the peak enstrophy values shown in Figure 18(a-b) can be obtained from the conditional enstrophies shown in Figure 18(c-d) via this relation. The conditional enstrophy plots reveal that, at low Atwood number, the fluid carrying the greatest amounts of enstrophy has density equal to the mean density at this position, which coincides with the average of the two free streams' densities. Since no pure (free stream) fluid has that density, fluid of such density is a mixture of the pure fluids, and it can be concluded that fluid of this density is associated with strong mixing. As Atwood number becomes large, the mean density at the plane of strongest enstrophy is less than ρ 0 and the local fluid densities associated with the strongest enstrophy magnitudes have yet lower values according to the conditional statistics. Since mixing is associated with the intense small-scale motions, this suggests that the lighter fluid is carrying the turbulence and these motions are also active in mixing with the heavier fluid. Conversely, the larger scales of turbulent motion lack strong velocity gradients and are thus associated with weaker enstrophy; the conditional enstrophy suggests that the larger scales are thus associated with the denser fluid. The νh/∆U 3 scaling of the enstrophies shown in Figure 18 is based on the arguments given in Rogers & Moser (1994) using the relation between enstrophy and dissipation as well as the property that integrated dissipation remains constant in self-similar growth.
The pdfs also shown in Figure 18 demonstrate that the densities of peak conditional average are frequently present for both the lowest A and A = 0.75 cases. Thus, the densities associated with strong enstrophy magnitudes are representative of fluid parcels playing dominant roles and not merely rare events. The A = 0.001 (essentially passive scalar) case indicates that fluid near the local mean densityρ carries the strongest enstrophy per unit volume and is most prevalent, for the position of strongest enstrophy that occurs near y = 0. At high Atwood number, fluid lighter than the local mean density carries most of the enstrophy and is also most prevalent, at the peak total enstrophy position that has drifted toward the light-fluid side (y < 0) relative to the initial interface.
While the conditional enstrophy statistics provides quantitative evidence that vorticity is concentrated in the light fluid, flow visualizations are consistent with this result and illustrate the asymmetries present in the flow. Surfaces indicating scarcely-mixed (nearly free-stream density) fluid are shown in Figure 17. Since these surfaces were initially parallel planes, the topologies of the surfaces at subsequent times are consequences of the motions that transport the fluid. The surfaces are symmetric in appearance for the lowest Atwood number case (Figure 17a). For A = 0.75, however, the higher-density red surface (Figure 17c) is smoother with larger-scale corrugations compared to the lower-density blue surface (Figure 17d) that indicates the presence of much finer scales of motion. These fine scales are associated with strong enstrophy. The visualization for elevated Atwood numbers is thus consistent with the the average enstrophy peak existing in the vicinity of y values at which the density is lower than the average of the two streams.

Conclusions
The present set of shear-driven variable-density mixing layer DNS spans a wide range of Atwood numbers. Since these simulations have reached self-similar growth at sufficient Reynolds numbers to be past the mixing transition, they form a comprehensive data set for evaluating the variable density effects on late-time turbulence dynamics. The results demonstrate that, as Atwood number is increased while keeping the average density of the two free streams constant, the most intense turbulence is sustained in lighter-thanaverage fluid during self-similar growth. This occurs both in the sense of the intense turbulent motions shifting to y-positions at which the mean density is lower and also in the sense of the the strongest small-scale turbulent motions preferentially concentrating in fluid of lighter-than-mean density at a given position.
The main Atwood number effects on the most basic statistics, the mean density and velocity profiles, can be explained by self-similar growth properties and flow physics arguments. Self-similar analyses of the mean mass conservation and mean streamwise momentum balance equations have shown that the peak cross-stream velocity occurs in the light fluid side, while the neutral point of streamwise velocity moves further to the same side, and the peakR 12 stress moves yet further into the light fluid side ( §4.2). The intense turbulent motions occur where production of turbulence is concentrated, which is where the mean velocity profile is steepest andR 12 magnitude is large. Since the intense turbulent motions are also associated with mixing that smooths the density profile, the mean density profile becomes shallower near the strong small-scale turbulence regions, while thickness growth of both the mean density and mean streamwise velocity interfaces preferentially occurs on the light fluid side. The alignment of the peak enstrophy with the mean streamwise velocity and density profiles confirms this behavior (Figure 18). The drift of the velocity and Reynolds stress profiles to the light fluid side, as well as the asymmetry of shallower decay on the light fluid side of the mean density profiles, strengthens in degree as Atwood number increases. These effects are robust with respect to statistical noise and occurs in mixing layers with streams of differing density produced by a single fluid with varied thermodynamic properties, both in low-speed (Almagro et al. 2017) and compressible but subsonic (Pantano & Sarkar 2002) regimes. These profiles for incompressible multi-species and low-speed varied thermodynamic properties cases are directly compared by Baltzer & Livescu (2020).
A prominent variable-density effect is the reduction in growth rates as Atwood number increases. Using the commonly-used momentum thickness quantity to measure growth rate, this reduction has been shown to be primarily associated with the density weighting in the definition. This reflects the momentum deficit (relative to free-stream) growing in progressively lighter density fluid with increasing Atwood number; the thickness growth produces smaller changes in momentum as the density in which the growth occurs decreases. When the thickness is defined in a manner not weighted by the fluid density relative to the average of the pure-fluid densities (such as based on momentum on a permass basis or a quantity such as h that considers only the velocity field), the thickness growth rates display much less reduction as Atwood number increases. This is consistent with the mixing layer being sustained in lighter-than-average density fluid and, to a first approximation, behaving as a single-density mixing layer with a smaller nominal density (i.e.,ρ where the turbulence is strongest rather than ρ 0 ). For an actual single-fluid mixing layer, the growth rate is dependent only on the velocity difference across the free-streams regardless of the density. However, when conventional momentum thickness growth rate is calculated for variable-density mixing layers using ρ 0 as the density scaling in the definition, there is a mismatch between ρ 0 and the actual density in which the turbulence is sustained. Since turbulence is sustained in increasingly lighter fluid relative to ρ 0 with increasing Atwood number, this definition indicates a growth rate reduction with Atwood number. When this dependency is removed, as in the case of the other thickness measures, the growth rate reductions are more modest but nonetheless significant. These latter effects indicate variable-density induced departures from idealized single-density behavior. Such reductions are principally associated with decreases in u 1 u 2 magnitude, as indicated by (5.11).
The low-Atwood number limit of variable-density mixing layers captures much of the mass flux (density-velocity correlation) behavior, though the density field approaches passive scalar behavior. At high Atwood number, the mass fluxes become asymmetric and peaked on the light fluid side (particularly for the streamwise component), in addition to shifting to the light fluid side from the origin. Normalizing by only ∆U while ρ 0 remains constant, the turbulent mass flux magnitudes generally increase with Atwood number but appear to decrease past A = 0.75. This may be due to the greater heavy fluid inertia damping out the turbulence (as suggested by the Reynolds stress magnitudes at increasing Atwood numbers) despite the strengthening maximum density fluctuations.
The shifting of turbulent motions to the light fluid side influences the density field evolution. The intense turbulent motions progressively shifting toward the lighter fluid stream produces the asymmetry in the mean density profile discussed above. Furthermore, only larger spatial scales of velocity fluctuations are present toward the heavy fluid stream, which explains density fluctuation behavior: While the larger scales of motions near the heavier fluid are much less effective at producing mixing, they are effective at transporting parcels of largely unmixed heavy fluid, thereby producing particularly strong density fluctuations near the heavy fluid side at high Atwood number. These large density contrasts in partially-mixed light fluid and mostly-unmixed heavy fluid also produce large density fluctuation ( ρ 2 ) magnitudes there. Conditional statistics support the picture of turbulence being sustained within relatively light density fluid and penetrating into higher density regions where it is damped by the greater fluid inertia at high Atwood number.
The widespread nature of variable-density multi-fluid mixing motivates further advancements in properly capturing and modeling the variable-density effects on the kinematic structure of turbulence. This is particularly true given these effects' importance for predicting mixing and more complicated phenomena that closely depend on mixing, such as reactions, that appear in a wide range of flows. While many properties of variabledensity mixing layers closely resemble those of single-density mixing layers, complex interactions between density field and velocity field must be captured to predict the flow. In particular, the intense turbulence migrates to locally light fluid that interacts through neighboring fluid both through advection and mixing; this process alters the mean velocity and density profile evolutions as well as detailed statistics of mixing. Capturing all of these phenomena in a consistent manner presents significant challenges for RANS-and LES-type models.
This work has been authored by employees of Triad National Security, LLC which operates Los Alamos National Laboratory under Contract No. 89233218CNA000001 with the U.S. Department of Energy/National Nuclear Security Administration. Computational resources were provided at Los Alamos National Laboratory through the Institutional Computing (IC) program and at Lawrence Livermore National Laboratory through the Advanced Simulation and Computation (ASC) program.
Declaration of Interests. The authors report no conflict of interest.
Mean specific volumev was replaced using the identityv = (1 + b) /ρ. The complete expression is The mean cross-stream velocity is determined by only the mean density profile and the profiles of the density fluctuations (and their correlations); this result does not explicitly depend of the mean streamwise velocity. If there are no spatial density fluctuations besides the y-dependent mean density profile, only the first term on the right-hand side is nonzero and (A 3) reduces toŨ 2 = −D y ymin (lnρ) ,ii dy. This expression is valid when the simulation is initialized. If density is constant, (A 3) dictates thatŨ 2 = 0, which is consistent with the divergence-free nature of incompressible constant-density flow. Away from the interface in the variable-density mixing layers, each free stream has uniform (but different) density, so the equations indicate thatŨ 2 will be constant in these regions. Since each slip wall dictates thatŨ 2 will be zero at the wall while these equations indicate that U 2 will be constant approaching each wall, it is concluded thatŨ 2 is zero away from the region of mean density gradient or active mixing of unequal-density fluids.

Appendix B. Self-similar analysis
Self-similar analysis similar to that performed by Pantano & Sarkar (2002) is now given for the present flow configuration. If self-similar growth occurs, each self-similar profile quantity is a function only of η (for a given Atwood number/flow configuration) rather than position y and time t (or thickness) independently. Considering first the mean mass conservation equation (4.1) and substituting the self-similar variable dependencies yields − η[dh/dt]dρ/dη + d(ρŨ 2 )dη = 0.
(B 1) In order to ensure that the terms in this equation depend on η only, the mixing layer thickness needs to vary linearly in time. Thus, self-similarity requires that where C (the growth rate) is a constant. IfŨ 2 is scaled by C, equation (B 1) becomes independent of C, which indicates growth rate, rather than ∆U , as the self-similar scaling ofŨ 2 . Next, the self-similar variable forms are applied to the Favre mean streamwise momentum equation (4.2). When the flow is self-similar, the viscous term has a small effect on the mean momentum, as the mean velocity profile is relatively shallow after the earliest times (whereas the viscous term continues to produce large contributions to the energy balance because these contributions are based on instantaneous gradients within the turbulent motions). In self-similar variables, (4.2) then becomes − η[dh/dt]d(ρŨ 1 )/dη + d(ρŨ 1Ũ2 )/dη + d(ρR 12 )/dη = 0.
(B 3) Assuming that the streamwise mean velocity scale is ∆U , this equation also indicates that the self-similar scaling forR 12 is ∆U dh/dt, since for this scaling the equation becomes independent of the flow.
When both the mean density and velocity profiles are initially specified and centered at y = 0, the mean density profile is monotonic, with dρ/dη > 0. As the flow evolves, it continues to vary between the same density values for light and heavy fluid. The self-similar growth behavior implies that a non-monotonic density profile would be very unlikely to maintain. The behavior of the density profiles suggests choosing the density scale, ρ 0 , as the average of the two density extremes, (ρ 1 + ρ 2 )/2, which also corresponds to the initial centerline density.
Equation (B 7) shows that d(ρR 12 )/dη is also zero at η = η 2 . Let then η = η 12 be the location where dR 12 /dη is zero. Using the mean configuration set-up considered here (such that dÛ 1 /dη > 0 and dρ/dη > 0), similar arguments as above can be made to show that η 12 is unique within the layer and η 12 < η 2 . Thus,R 12 < 0 and has its largest magnitude further to the light fluid side thanÛ 2 . Analogous conclusions were previously drawn by Pantano & Sarkar (2002) from these equations for their configuration.