2.1. Governing equations

In this paper we primarily focus on what we call the single-species two-density RT model, which is governed by the fully compressible Navier–Stokes equations. In Appendix B, we repeat our main analysis on the more widely used two-species incompressible Navier–Stokes model (Cook & Dimotakis Reference Cook and Dimotakis2001; Livescu & Ristorcelli Reference Livescu and Ristorcelli2007) for comparison.

The dynamics is governed by the continuity equation (2.1), the momentum equation (2.2) and the total energy equation (2.3); the IE (2.4) and KE (2.5) equations are also shown for convenience as follows:

(2.1)\begin{gather} \partial_t\rho+\partial_j(\rho u_j)=0, \end{gather}

(2.2)\begin{gather}\partial_t(\rho u_i)+\partial_j(\rho u_i u_j)={-}\partial_i P+\partial_j\sigma_{ij}+\rho F_i, \end{gather}

(2.3)\begin{gather} \partial_t (\rho E) + \partial_j (\rho E u_j)={-}\partial_j (Pu_j)+\partial_j \left[2\mu u_i \left(S_{ij}-\frac{1}{d} S_{kk}\delta_{ij}\right)\right]\nonumber\\ -\partial_j q_j +\rho u_i F_i, \end{gather}

(2.4)\begin{gather} \partial_t (\rho e)+\partial_j(\rho e u_j)={-}P\partial_j u_j+2\mu\left(|S_{ij}|^{2}-\frac{1}{d}|S_{kk}|^{2}\right)-\partial_j q_j, \end{gather}

(2.5)\begin{gather} \partial_t\left(\rho \frac{|\boldsymbol{u}|^{2}}{2}\right)+\partial_j \left\{\left(\rho \frac{|\boldsymbol{u}|^{2}}{2}+P\right)u_j-2\mu\left(u_i S_{ij}-\frac{1}{d}u_j S_{kk}\right)\right\} \nonumber\\ =P\partial_j u_j-2\mu\left(|S_{ij}|^{2}-\frac{1}{d}|S_{kk}|^{2}\right)+\rho u_i F_i. \end{gather}

Here, $\rho$ is the density field, $\boldsymbol {u}$ is the velocity field, $P$ is pressure, $e$ is the specific IE (or IE per unit mass), $E=e+|\boldsymbol {u}|^{2}/2$ is total energy per unit mass. Here $S_{ij}=(\partial _i u_j+\partial _j u_i)/2$ is the strain rate tensor, while $\sigma _{ij}=2\mu (S_{ij}-d^{-1}S_{kk}\delta _{ij})$ is the deviatoric viscous stress tensor in $d$ space dimensions, where $\mu$ is dynamic viscosity. Here $\boldsymbol {F}$ is the external acceleration, $\boldsymbol {q}$ is the heat flux satisfying Fourier's law $\boldsymbol {q}=-\kappa \boldsymbol {\nabla } T$, where $\kappa$ is thermal conductivity and $T$ is temperature. The ideal gas law $P=\rho R T$ is adopted to close the set of equations, with $R$ the specific gas constant.

The governing equations are solved in 2-D and 3-D rectangular domains in Cartesian coordinates. The boundary conditions are periodic in the horizontal directions and no-slip rigid walls in the vertical $z$-direction. Following Bian et al. (Reference Bian, Aluie, Zhao, Zhang and Livescu2020), the density is initialized such that $\rho _0(z)$ is uniform in the $z$-direction on either side of the interface. Hydrostatic equilibrium requires that away from the interface, the initial (background) pressure varies as $P_0\sim \rho _0 g z$. For our single fluid case here, using the ideal gas equation of state yields that the background temperature gradient is constant and equal on both sides of the interface. Thus, the initial conditions represent a particular case of the analysis of Gerashchenko & Livescu (Reference Gerashchenko and Livescu2016), with $\textrm {d} T_0/\textrm {d}z= -g/R$. Away from the interface, when the thermal conductivity coefficient is constant as in our case here, a constant temperature gradient implies that the heat conduction term vanishes in the energy equation, so that the initial conditions are also in thermal equilibrium.

In this work, we choose viscosity $\mu$ to be spatially constant. In most practical applications, $\mu$ varies with temperature, sometimes strongly such as in ICF. Previous work by Aluie (Reference Aluie2013) and Zhao & Aluie (Reference Zhao and Aluie2018) has shown that the effects of a variable viscosity are still guaranteed to be localized to the smallest scales if the Reynolds number remains sufficiently high, thereby allowing for an inertial range of scales to exist, which is the main focus of this paper. We also choose the Prandtl number $Pr=c_p \mu /\kappa =1$, where again $\kappa$ is spatially constant, which is often not the case, especially in ICF (e.g. Zhang et al. Reference Zhang, Betti, Gopalaswamy, Yan and Aluie2018a). Our choices are similar to those adopted in previous studies of the fundamental hydrodynamics of RT (e.g. Cabot & Cook Reference Cabot and Cook2006; Reckinger, Livescu & Vasilyev Reference Reckinger, Livescu and Vasilyev2016; Wieland et al. Reference Wieland, Hamlington, Reckinger and Livescu2019). They are meant to simplify the problem to focus on the fundamental study of the cascades in RT flows, although we acknowledge ongoing research on the potentially important roles of temperature dependent viscosity and heat diffusion in some applications (e.g. Weber et al. Reference Weber, Clark, Cook, Busby and Robey2014).

2.2. Numerical method

The pseudospectral method (Patterson & Orszag Reference Patterson and Orszag1971) is adopted in the periodic directions, in which spatial derivatives and time integration are performed in Fourier space, while nonlinear terms are calculated in physical space. The fast Fourier transform provides an efficient way to convert the fields between physical and Fourier spaces. Aliasing errors in nonlinear terms are removed by the 2/3 dealiasing rule (Patterson & Orszag Reference Patterson and Orszag1971). For the vertical direction with wall boundaries, we use a compact finite difference scheme that is sixth order inside the domain and fifth order at the boundaries (Lele Reference Lele1992; Brady & Livescu Reference Brady and Livescu2019). In addition, the compact filtering scheme is performed in the inhomogeneous direction to filter out unphysical high-wavenumber oscillations (Lele Reference Lele1992; Brady & Livescu Reference Brady and Livescu2019). A fourth-order explicit Runge–Kutta method is used to integrate in time.

At the interface, as density changes between the heavy and light regions, the temperature gradient can no longer be constant. The energy equation then has non-zero time derivative at the initial time, $\partial _t (\rho e) = \partial _j (\kappa \partial _j T_0)$, which results in the generation of acoustic waves at initial time that can reflect off the top and bottom walls. We have carried out controlled tests by including a sponge layer at the top and bottom boundaries to damp these transient waves (Bian et al. Reference Bian, Aluie, Zhao, Zhang and Livescu2020), and could not detect any effect on the RT evolution. Even without the sponge layer, these initial acoustic waves are quickly damped at very early times, long before the flow becomes turbulent.

2.3. Analysis method

In fluid dynamics, coarse graining provides a natural and versatile framework to understand scale interactions (Leonard Reference Leonard1975; Meneveau & Katz Reference Meneveau and Katz2000; Eyink Reference Eyink2005). For any field $\boldsymbol {a}(\boldsymbol {x})$, a coarse-grained or (low-pass) filtered field, which contains modes at scales $\gtrsim \ell$, is defined in $n$-dimensions as

(2.6)\begin{equation} \bar{\boldsymbol{a}}_\ell(\boldsymbol{x}) = \int\,\textrm{d}^{n}\boldsymbol{r}\ G_\ell(\boldsymbol{r})\, \boldsymbol{a}(\boldsymbol{x}+\boldsymbol{r}), \end{equation}

where $G(\boldsymbol {r})$ is a normalized convolution kernel and $G_\ell (\boldsymbol {r})= \ell ^{-n} G(\boldsymbol {r}/\ell )$ is a dilated kernel with width comparable to $\ell$. The scale decomposition in (2.6) decomposes the field $\boldsymbol {a}(\boldsymbol {x})$ into a large-scale (${\gtrsim }\ell$) component $\bar {\boldsymbol {a}}_\ell$, and a small-scale (${\lesssim }\ell$) component, captured by the residual $\boldsymbol {a}'_\ell =\boldsymbol {a}-\bar {\boldsymbol {a}}_\ell$. The coarse-graining framework is a general method to analyse unsteady, anisotropic and inhomogeneous flows, such as the RT instabilities studied here, which conventional turbulence analysis techniques are unable to handle. More extensive discussions of the framework and its utility can be found in many references (Piomelli et al. Reference Piomelli, Cabot, Moin and Lee1991; Germano Reference Germano1992; Meneveau Reference Meneveau1994; Eyink Reference Eyink1995; Aluie & Eyink Reference Aluie and Eyink2009; Rivera, Aluie & Ecke Reference Rivera, Aluie and Ecke2014; Fang & Ouellette Reference Fang and Ouellette2016). In what follows, the subscript $\ell$ may be dropped from filtered variables when there is no risk of ambiguity.

Here, we use a Gaussian filter kernel at scale $\ell$ (Piomelli et al. Reference Piomelli, Cabot, Moin and Lee1991),

(2.7)\begin{equation} G_\ell({|\boldsymbol{x}|})=\left(\frac{6}{{\rm \pi} \ell^{2}}\right)^{n/2}\exp\left({-\frac{6}{\ell^{2}}|\boldsymbol{x}|^{2}}\right), \end{equation}

in dimensions $n=1,2,3$. Its Fourier transform is

(2.8)\begin{equation} \hat{G}_\ell({|\boldsymbol{k}|})=\exp\left({-\frac{\ell^{2}}{24}|\boldsymbol{k}|^{2}}\right). \end{equation}

Figure 1 plots the Gaussian kernel in $k$-space and compares it with sharp-spectral truncation. Note that a Fourier analysis is not possible in the non-homogeneous (vertical) direction. A detailed discussion on the different kernels can be found in Rivera et al. (Reference Rivera, Aluie and Ecke2014) (see figures 12 and 13 in that paper and associated discussion). The Gaussian kernel form in (2.7) has been used in several prior studies (Piomelli et al. Reference Piomelli, Cabot, Moin and Lee1991; Wang et al. Reference Wang, Wan, Chen and Chen2018) due to its numerical discretization advantages (see John Reference John2012, p. 30).

Figure 1. The Gaussian kernel $\hat {G}_\ell ({|\boldsymbol {k}|})$ in (2.8) compared with a sharp-spectral cutoff, ${H}_\ell ({|\boldsymbol {k}|})$, in Fourier space, with cutoff wavenumber $k_\ell =100$ and $\ell = L_z/k_\ell =6.4/100$. The horizontal axis has a logarithmic scale.

In simulations with non-periodic boundary conditions, such as our RT flow with rigid walls at the top and bottom, filtering near the boundary requires a choice for the fields beyond the boundary. In this paper, the choice is to extend the domain beyond the physical boundaries with values compatible with the boundary conditions (Aluie, Hecht & Vallis Reference Aluie, Hecht and Vallis2018). Filtering near the wall is then performed on this extended data. For our RT problem specifically, the velocity is kept at zero beyond the walls, the density field is kept constant (zero normal gradient), and the extended pressure field satisfies the hydrostatic condition $\textrm {d}P/\textrm {d}z=-\rho g$. The coarse-graining operation can then be performed at every point in the flow domain without requiring complicated inhomogeneous filters, which do not commute with spatial derivatives.

Analysing the energy cascade and inertial range dynamics in constant-density turbulence centres on studying the large-scale KE, $|\bar {\boldsymbol {u}}_\ell |^{2}/2$. Since KE in VD flows, $\rho |\boldsymbol {u}|^{2}/2$, is non-quadratic, its large-scale definition is not as straightforward. Several different definitions have been proposed in the literature, including $\bar {\rho }_\ell |\bar {\boldsymbol {u}}_\ell |^{2}/2$ (Chassaing Reference Chassaing1985; Bodony & Lele Reference Bodony and Lele2005; Burton Reference Burton2011; Karimi & Girimaji Reference Karimi and Girimaji2017) and $|\overline {(\sqrt {\rho }\boldsymbol {u})}_\ell |^{2}/2$ (Kida & Orszag Reference Kida and Orszag1990; Cook & Zhou Reference Cook and Zhou2002; Wang et al. Reference Wang, Wan, Chen and Chen2018). However, as was shown in Zhao & Aluie (Reference Zhao and Aluie2018), the above two definitions can fail to satisfy the ‘inviscid criterion’ in flows with large density variations, leading to viscous contamination at all scales, and preventing the unravelling of inertial range dynamics. The ‘inviscid criterion’ requires that viscous contributions be negligible at large scales. Aluie (Reference Aluie2013) proved mathematically, and Zhao & Aluie (Reference Zhao and Aluie2018) showed numerically, that the Favre decomposition, $(\bar {\rho }_\ell |\tilde {\boldsymbol {u}}_\ell |^{2})/2$, where

(2.9)\begin{equation} \tilde{\boldsymbol{u}}=\overline{\rho\boldsymbol{u}}_\ell/\bar{\rho}_\ell \end{equation}

is the Favre filtered velocity (Favre, Gaviglio & Dumas Reference Favre, Gaviglio and Dumas1958; Aluie Reference Aluie2013; Eyink & Drivas Reference Eyink and Drivas2018), satisfies the inviscid criterion even in the presence of arbitrarily large density variations. In the current paper, we shall adopt the Favre decomposition when analysing scale dynamics.

2.4. Simulation results

Two-dimensional and 3-D RT simulations are performed to analyse the energy budget and cascade in RT turbulence. These simulations, summarized in table 1, are performed at low and high Atwood numbers, with grid resolution up to $4096\times 8192$ in two dimensions and $1024\times 1024\times 2048$ in three dimensions. The Atwood number $\mathcal {A}$ is defined by $\mathcal {A}=(\rho _h-\rho _l)/(\rho _h+\rho _l)$, where $\rho _h, \rho _l$ are the densities of the initial uniform heavy and light fluids. The results shown below are primarily based on the largest simulations 2D4096 and 3D1024, while other simulations shown in the table are for numerical convergence and verification purposes.

Table 1. Parameters of the 2-D RT and 3-D RT simulations conducted in this paper. All simulations have a spatially uniform dynamic viscosity $\mu$ and Prandtl number $Pr=c_p \mu /\kappa =1$, with thermal conductivity $\kappa$ and specific heat capacity $c_p$ at constant pressure. Here $L_x,L_y,L_z$ are the domain lengths in three directions, and $\mathcal {A}$ is the Atwood number. Gravitational acceleration $g=1$ is in the $-z$ direction. The mesh Grashof number is $Gr=2\mathcal {A}g\langle \rho \rangle ^{2} \Delta x^{3}/\mu ^{2}$, the Kolmogorov scale is $\eta =\mu ^{3/4}/(\epsilon ^{1/4}\langle \rho \rangle ^{3/4})$ and the Reynolds number is $Re=\langle |\boldsymbol {u}|^{2}\rangle ^{1/2} L_x\langle \rho \rangle /\mu$, where $\langle \cdot \rangle$ denotes the spatial mean value. Here $Ma=\sqrt {g L_x}/\sqrt {\gamma R T_0}$ denotes the interface Mach number (Reckinger et al. Reference Reckinger, Livescu and Vasilyev2016), where $R$ is the specific gas constant, $\gamma$ is the ratio of specific heats and $T_0$ is the initial temperature at the interface. In the table, the Reynolds number and the Kolmogorov scale are calculated at dimensionless time $\hat {t}=t/\sqrt {{L_x}/{\mathcal {A}g}}=4.0$.

The initial conditions for simulations shown in table 1 are as follows. The initial density field is composed of uniform heavy fluid on top with $\rho _h=1$, and uniform light fluid with density $\rho _l$ at the bottom. The initial pressure field satisfies hydrostatic balance $\textrm {d}P / \textrm {d}z = -\rho g$, with pressure at the bottom set to be 10 for all simulations in table 1, except 2D1024lowM, where the pressure is set to 100 at the bottom to yield lower compressibility levels. In Appendix A, we show that our results do not depend on the compressibility level of RT flows. The velocity fields are initially zero with perturbations added to the interface of vertical velocity. The interfacial velocity perturbations are imposed in wavenumber space, and are non-zero only within a band of high wavenumbers. In three dimensions, the velocity amplitude is perturbed in a shell $k\in [32, 128]$, and the magnitude is set to be proportional to $\exp ({-({1}/{c})|k_x^{2}+k_y^{2}-80^{2}|})$, where $c = 22.63$ is a normalization constant to further limit the range of effectively perturbed wavenumbers. While in two dimensions, the velocity amplitude is perturbed in an annulus $k\in [72, 128]$, with magnitude proportional to $\exp ({-({1}/{c})|k_x^{2}-100^{2}|})$, where $c$ is the same constant as in three dimensions. Since the perturbation wavenumber band is slightly different between 2D4096 and 3D1024, for verification purposes, we also perform additional 2-D simulations 2D1024 and 2D2048 (see table 1) with initial perturbations identical to those in 3D1024. Analysis results based on these 2-D simulations are similar to the 2D4096 case, and are shown in the supplementary material, section D.

Visualizations of density snapshots are shown in figure 2(a) for the 2D4096 case and figure 2(b) for the 3D1024 case, both at dimensionless time $\hat {t} = t/\sqrt {{L_x}/{\mathcal {A}g}}= 4.0$, where $\sqrt {{L_x}/{\mathcal {A}g}}$ is the characteristic time scale and $g$ is the gravitational acceleration. Both figures show complex flow patterns and structures with a wide spectrum of sizes. It is also easy to see that the 3-D field develops much finer scale structures than the 2-D field, a property which we will quantify using spectra below. We now present some basic properties of our simulations.

Figure 2. Density visualizations from the simulations 2D4096 and 3D1024 in table 1. The dimensionless time for both is $\hat {t} = t/\sqrt {{L_x}/{\mathcal {A}g}}= 4.0$. The mixing width, which is defined and analysed in § 2.5, is $h\approx 1.69$ for 2-D RT, and $h\approx 1.1$ for 3-D RT. The density fields are turbulent with the presence of a wide range of scales. (a) This is the 2D4096 snapshot and (b) is the 3D1024 snapshot.

2.5. Mixing width

The mixing width $h(t)$ is an important quantity characterizing the total width of the mixing layer between bubble and spike fronts in RT flows. It is used to quantify the penetration depth of RT bubbles and spikes into the unperturbed flow. When there are sufficient multimode interactions, turbulent RT flows become self-similar at later times (Birkhoff Reference Birkhoff1955; Dimonte et al. Reference Dimonte2004; Ristorcelli & Clark Reference Ristorcelli and Clark2004). During the self-similar stage, the mixing width is well-described by $h(t)=\alpha \mathcal {A} g t^{2}$, a quadratic growth in time which is proportional to a constant coefficient $\alpha$.

The value for $\alpha$ has been the subject of vigorous research and discussion within the community. It is believed to depend on whether initial conditions render the RT flow in the bubble competition regime or in the bubble merger regime (e.g. Ramaprabhu et al. Reference Ramaprabhu, Dimonte and Andrews2005). Numerical simulations of 3-D RT in the bubble merger regime generally report values $\alpha \approx 0.02$ to $0.03$ (see table 6.1 in Zhou (Reference Zhou2017a) for a summary), although values larger than $0.05$ have been reported in 3-D simulations using front tracking (e.g. Glimm et al. Reference Glimm, Sharp, Kaman and Lim2013). Experiments, on the other hand, generally report much larger values $\alpha \approx 0.05$ to $0.06$ (Zhou Reference Zhou2017a), which is thought to be a result of unintended large-scale perturbations (Dimonte et al. Reference Dimonte2004).

While we report $\alpha$ values from our simulations here, its precise value is far from being the focus of our paper. Our interest in $\alpha$ here is rather more limited to comparing its values from 2-D simulations relative to 3-D simulations within the same modelling framework. For reasons that remain unclear and subject to various speculations, the parameter $\alpha$ in 2-D simulations has often been found larger than in 3-D simulations (e.g. Youngs Reference Youngs1991,  Reference Youngs1994; Dimonte et al. Reference Dimonte2004; Cabot Reference Cabot2006). Such enhanced RT growth in 2-D simulations has been somewhat of a puzzle given that both analytical models and simulations of single-mode RT predict that it is growth in three dimensions that should be faster (Layzer Reference Layzer1955; Goncharov Reference Goncharov2002; Bian et al. Reference Bian, Aluie, Zhao, Zhang and Livescu2020). A goal of this paper is to shed light on the puzzle.

We calculate the mixing width using the formula (Cabot & Cook Reference Cabot and Cook2006; Zhang et al. Reference Zhang, Betti, Yan, Zhao, Shvarts and Aluie2018b),

(2.10)\begin{equation} h(t)=\frac{2}{\rho_h-\rho_l}\int_{-\infty}^{\infty} \min\left(\langle \rho\rangle_{xy}(z)-\rho_l, \rho_h-\langle \rho\rangle_{xy}(z)\right)\,\textrm{d}z, \end{equation}

where $\langle \rho \rangle _{xy}(z)$ denotes a horizontal average of the density field, which is a function of $z$ and $t$, and $\rho _h, \rho _l$ denote the original heavy and light fluid densities. After obtaining $h(t)$, $\alpha$ can be calculated using one of two methods. The first is by directly calculating $h(t)/\mathcal {A}gt^{2}$, while the second is by carrying a best linear fit between $\sqrt {h(t)}$ and $\sqrt {\mathcal {A}g}t$ in the self-similar regime (Olson & Jacobs Reference Olson and Jacobs2009). Figure 3 plots $\sqrt {h(t)}$ versus time from the 2D4096 and 3D1024 data, and a linear fit in the self-similar regime is marked with a dashed line. The inset of the figures corresponds to the first approach: $\alpha =h/\mathcal {A}gt^{2}$ is plotted against time, while the $\alpha$ value obtained by the linear fit is marked with a horizontal dashed line. We obtain the same results from these two methods. Figure 3 shows that there is a substantial time range of self-similarity in both two dimensions and three dimensions, in which there is an excellent linear fit between $\sqrt {h(t)}$ and $t$. The values of $\alpha$ we obtain (0.036 in two dimensions and 0.021 in three dimensions) fall within the range reported from simulations in the literature (Dimonte et al. Reference Dimonte2004; Cabot & Cook Reference Cabot and Cook2006). Note that $\alpha$ from the 2-D RT is almost twice as large as that from 3-D RT.

Figure 3. Evolution of the square root of mixing width $\sqrt {h(t)}$ versus time in two dimensions and three dimensions, for which self-similarity analysis predicts a linear relation at late time (Olson & Jacobs Reference Olson and Jacobs2009). The slope of the straight line is 0.135 in (a) and 0.104 in (b), which corresponds to $\alpha =0.036$ and $\alpha =0.021$, respectively. The ‘+’ markers in both panels correspond to dimensionless time $\hat {t}=t/\sqrt {{L_x}/{\mathcal {A}g}}= 2, 3, 4$. The inset gives the compensated plot $\alpha =h(t)/(\mathcal {A}gt^{2})$ versus time, in which the horizontal lines correspond to the $\alpha$ values obtained by the linear fit. (a) The 2-D mixing width and (b) the 3-D mixing width.

2.6. Overall energy balance

We shall now discuss the bulk balance between PE, KE and IE in our simulated flows. Temporal evolution of the spatial mean of each of the budget terms in the KE equation (2.5) is shown in figures 4(a) and 4(c) for two dimensions and three dimensions, where $D=S_{ij} (2\mu S_{ij}-({2}/{d}) \mu S_{kk} \delta _{ij})$ is the dissipation term, and $\epsilon ^{{inj}}=-\rho g u_z$ is the energy injection. Correspondingly, figures 4(b) and 4(d) show the temporal evolution of bulk (over volume $V$ of the domain) PE $\delta \text {PE}$, bulk IE $\delta \text {IE}$ and bulk KE $\delta \text {KE}$, in two dimensions and three dimensions. The bulk energy is measured relative to its initial value:

(2.11)\begin{gather} \delta \text{PE}(t)=\int [\rho(\boldsymbol{x},t)-\rho(\boldsymbol{x},0)]gz\,\textrm{d}V={-}\int_0^{t}\,\textrm{d}t \int\,\textrm{d}V \epsilon^{{inj}}, \end{gather}

(2.12)\begin{gather}\delta \text{IE}(t) = \int [\rho e(\boldsymbol{x},t)-\rho e(\boldsymbol{x},0)]\,\textrm{d}V=\int_0^{t}\,\textrm{d}t \int\,\textrm{d}V (D-P\boldsymbol{\nabla} \boldsymbol{\cdot} \boldsymbol{u}), \end{gather}

(2.13)\begin{gather}\delta \text{KE}(t) =\int \left[\tfrac{1}{2}\rho \boldsymbol{u}^{2}(\boldsymbol{x},t)-\tfrac{1}{2}\rho \boldsymbol{u}^{2}(\boldsymbol{x},0)\right]\,\textrm{d}V=\int_0^{t}\,\textrm{d}t\int\,\textrm{d}V \partial_t \left(\tfrac{1}{2}\rho |\boldsymbol{u}|^{2}\right). \end{gather}

The second equality in (2.11) is due to the relation $\int _V \rho g u_z\,\textrm {d}V = ({\textrm {d}}/{\textrm {d}t})\int _V \rho g z\,\textrm {d}V$ (Cabot & Cook Reference Cabot and Cook2006), while (2.12) can be obtained by time integrating the IE equation (2.4).

Figure 4. Temporal evolution of terms in the KE budget, and overall energy balance in two dimensions (a,b) and three dimensions (c,d). Here $D=-\sigma _{ij} S_{ij}$ is the KE dissipation rate, $\epsilon ^{{inj}}=\rho u_i g_i$ is the energy injection rate. Time $\hat {t}= t/\sqrt {{L_x}/{\mathcal {A}g}}$ in the panels is dimensionless.

Plot of 2-D instantaneous energy budget terms in figure 4(a) shows the time derivative of KE and the energy injection rate keep increasing in time, while the average value of $P\boldsymbol {\nabla }\boldsymbol {\cdot } \boldsymbol {u}$ is positive so that IE is converted to KE via this term. The overall viscous dissipation term $-\langle D\rangle$ is negative and converts KE to IE, but the amount of conversion is smaller in magnitude relative to $\langle P \boldsymbol {\nabla }\boldsymbol {\cdot } \boldsymbol {u}\rangle$ in two dimensions. Thus, in figure 4(b) both the overall IE and PE decrease and are converted to KE. The overall bulk energy balance satisfies $\delta \text {PE}+\delta \text {KE} + \delta \text {IE} = 0$, as is marked by the ‘residue’ term in the figure.

The corresponding 3-D results are shown in figures 4(c) and 4(d). The injection rate and time derivative of KE keeps increasing, while the dissipation and pressure dilatation terms convert between KE and IE. One difference from two dimensions is that in three dimensions, the dissipation term contributes more than $P\boldsymbol {\nabla }\boldsymbol {\cdot } \boldsymbol {u}$, and the overall energy conversion is from KE to IE, opposite to the 2-D case. Total energy conservation, $\delta \text {PE}+\delta \text {KE} + \delta \text {IE} = 0$, still holds in three dimensions. This conservation holds in all simulations in table 1. In § 3, we will present a more refined balance relation applied to the coarse-grained quantities, where the balance of energy at different length scales is established.

The rate-governing quantities in our figure 4 seem more sensitive to fluctuations around the self-similar trajectory than the mixing width $h(t)$ in figure 3, which has often been the main focus of previous RT studies. The fluctuations may be either due to intrinsic variability or due to a transient response to initial conditions, the removal of which has been recently discussed in Horne & Lawrie (Reference Horne and Lawrie2020). In Appendix A for the high interface Mach number case, the corresponding figure shows similar fluctuations under different initial and compressibility conditions suggests it may be the former, yet a careful investigation of this observation is beyond the scope of our paper. Nevertheless, we shall show below that energy transfer across scales exhibits clear self-similar evolution despite the presence of these fluctuations in the rate-governing quantities.

One may notice that mean injection, $\langle \epsilon ^{{inj}}\rangle \equiv -\langle \rho g u_z\rangle$, shown in figures 4(a) and 4(c), grows more rapidly in two dimensions than in three dimensions. There are two (related) ways to understand this. (i) The mixing width growth rate $\alpha$ in two dimensions is larger than in three dimensions and, thus, the mixing layer width should be larger in two dimensions than in three dimensions at any physical time. This implies that relatively more PE is released in two dimensions due to the rising light fluid and sinking heavy fluid, leading to a higher mean injection $\langle \epsilon ^{{inj}}\rangle$. (ii) From the balance relation $\delta \text {PE}+\delta \text {KE} + \delta \text {IE} = 0$, we have $\delta \text {IE}<0$ in two dimensions, which is converted to KE, while $\delta \text {IE}>0$ in three dimensions. This implies a larger mean $\delta \text {KE}$ is expected in two dimensions, and thus $u_z$ is expected to be larger in two dimensions than in three dimensions. This leads to a larger $\langle \epsilon ^{{inj}}\rangle = -\langle \rho g u_z\rangle$ in two dimensions. We shall show below that there is a deeper physical explanation behind these differences between 2-D RT and 3-D RT arising from differences in the energy pathways across scales. Given the difference in $\langle \epsilon ^{{inj}}\rangle$ between the RT flows in two dimensions and three dimensions, we shall normalize by it to allow for a fair comparison of the budget terms between two dimensions and three dimensions.

2.7. ‘Filtering’ spectra

For completeness, we now present the spectra of density, velocity and KE using our RT data. In homogeneous constant-density flows, the energy spectrum is obtained by taking the Fourier transform of the velocity field and summing over shells in wavenumber space: $E(k) = ({1}/{2})\varSigma _{k-0.5 < |\boldsymbol {k}| \leq k+0.5}|\hat {u}(\boldsymbol {k})|^{2}$. Obtaining the spectra for variable density RT flows is difficult because (i) the domain is not periodic in the vertical direction so the Fourier transform is not straightforward, and (ii) the non-quadratic nonlinearity of KE poses another problem.

Previous works in RT (Cook & Zhou Reference Cook and Zhou2002; Cabot et al. Reference Cabot, Schilling and Zhou2004) obtained the spectra for density and velocity by Fourier transforming in the horizontal periodic directions and then taking an average in the vertical direction. By taking Fourier transforms only in the horizontal directions, spectral information in the vertical direction is completely ignored.

Unlike in constant density turbulence, where the spectrum of energy is the same as that of velocity, in VD the two quantities are different. Moreover, when KE is non-quadratic, there are many possible ways to construct a ‘spectrum’ in $k$-space having units of energy per wavenumber. These different ‘energy spectra’ correspond to the Fourier transform (or scale decomposition) of different combinations of variables and are not equivalent. To avoid the difficulty pertaining to its non-quadratic nature, previous studies constructed a KE spectrum by introducing a new variable $w_i=\sqrt {\rho }u_i$ and analysing it in Fourier space, $\hat {w}_i(k_x,k_y)$. However, as emphasized in Zhao & Aluie (Reference Zhao and Aluie2018), the evolution of $w_i$ is fundamentally different from the evolution of momentum, $(\rho u_i)$, especially when each of these variables is being decomposed in scale. A length scale $\sim 1/k$ does not exist on its own but is associated with the variable being decomposed. Zhao & Aluie (Reference Zhao and Aluie2018) showed that this new variable $w_i=\sqrt {\rho }u_i$ fails to satisfy the ‘inviscid criterion’ in VD flows, which means that dissipation can contaminate its large-scale dynamics in the presence of significant density variations. Thus, $\sqrt {\rho }u_i$ is not an appropriate quantity to analyse because it is inconsistent with inertial-range dynamics. Moreover, as was observed in Zhao & Aluie (Reference Zhao and Aluie2018), a putative power-law scaling of a quantity alone, such as $|\widehat {\sqrt {\rho }\boldsymbol {u}}|^{2}(k)$, is not indicative of an inertial cascade (see supplementary material, section C).

Here we adopt a new approach to obtain the spectra using spatial coarse graining. For any field $\boldsymbol {u}$, the filtering spectrum is defined as (Sadek & Aluie Reference Sadek and Aluie2018)

(2.14)\begin{equation} E_{\boldsymbol{u}}(k_\ell)\equiv \frac{\textrm{d}}{\textrm{d}k_\ell}\langle |\bar{\boldsymbol{u}}_\ell(\boldsymbol{x})|^{2}\rangle /2, \end{equation}

where $k_\ell =L/\ell$, $L$ is the domain size of interest, $\ell$ is the scale we are probing and $\langle \cdot \rangle$ stands for spatial averaging. The filtering spectrum of KE is

(2.15)\begin{equation} E_{{KE}}(k_{\ell})\equiv \frac{\textrm{d}}{\textrm{d}k_\ell}\langle \bar{\rho}_\ell|\tilde{\boldsymbol{u}}_\ell(\boldsymbol{x})|^{2}\rangle /2. \end{equation}

Sadek & Aluie (Reference Sadek and Aluie2018) showed that this filtering spectrum guarantees energy conservation, similar to the Plancherel theorem when working with Fourier spectra. The filtering spectrum allows us to measure the energy content at scale $\ell$ by probing the change in $\langle \bar {\rho }_\ell |\tilde {\boldsymbol {u}}_\ell (\boldsymbol {x})|^{2}\rangle /2$ as $\ell$ is varied. The calculation can be done in physical space without any Fourier transforms. The obvious advantage is that definition (2.14) allows us to extract the spectrum in inhomogeneous flows such as RT. Moreover, the filtering spectrum generalizes easily to non-quadratic quantities, such as KE, $\rho \boldsymbol {u}^{2}/2$, in VD flows. The Fourier spectrum, on the other hand, is limited to quadratic quantities and requires treating KE as such, even when this is inconsistent with the dynamics.

Figure 5 shows the spectra of density, velocity and KE for 2-D RT and 3-D RT flows calculated with the above definition. The density spectra follow a putative $k_\ell ^{-5/3}$ scaling over an intermediate wavenumber range in both two dimensions and three dimensions. The velocity spectrum in two dimensions is steeper than $k_\ell ^{-5/3}$ but shallower than $k_\ell ^{-3}$ as would have been expected from 2-D homogeneous constant-density turbulence (e.g. Boffetta & Ecke Reference Boffetta and Ecke2012). In three dimensions the filtering spectrum is slightly shallower than a $k_\ell ^{-5/3}$ scaling over a decade. The KE spectra follow similar scalings as the velocity field. Both velocity and KE spectra attain a peak value at low wavenumbers, where most of the energy content resides. In contrast, the density spectra show the importance of much larger scales due to density variation in the vertical. Indeed, the density jumps between heavy and light fluids are expected to yield a spectrum ${\sim }k_\ell ^{-2}$, which is consistent with the scaling at the largest scales. Plots in figure 5 demonstrate the utility of the coarse-graining approach to probe the spectral content in all directions, including scales along the inhomogeneous (vertical) direction, as discussed in Sadek & Aluie (Reference Sadek and Aluie2018). For completeness, we show in the supplementary material section C the filtering spectra for three different scale decompositions of KE.

Figure 5. Density, velocity and KE ‘filtering’ spectra versus filtering wavenumber $k_\ell =L_z/\ell$ (horizontal axes) for 2-D RT and 3-D RT turbulence at dimensionless time $\hat {t}=4$, visualized in figure 2. The spectra are calculated using a Gaussian filter. A $k_\ell ^{-5/3}$ scaling is shown for reference, and for 2-D velocity and KE spectra an additional $k_\ell ^{-3}$ scaling is shown. Note that $k_\ell$ is not a Fourier wavenumber, just a proxy for scale. The panels show: (a) density spectrum (2-D); (b) velocity spectrum (2-D); (c) KE spectrum (2-D); (d) density spectrum (3-D); (e) velocity spectrum (3-D); (f) KE spectrum (3-D).

3.1. Coarse-grained energy budget equations

Within the Favre decomposition, the budgets for large- and small-scale KE, IE and PE are

(3.1)\begin{gather} \partial_t\bar{\rho}_\ell \frac{|\tilde{\boldsymbol{u}}_\ell |^{2}}{2}+\boldsymbol{\nabla}\boldsymbol{\cdot} \boldsymbol{J}_\ell ={-}\varPi_\ell-\varLambda_\ell+\bar{P}_\ell \boldsymbol{\nabla}\boldsymbol{\cdot}\bar{\boldsymbol{u}}_\ell-D_\ell+\epsilon_\ell^{{inj}}, \end{gather}

(3.2)\begin{gather}\partial_t \frac{\bar{\rho}_\ell \tilde{\tau}_\ell(u_i, u_i)}{2} + \boldsymbol{\nabla} \boldsymbol{\cdot} \boldsymbol{J}_\ell^{{small}} = \varPi_\ell +\varLambda_\ell +\bar{\tau}_\ell(P,\boldsymbol{\nabla}\boldsymbol{\cdot} \boldsymbol{u})-D_\ell^{{small}}+\epsilon_\ell^{{small}}, \end{gather}

(3.3)\begin{gather}\partial_t \overline{\rho e}_\ell + \boldsymbol{\nabla} \boldsymbol{\cdot} [\overline{\rho e \boldsymbol{u}} -\overline{\kappa \boldsymbol{\nabla} T}] ={-}\bar{P}_\ell \boldsymbol{\nabla} \boldsymbol{\cdot} \bar{\boldsymbol{u}}_\ell - \bar{\tau}_\ell(P, \boldsymbol{\nabla}\boldsymbol{\cdot} \boldsymbol{u})+D_\ell^{{int}}, \end{gather}

(3.4)\begin{gather}\partial_t \bar{\mathcal{P}}_\ell + \boldsymbol{\nabla}\boldsymbol{\cdot} \overline{\mathcal{P} \boldsymbol{u}} ={-}\epsilon_\ell^{{inj}}, \end{gather}

where $e$ is specific IE, $\mathcal {P}=\rho g z$ is PE density, $\boldsymbol {J}_\ell$ is spatial transport of KE, $\varPi _\ell$ and $\varLambda _\ell$ are KE fluxes across scale $\ell$, $-\bar {P}_\ell \boldsymbol {\nabla } \boldsymbol {\cdot } \bar {\boldsymbol {u}}_\ell$ is large-scale pressure dilatation, $D_\ell$ is viscous dissipation and $\epsilon _\ell ^{{inj}}$ is KE injection (hence, the ‘inj’ superscript) due to gravity at scales larger than $\ell$. These terms are defined as

(3.5)\begin{align} \left. \begin{aligned} \varPi_\ell(\boldsymbol{x}) & ={-}\bar{\rho} \partial_j \tilde{u}_i \tilde{\tau}(u_i,u_j); \ \varLambda_\ell(\boldsymbol{x})=\frac{1}{\bar{\rho}}\partial_j \bar{P} \bar{\tau}(\rho,u_j);\ \epsilon_\ell^{{inj}}(\boldsymbol{x})=\tilde{u}_i\bar{\rho}\tilde{g}_i; \\ D_\ell(\boldsymbol{x}) & =\partial_j \tilde{u}_i \left[2\overline{\mu S_{ij}}-\frac{2}{d}\overline{\mu S_{kk}}\delta_{ij}\right]; \ J_j(\boldsymbol{x})=\bar{\rho} \frac{|\tilde{\boldsymbol{u}}|^{2}}{2}\tilde{u}_j+\bar{P}\bar{u}_j+ \tilde{u}_i \bar{\rho} \tilde{\tau}(u_i,u_j)-\tilde{u}_i\bar{\sigma}_{ij}; \\ D_\ell^{{small}}(\boldsymbol{x}) & =\overline{\partial_j u_i \sigma_{ij}}-\partial_j \tilde{u}_i \overline{\sigma}_{ij}; \ \epsilon_\ell^{{small}}=\bar{\rho}\tilde{\tau}(u_i,g_i); \ D_\ell^{{int}}=\overline{\partial_j u_i \sigma_{ij}} = D_\ell + D_\ell^{{small}}; \\ J_j^{{small}}(\boldsymbol{x}) & =\frac{\bar{\rho}\tilde{\tau}(u_i,u_i)}{2} \tilde{u}_j+\frac{1}{2} \bar{\rho}\tilde{\tau}(u_i,u_i,u_j)+\bar{\tau}(P,u_j)-(\overline{u_i\sigma_{ij}}-\tilde{u}_i\bar{\sigma}_{ij}); \end{aligned} \right\} \end{align}

in which $\tau (\cdot,\cdot )$ is the generalized second-order moment (Germano Reference Germano1992): $\tilde {\tau }(u_i,u_j)=\widetilde {u_i u_j}-\tilde {u}_i\tilde {u}_j$ and $\bar {\tau }(\rho,u_j)=\overline {\rho u_j}-\bar {\rho }\ \bar {u}_j$. Note that our small-scale KE, defined as $\bar {\rho }(\widetilde {|\boldsymbol {u}|^{2}} - |\tilde {\boldsymbol {u}}|^{2})/2$, is sometimes referred to as ‘subscale’ or ‘subgrid’ KE in the LES literature.

Coarse-grained PE density equation (3.4) can be derived by noting that

(3.6)\begin{align} \frac{\textrm{D} \mathcal{P}}{\textrm{D}t} &= \frac{\partial \mathcal{P}}{\partial t} + \boldsymbol{u} \boldsymbol{\cdot} \boldsymbol{\nabla}\mathcal{P}\nonumber\\ &= g z \frac{\partial\rho}{\partial t} + \rho g u_z + g z \boldsymbol{u}\boldsymbol{\cdot} \boldsymbol{\nabla}\rho\nonumber\\ &={-}\mathcal{P} \boldsymbol{\nabla}\boldsymbol{\cdot} \boldsymbol{u} + \rho g u_z. \end{align}

We analyse the energy pathways by inspecting the budget equations (3.1)–(3.4) using our simulations in table 1. In particular, we shall demonstrate the existence of an inertial range in such flows. Note that we do not formulate budgets for small-scale IE and PE, since it is unclear if these two quantities cascade (Eyink & Drivas Reference Eyink and Drivas2018).

3.2. Energy pathways in RT flows

This section presents some of the main results of our paper, which we summarize in § 3.2.5 for the reader's convenience. Here, we refine our analysis of the global energy balances in figure 4 at different scales via coarse-grained budgets. The types of energy involved in the scale transfer, as shown in (3.1)–(3.4), are PE, IE, large-scale KE and also small-scale KE. Differences in the energy pathways between two dimensions and three dimensions will be emphasized during our presentation.

3.2.1. Conversion between different forms of energy

Figure 6 summarizes the RT energy pathways we discuss in this section. In RT flows, the ultimate energy source is the release of PE from the initial unstable density stratification. From (3.4), large-scale PE is a source for large-scale KE via $\epsilon ^{{inj}}$ in (3.1). Small-scale KE injection $\epsilon ^{{small}}_\ell = \bar {\rho }\tilde {\tau }(u_i, g_i)$ is identically zero for RT flows since the gravitational acceleration $\boldsymbol {g}$ is constant in space. Thus, PE is available directly to large-scale KE. Meanwhile, KE is also linked to IE via pressure dilatation and viscous dissipation. The term $\bar {P}_\ell \boldsymbol {\nabla }\boldsymbol {\cdot } \bar {\boldsymbol {u}}_\ell$ appears as a source in the large-scale KE equation (3.1) and a corresponding term $\bar {\tau }(P, \boldsymbol {\nabla }\boldsymbol {\cdot } \boldsymbol {u})$ is a source in the small-scale KE budget equation (3.2), while both appear as sinks in the IE equation (3.3). Similarly, dissipation terms $D_\ell$ and $D_\ell ^{{small}}$ are sinks in large- and small-scale KE budgets, respectively, and sources in the IE budget: $D_\ell ^{{int}} = D_\ell + D_\ell ^{{small}}$. These channels, summarized in figure 6, are responsible for converting energy between different forms. In addition to these conversion channels, the transfer between large- and small-scale KE is now presented.

Figure 6. Schematic of energy pathways in RT flows between reservoirs of KE, PE and IE. The flow is further decomposed into its constituent scales within the KE reservoir, where a vertical dashed line denotes an arbitrary inertial scale $\ell$ we are probing. Kinematically possible pathways are depicted with dashed arrows, while those that are dynamically manifested in 3-D RT are depicted with solid arrows (see §§ 3.2.3–3.2.4). Gravitational PE is converted into KE at the largest scale by $\epsilon ^{{inj}}$. This energy is then transferred to inertial scales by fluxes $\varLambda$ and $\varPi$. On average, $\varLambda$ acts as a conduit for PE, transferring KE to smaller inertial scales in both two dimensions and three dimensions. On average, $\varPi$ transfers energy downscale in three dimensions and upscale in two dimensions. The KE is coupled to IE at largest scales via pressure dilatation, and at the smallest viscous scales via dissipation.

3.2.3. The inertial range in turbulent RT flows

Having identified the pathways, we now present a quantitative analysis. We calculate all terms in the large-scale KE budget (3.1) as a function of scale $\ell$, which allows us to test and validate ideas on energy conversion and the cascade.

The results from the 2D4096 simulation (see table 1) are shown in figure 7, while the 3D1024 results are in figure 8. These figures plot the domain-averaged terms in (3.1) as a function of filtering wavenumber (inverse of scale) $k_\ell =L_z/\ell$. Note that our ‘filtering wavenumber’ $k_\ell$ does not involve any Fourier transform, and we only use it to refer to the corresponding scales, as is conventional in the turbulence literature. For large $k_\ell$ (small scales), these terms approach their unfiltered values. Therefore, they represent cumulative processes acting at all scales larger than $\ell$ and not just the contribution from scale $\ell$ itself. For example, if $\langle \partial _t \overline {\text {KE}}\rangle \equiv \langle \partial _t (\bar {\rho }_\ell ({|\tilde {\boldsymbol {u}}_\ell |^{2}}/{2}))\rangle$ increases in value over the range $[k_1,k_2]$ ($k_1< k_2$), then the KE content of scales within the band $[L_z/k_2,L_z/k_1]$ is growing. The only exceptions are the flux terms, $\varLambda _\ell$ and $\varPi _\ell$, which are hybrid quantities involving coarse-grained fields ($\boldsymbol {\nabla }\bar {P}_\ell$ or $\tilde {\boldsymbol {S}}_\ell$) acting against subscale contributions ($\bar {\tau }_\ell (\rho,\boldsymbol {u}$) or $\tilde {\tau }_\ell (\boldsymbol {u},\boldsymbol {u})$) and, therefore, represent the energy transfer across the partitioning scale $\ell$ (Aluie Reference Aluie2011,  Reference Aluie2013).

Figure 7. Kinetic energy processes as a function of scale in the 2D4096 simulation at dimensionless time $\hat {t}=4.0$. Filtering wavenumber $k_\ell =L_z/\ell$ is a proxy for length scale. Plots are normalized by $\langle \epsilon ^{{inj}} + P\boldsymbol {\nabla }\boldsymbol {\cdot } \boldsymbol {u}\rangle$, the (unfiltered) available mean source of KE. All mean budget terms in the filtered large scale KE equation (3.1) are plotted as a function of $k$. The time derivative of filtered KE $\partial _t \overline {\text {KE}}\equiv \partial _t (\bar {\rho }_\ell ({|\tilde {\boldsymbol {u}}_\ell |^{2}}/{2}))$ is also shown for this unsteady RT flow.

Figure 8. Same as in figure 7 but applied to the 3D1024 data at the same dimensionless time $\hat {t}=4.0$. Plots are normalized by $\langle \epsilon ^{{inj}} + P\boldsymbol {\nabla }\boldsymbol {\cdot } \boldsymbol {u}\rangle$ at that time. Main differences from two dimensions are that here $\langle \varPi \rangle > 0$ (blue), transferring energy to smaller scales and resulting in much higher dissipation (magenta). Also, small scales at $k_\ell >100$ are still growing in time (black) unlike in figure 7 where $\langle \partial _t \overline {\text {KE}}\rangle$ (a cumulative quantity with contributions from all scales $< k_\ell$) becomes $k_\ell$-independent at large $k_\ell$, indicating that there is negligible growth at $k_\ell >100$ with most of the temporal KE increase happening at very larger scales in two dimensions (see figure 7, jump in black curve over $k\in [4,30]$).

Figures 7 and 8 show results from the 2D4096 and 3D1024 simulations, respectively. Spatial mean injection $\langle \epsilon ^{{inj}}\rangle$ shown in the figures is constant in $\ell$, which is due to the spatially uniform gravitational acceleration used here (Aluie Reference Aluie2013). The reason becomes clear when considering that $\langle \epsilon ^{{inj}}\rangle = \langle \bar {\rho } \tilde {u}_i \tilde {g}_i\rangle = -\langle \overline {\rho u_z} \rangle g=-\langle \rho u_z \rangle g$ is independent of $\ell$. This implies that mean injection from PE is released directly only at the largest scale (Aluie Reference Aluie2013). To elaborate, consider two scales, one with length $L$, the largest in the system, and another with a smaller length scale $\ell$. Since mean injection at the two distinct scales is the same, injection within the band of scales $(\ell, L)$ is zero. Taking the limit $\ell \rightarrow 0$, it is obvious that no mean injection occurs at any scale below $L$. Thus, mean injection occurs only at the largest length scale, that of the domain size. Subsequent transfer to smaller scales occurs via the fluxes $\varLambda$ and $\varPi$, as we shall discuss below.

Similar to injection, mean pressure dilatation $\langle -\bar {P}\boldsymbol {\nabla }\boldsymbol {\cdot }\bar {\boldsymbol {u}}\rangle$ is almost flat as a function of $\ell$ in figures 7 and 8. This is consistent with predictions in Aluie (Reference Aluie2011) and simulations by Aluie et al. (Reference Aluie, Li and Li2012), where the pressure dilatation was shown to operate over a transitional scale range of limited extent. Hence, mean pressure dilatation only affects the large-scale dynamics. Dissipation $\langle D_\ell \rangle$ within the Favre scale decomposition has already been proved to be negligible at large scales (Aluie Reference Aluie2013; Zhao & Aluie Reference Zhao and Aluie2018), which is again confirmed here in figures 7 and 8. Thus, in well-developed RT turbulence, there is a range of length scales which is immune from injection and pressure dilatation at large scales, and also from dissipation at small scales. This range of length scales is the ‘inertial range’ in the spirit of Kolmogorov's turbulence theory. Within this inertial range, KE dynamics is governed by the two fluxes $\varPi$ and $\varLambda$.

3.3. Spatial distribution of fluxes

The spatial distribution of values for fluxes $\varLambda$ and $\varPi$ is shown by their p.d.f.s in figure 9, for 2-D RT and 3-D RT. Results at three scales $\ell =L_z/8, L_z/32, L_z/128$ are shown. All p.d.f.s are non-Gaussian, with heavy tails and a stronger departure from Gaussianity at smaller scales. In two dimensions, $\varPi$ is skewed towards negative values, but in three dimensions it is positively skewed, reflecting the inverse versus forward net energy cascade by $\varPi$ in two dimensions versus three dimensions. The p.d.f. of $\varLambda$ is always skewed towards positive values. In particular, at large scales (e.g. $\ell =L_z/8$), $\varLambda _\ell$ is positive everywhere in the domain as is seen in figures 9(b) and 9(d). As we shall discuss in a subsequent paper, negative values of $\varLambda _\ell$ are correlated with spiralling regions of the flow, which are part of the mushroom structure that appears in late time single-mode RT. Note that the p.d.f. of $\varLambda$ in two dimensions (figure 9a) contains both positive and negative values when $\ell$ is below or equal to $L_z/32$, while the p.d.f. of $\varLambda$ in three dimensions (figure 9d) is almost always positive even at $\ell =L_z/32$. This observation is in accord with visualizations of density fields in figures 2(a) and 2(b) where 3-D RT exhibits much finer structures than in 2-D RT, such that the spiralling structures are on average smaller in three dimensions.

Figure 9. Probability density functions (p.d.f.s) of $\varPi$ and $\varLambda$ from simulations 2D4096 and 3D1024 at dimensionless time $\hat {t}=4$. Results for three different scales $\ell =L_z/8, L_z/32, L_z/128$ are shown. The panels are for: (a) the p.d.f. of $\varPi$ in two dimensions; (b) the p.d.f. of $\varLambda$ in two dimensions; (c) the p.d.f. of $\varPi$ in three dimensions; (d) the p.d.f. of $\varLambda$ in three dimensions.

3.4. Temporal self-similarity of flux terms

So far, we have discussed the fluxes as a function of scale at one time instant. In this section we investigate their temporal evolution. Figure 10 shows mean $\varPi$ and $\varLambda$ normalized by $\langle \epsilon ^{{inj}} + P\boldsymbol {\nabla }\boldsymbol {\cdot } \boldsymbol {u}\rangle$ at dimensionless times $\hat {t}=2,3,4$ using the 2D4096 and 3D1024 RT data (see table 1). As RT turbulence develops, the magnitude of normalized $\varLambda$ decreases, whereas the normalized $\varPi$ increases in magnitude. Note, however, that the absolute values of both fluxes keep increasing (not shown). A decrease in normalized $\varLambda$ indicates that as KE at scale $\ell$ saturates in time, the relative importance of $\varLambda$, which is responsible for redistributing injected energy from the potential release, also decreases. In contrast, turbulence levels in RT flows increase in time, leading to an increasing relative magnitude of $\langle \varPi \rangle$. The scale corresponding to the peak magnitude of $\varPi$, indicated by $k_{{peak}}$, shifts to smaller $k_\ell$ (or larger scales) over time. For example, in two dimensions, $k_{{peak}}$ shifts from $\approx 64$ at $\hat {t}=2$, to $\approx 10$ at $\hat {t}=4$, while in three dimensions, $k_{{peak}}$ shifts from approximately $\approx 256$ at $\hat {t}=2$ to $\approx 100$ at $\hat {t}=4$. Thus, as expected, turbulence encroaches onto larger scales as RT evolves, so that vortex stretching (in three dimensions) or mergers (in two dimensions) which drive the $\varPi$ flux are more pronounced at these scales.

Figure 10. Temporal evolution of fluxes in two dimensions and three dimensions at dimensionless time $\hat {t}=2,3,4$. The filtering wavenumber is $k_\ell =L_z/\ell$. Plots are normalized by $\langle \epsilon ^{{inj}} + P\boldsymbol {\nabla }\boldsymbol {\cdot } \boldsymbol {u}\rangle$ at the corresponding time: (a) normalized $\langle \varLambda \rangle$ evolution in two dimensions; (b) normalized $\langle \varPi \rangle$ evolution in two dimensions; (c) normalized $\langle \varLambda \rangle$ evolution in three dimensions; (d) normalized $\langle \varPi \rangle$ evolution in three dimensions.

The above observations indicate possible self-similar (in time) behaviour of turbulent RT fluxes. If we plot the scales corresponding to the peak magnitude of $\langle \varPi \rangle$, which we denote as $\ell _{\varPi, {peak}}$, against dimensionless time $\hat {t}$, we obtain a quadratic scaling in time as shown in figure 11(a,d), similar to the mixing width $h(t)$ in figure 3. Figure 11(a,d) indicate that, in both two dimensions and three dimensions, $\sqrt {\ell _{\varPi, {peak}}} \propto \hat {t} \propto t$, with the coefficients

(3.7a,b)\begin{equation} \ell_{\varPi, {peak}}^{2D} = \alpha_{\rm \pi}^{2D} \mathcal{A} g t^{2} \approx 0.37 h^{2D}(t), \quad \ell_{\varPi, {peak}}^{3D} = \alpha_{\rm \pi}^{3D} \mathcal{A} g t^{2} \approx 0.067 h^{3D}(t), \end{equation}

where $h^{2D}(t), h^{3D}(t)$ are the mixing widths from the 2D4096 and 3D1024 RT data. Note from figure 10 (also figure 11b,e) how $\langle \varPi \rangle$ in two dimensions is skewed to larger scales, unlike in three dimensions, peaking at approximately half the size of the mixing width, $\ell _{\varPi, {peak}}^{2D}\approx h^{2D}/2$ (3.7a,b). This is consistent with the net upscale transfer by $\varPi$ driving the mixing layer growth as discussed in § 3.2.4. In comparison, $\langle \varPi \rangle$ in three dimensions is symmetric, peaking at a scale $O(10)$ smaller than the mixing width, $\ell _{\varPi, {peak}}^{3D}\approx h^{3D}/10$.

Figure 11. Temporal self-similarity of turbulent RT fluxes. Panels (a–c) show results for 2-D RT, (d–f) show 3-D RT. Panels (a,d) plot length scale $\ell _{\varPi, {peak}}$ associated with the peak of $\varPi$ in two dimensions and three dimensions, versus dimensionless time $\hat {t}$. Panels (b,e) show rescaled $\langle \hat {\varPi }\rangle$ in (3.8). Panels (c,f) show rescaled $\langle \hat {\varLambda }\rangle$.

The temporal scaling of $\ell _{\varPi, {peak}}$ in figure 11(a,d) is evidence that the evolution of mean fluxes $\langle \varPi \rangle, \langle \varLambda \rangle$ is self-similar. Recall that self-similarity implies that functions of length and time, such as $\langle \varPi \rangle$ and $\langle \varLambda \rangle$, should collapse to the same curve when properly rescaled (Ristorcelli & Clark Reference Ristorcelli and Clark2004). Specifically, the rescaled fluxes should satisfy

(3.8)\begin{equation} \left. \begin{aligned} \langle\hat{\varPi}\rangle(\hat{k}_\ell)\equiv f_{\rm \pi}(k_\ell\,t^{2}) & = F_{\rm \pi}^{{-}1}(t)\langle\varPi\rangle\left(\hat{k}_\ell\frac{L_z}{h},t\right),\\ \langle\hat{\varLambda}\rangle(\hat{k}_\ell)\equiv f_\lambda(k_\ell\,t^{2}) & =F_\lambda^{{-}1}(t)\,\langle\varLambda\rangle\left(\hat{k}_\ell\frac{L_z}{h},t\right), \end{aligned} \right\} \end{equation}

where $\hat {k}_\ell$ is just filtering wavenumber $k_\ell =L_z/\ell$, rescaled by the mixing width $h(t)$, $\hat {k}_\ell \equiv k_\ell \,h/L_z$. Here $F_{\rm \pi} (t)$, $F_\lambda (t)$ are temporal scaling functions, and $f_{\rm \pi} (\cdot ), f_\lambda (\cdot )$ are the similarity functions. If mean fluxes $\langle \varPi \rangle, \langle \varLambda \rangle$ are temporally self-similar, then $\langle \hat {\varPi }\rangle, \langle \hat {\varLambda }\rangle$ should collapse onto the same curve.

Figure 11(b,e) for $\langle \hat {\varPi }\rangle$ in two dimensions and three dimensions, and figure 11(c,f) for $\langle \hat {\varLambda }\rangle$ in two dimensions and three dimensions indicate that they indeed collapse onto each other. In figure 11(b,e), the horizontal axis is rescaled by the mixing width $h(t)/L_z$. To account for the temporal function $F_{\rm \pi} (t)$, we rescaled the plots at $\hat {t}=3, 4$ so that their maxima match with that at $\hat {t}=2$. The same rescaling amplitudes for $\langle \hat {\varPi }\rangle$ are also applied to $\langle \hat {\varLambda }\rangle$, as is shown in figure 11(c,f), which indicates that $F_{\rm \pi} (t)\approx F_\lambda (t)$. We have found that $F_{\rm \pi} (t)$ and $F_\lambda (t)$ follow a scaling close to $t^{3}$ both in two dimensions and three dimensions (not shown). This $t^{3}$ growth of the fluxes can be understood from the cumulative amount of PE released, which follows $\delta P\sim h^{2}\sim t^{4}$ (Dimonte et al. Reference Dimonte2004), implying that the rate of KE injection should scale as $\langle \epsilon ^{{inj}}\rangle \sim t^{3}$. This energy is eventually channelled by the fluxes $\varLambda$ and $\varPi$, which grow in magnitude accordingly.

While plots of the fluxes show good collapse in figure 11, it is worth discussing the relatively poor collapse of $\langle \hat {\varPi }\rangle$ in three dimensions (figure 11e) at scales smaller than that of the peak. The reason $\langle \hat {\varPi }\rangle$ decays at those scales is due to viscous dissipation removing the KE available to cascade. As the RT flow evolves and there is more KE cascading, the dissipation scale $\ell _d$ should become smaller as $\ell _d\sim t^{-1/4}$ as shown in Ristorcelli & Clark (Reference Ristorcelli and Clark2004) by invoking Kolmogorov's turbulence phenomenology. This is reflected in $\langle \hat {\varPi }\rangle$ extending to smaller scales at later times in figure 11(e). Therefore, the dissipation dynamics is not expected to follow the temporal self-similarity of inertial scales captured in (3.8). Note that the $\langle \hat {\varPi }\rangle$ cascade in two dimensions, being upscale, is not arrested by viscous dissipation at the smallest scales, which explains the much better collapse of $\langle \hat {\varPi }\rangle$ in figure 11(b).

4.1. Directionally split anisotropy at scales

Here, we investigate anisotropy of RT flows at different scales in physical space. To quantify the contribution to the KE budget from different directions, we first formulate the budgets along each direction separately, and then diagnose the constituent terms from numerical data. We call this ‘directionally split’ analysis of anisotropy, which is different from the analysis of scale anisotropy done in supplementary material, section B. Similar to the filtered total KE equation (3.1), we can get the filtered KE equations in separate directions,

(4.1)\begin{equation} \partial_t\bar{\rho}_\ell \frac{|\tilde{\boldsymbol{u}}_\xi |^{2}}{2}+\boldsymbol{\nabla}\boldsymbol{\cdot} \boldsymbol{J}_\ell^{\xi} ={-}\varPi^{\xi}_\ell-\varLambda^{\xi}_\ell+\bar{P}_\ell \partial_\xi\bar{u}_\xi-D^{\xi}_\ell+\epsilon_\ell^{inj,\xi}, \end{equation}

where $\xi$ is a symbol denoting a single direction $x$, $y$ or $z$. Note that since we are not concerned with scale anisotropy here, the filtering kernel is still isotropic as before, $G(\boldsymbol {r}) = G(|\boldsymbol {r}|)$. The budget terms in the above equation are defined by

(4.2)\begin{equation} \left. \begin{aligned} & \varPi^{\xi}_\ell(\boldsymbol{x})={-}\bar{\rho} \partial_j \tilde{u}_\xi \tilde{\tau}(u_\xi,u_j), \quad \varLambda^{\xi}_\ell(\boldsymbol{x})=\frac{1}{\bar{\rho}}\partial_\xi \bar{P} \bar{\tau}(\rho,u_\xi), \quad D^{\xi}_\ell(\boldsymbol{x})=\partial_j \tilde{u}_\xi \bar{\sigma}_{\xi j}, \\ & J^{\xi}_j(\boldsymbol{x})=\bar{\rho} \frac{\tilde{u}_\xi^{2}}{2}\tilde{u}_j+\bar{P}\bar{u}_\xi\delta_{\xi j}+\tilde{u}_\xi \bar{\rho} \tilde{\tau}(u_\xi,u_j)-\tilde{u}_\xi\bar{\sigma}_{\xi j}, \quad \epsilon^{inj,\xi}_\ell(\boldsymbol{x})=\tilde{u}_\xi\bar{\rho}\tilde{F}_\xi, \end{aligned} \right\} \end{equation}

whose forms are similar to their counterparts in the total KE budget equation (3.1). The sum of these equations in all directions reduces to the full budget equation (3.1). Note that in the above definitions, the subscript $\xi$ represents only one direction, and $\xi$ itself does not follow the Einstein summation convention.

Before delving into details of the budgets in separate directions, we shall first check the coarse-grained KE in each direction. This is shown as a function of filtering wavenumber $k_\ell =L_z/\ell$ in figure 12 for 2-D RT and 3-D RT data. In two dimensions, mean KE is approximately isotropic at all scales, with the horizontal component slightly greater than the vertical at large scales. This may seem counterintuitive since KE is injected in the vertical direction. In contrast, we see that in 3-D, horizontal KE in the $x$ and $y$ directions is almost equal, but is more than three times smaller than KE in the vertical at large scales. This indicates that in the horizontal plane, 3-D RT turbulent flow is isotropic as expected, but the vertical flow is much stronger than the horizontal at large scales, consistent with results of Cabot & Zhou (Reference Cabot and Zhou2013) based on Fourier spectra in the homogeneous direction. Thus, we have different behaviours for 2-D RT and 3-D RT. In 2-D RT, the flow seems to be close to isotropic, while it is highly anisotropic in 3-D RT. Note that in figure 12, derivatives of the plots with respect to $k_\ell$ give rise to the filtering spectra defined in (2.14). At small scales, the $k_\ell$-derivatives of both 2-D and 3-D filtered KE are small and comparable in all directions. At large scales, the filtering spectrum indicates that the flow is approximately isotropic in two dimensions (figure 12a) and highly anisotropic in three dimensions (figure 12b), where the vertical large-scale flow contains more energy than the horizontal large-scale flow. We will now investigate this difference carefully by analysing the energy budgets in separate directions.

Figure 12. Directional anisotropy of RT flows. Coarse-grained KE of the flow in different directions as a function of scale in two dimensions and three dimensions. Their sum over all directions is also shown for reference. Data is at dimensionless time $\hat {t}=4.0$. These are cumulative quantities, showing energy at all scales larger than $\ell =L_z/k_\ell$. Energy content at any $k_\ell$ (i.e. the spectrum) is obtained from the $k_\ell$-derivative of these plots. Panel (a) is for the 2-D result and (b) for the 3-D result.

In figure 13 we show mean values of the normalized budget terms in (4.2) for the 2D4096 simulation at time $\hat {t}=4.0$. Budget terms corresponding to the evolution of horizontal KE, vertical KE and the full KE are shown. In figure 13(a), $\varPi _x$ is larger in magnitude than $\varPi _z$, in accordance with the larger magnitude of KE in the horizontal direction in 2-D RT. The $\varLambda$ term resides mostly in the vertical component $\varLambda _z$, while $\varLambda _x$ is almost zero. This indicates that $\varLambda$ is mainly responsible for transferring injected energy to vertical scales, and does not redistribute KE among different directions.

Figure 13. The 2-D RT directionally split budget terms for filtered KE using the 2D4096 data at $\hat {t}=4$. Filtering wavenumber $k_\ell =L_z/\ell$. Plots are normalized by $\langle \epsilon ^{{inj}} + P\boldsymbol {\nabla }\boldsymbol {\cdot } \boldsymbol {u}\rangle$. Panels show deformation work, baropycnal work, pressure dilatation and dissipation.

The next term in figure 13 is the (negative) pressure dilatation, which is a cumulative quantity, similar to filtered KE in figure 12, meaning that its contribution at scale $\ell$ is determined by its derivative with respect to filtering wavenumber $k_\ell =L_z/\ell$. That is, its contribution at scale $\ell$ equals $({\partial }/{\partial k})\langle \bar {P}_\ell \boldsymbol {\nabla }\boldsymbol {\cdot }\bar {\boldsymbol {u}}_\ell \rangle \rvert _{k=L_z/\ell }$. Therefore, as indicated by figure 13(c), pressure dilatation acts as a source for mean horizontal KE, since $\langle \bar {P} \partial _x \bar {u}_x\rangle$ increases with respect to scale, so that its contribution to horizontal KE at scale $\ell =L_z/k_\ell$ measured by the derivative is positive. For the mean vertical KE, the vertical component $\langle \bar {P} \partial _z \bar {u}_z\rangle$ at the largest scale is positive and acts as a source, but at all subsequent scales its mean value decreases with increasing filtering wavenumber, $({\textrm {d}}/{\textrm {d}k})\langle \bar {P} \partial _z \bar {u}_z\rangle < 0$, indicating it is a sink term at these scales. The combined value in all directions gives rise to a pressure dilatation that is almost constant on average with respect to filtering wavenumber. Therefore, pressure dilatation acts to mix up horizontal and vertical velocities and isotropize the KE in different directions, similar to the role of pressure strain in constant-density turbulence (Pope Reference Pope2001). Considering dissipation in figure 13, the two components have comparable magnitude, but since there is little KE residing at small scales due to the fast decay of spectrum and a net upscale cascade by $\varPi$, the magnitude of dissipation in two dimensions is far smaller than other budget terms.

The directionally split budget terms for 3-D RT are shown in figure 14. In all four panels, the $x$ and $y$ components are equal, indicating isotropy in the horizontal plane. However, anisotropy in the vertical direction is more pronounced than in two dimensions. Furthermore, $\langle \varPi _z\rangle$ in the vertical direction is much larger than the horizontal components, contrary to the 2-D case where the horizontal $\langle \varPi _x\rangle$ is larger in magnitude. Also in three dimensions, $\langle \varPi _x\rangle$ and $\langle \varPi _y \rangle$ peak at smaller scales compared with $\langle \varPi _z\rangle$, indicating a strengthening of the cascade in the horizontal direction at smaller scales, where it becomes more comparable to $\langle \varPi _z \rangle$. Mean $\varLambda$ is only significant in the vertical direction, with almost zero values in the horizontal directions, similar to two dimensions. Also, pressure dilatation is the source of KE input in the horizontal directions, similar to two dimensions. Dissipation is more isotropic in all directions, consistent with a trend towards isotropy of $\langle \varPi \rangle$ at these small scales, although the vertical component attains a slightly higher value.

Figure 14. Same as in figure 13 but for 3-D RT using the 3D1024 data at $\hat {t}=4$.

In summary, the directionally split energy budgets in 2-D RT and 3-D RT indicates that mean injection and its conduit $\varLambda$ contribute only to the vertical component of KE. Horizontal KE is forced by the horizontal component of pressure dilatation in both two dimensions and three dimensions. Dissipation, which is at the smallest scales, is approximately isotropic. However, significant differences between two dimensions and three dimensions appear in the fluxes $\langle \varPi _x\rangle$ and $\langle \varPi _z\rangle$ at inertial scales.

In the earlier discussion of this section, we observed a sharp contrast between two dimensions and three dimensions in figure 12, which indicated that the 2-D flow is surprisingly isotropic at large scales, while it is highly anisotropic in three dimensions due to the expected stronger vertical velocity. The reason can be gleaned from the foregoing results, especially those in figures 13 and 14. The primary difference between the directionally split budgets in two dimensions and three dimensions is due to the cascade $\varPi$. In two dimensions, figure 13(a) shows that $\varPi _x$ is an energy source for mean horizontal large-scale KE, and similarly, $\varPi _z$ is a source for vertical large-scale KE. However, the magnitude of $\langle \varPi _x\rangle$ is approximately two to three times greater than $\langle \varPi _z\rangle$, indicating that large-scale KE gains more energy from $\varPi$ in the horizontal direction than in the vertical, offsetting anisotropy due to $\varLambda$, which energizes the vertical flow. Thus, the inverse cascade $\varPi$ acts alongside pressure dilatation to isotropize the large-scale flow, where most of the KE resides. In contrast, results from 3-D RT indicate that $\langle \varPi _x\rangle$ and $\langle \varPi _y\rangle$ are positive, acting as an energy sink for horizontal large-scale KE. Therefore, the large-scale horizontal flow in three dimensions has only pressure dilatation as a source of energy, unlike the vertical flow that is being driven by $\varLambda$, which ultimately leads to the large anisotropy in 3-D RT.

While the anisotropy analysed above is due to horizontal and vertical components of the flow, another measure of anisotropy is that of length scales by using anisotropic filtering kernels as we do in the supplementary material, section B.