2.1. Evolution equations

In the following, we consider two fluids initially separated by a horizontal perturbed interface, with the heavier fluid on top of the lighter one in a downward gravity field $\boldsymbol g = - g {\boldsymbol n_3}$, with $\boldsymbol n_3$ the vertical unit vector. The incompressible Navier–Stokes equations within the Boussinesq and MHD approximations read

(2.1a)$$\begin{gather} \frac{\partial \boldsymbol U}{\partial t} + (\boldsymbol U \boldsymbol{\cdot} \boldsymbol \nabla) \boldsymbol U ={-} \boldsymbol \nabla P - 2 \mathcal{A} g C {\boldsymbol n_3} + (\boldsymbol \nabla \times \boldsymbol B ) \times \boldsymbol B + \nu \nabla^2 \boldsymbol U, \end{gather}$$

(2.1b)$$\begin{gather}\frac{\partial \boldsymbol B}{\partial t} + (\boldsymbol U \boldsymbol{\cdot} \boldsymbol \nabla) \boldsymbol B = (\boldsymbol B \boldsymbol{\cdot} \boldsymbol \nabla) \boldsymbol U + \eta \nabla^2 \boldsymbol B, \end{gather}$$

(2.1c)$$\begin{gather}\frac{\partial C}{\partial t} + (\boldsymbol U \boldsymbol{\cdot} \boldsymbol \nabla) C = \kappa \nabla^2 C, \end{gather}$$

(2.1d)$$\begin{gather}\boldsymbol \nabla \boldsymbol{\cdot} \boldsymbol U = 0, \end{gather}$$

(2.1e)$$\begin{gather}\boldsymbol \nabla \boldsymbol{\cdot} \boldsymbol B = 0, \end{gather}$$

where $\boldsymbol U$ is the total velocity field, $\boldsymbol B$ the magnetic field scaled as a velocity ($\boldsymbol B = {\boldsymbol H}/\sqrt {\rho _0 \mu _0}$, where ${\boldsymbol H}$ is the true magnetic field, $\mu _0$ the vacuum magnetic permeability and $\rho _0=(\rho _1 + \rho _2)/2$ the reference density) and $C$ the concentration of the heavy fluid. Additionally, $P$ is the reduced pressure, and $\nu$, $\eta$ and $\kappa$ are the kinematic viscosity, magnetic diffusivity and molecular diffusivity, respectively.

For the Rayleigh–Taylor configuration, we assume homogeneity of the statistics in horizontal planes, so that we can decompose any field $A$ into its horizontal average $\langle A \rangle$ and its fluctuating part $a=A- \langle A \rangle$. The volume average is denoted by $\bar {\star }$ and is obtained by an additional integration over the vertical direction, weighted by the mixing zone length $L$ defined in (2.4). For example, the volume average of the correlation between two fluctuating fields $a_1$ and $a_2$ read

(2.2)\begin{equation} \overline{a_1 a_2}(t) = \frac{1}{L} \int \langle a_1 a_2 \rangle (x_3, t) \, \mathrm{d}\kern 0.06em x_3.\end{equation}

Note that fluctuations decrease quickly outside the mixing layer, so that integration along the whole vertical direction amounts to averaging upon the mixing zone width for second-order correlations.

In the particular case of the magnetic Rayleigh–Taylor instability, we consider the additional presence of a vertical uniform mean magnetic field, so that the total magnetic field reads $\boldsymbol B = B_0 {\boldsymbol n_3} + \boldsymbol b$ (with $\boldsymbol B$ normalized as a velocity, $B_0$ can be interpreted as the mean Alfvén velocity). It is also convenient to express the Lorentz force $(\boldsymbol \nabla \times \boldsymbol B ) \times \boldsymbol B = - \boldsymbol \nabla (B^2/2) + (\boldsymbol B \boldsymbol {\cdot } \boldsymbol \nabla ) \boldsymbol B$ as the sum of magnetic pressure, which is included in the total pressure, and magnetic tension. In addition, note that in this configuration, there are no induced mean magnetic fields, so that $B_0$ remains constant, see Appendix A.1. After some algebra, the equations for the fluctuations read, in the absence of mean velocity and magnetic fields,

(2.3a)$$\begin{gather} \frac{\partial u_i}{\partial t} + \frac{\partial u_j u_i}{\partial x_j} ={-} \frac{\partial \varPi}{\partial x_i} - 2 \mathcal{A} g c \delta_{i3} + \frac{\partial b_i b_j}{\partial x_j} + B_0 \frac{\partial b_i}{\partial x_3} + \frac{\partial}{\partial x_3} \langle u_i u_3 - b_i b_3 \rangle + \nu \nabla^2 u_i, \end{gather}$$

(2.3b)$$\begin{gather}\frac{\partial b_i}{\partial t} + \frac{\partial u_j b_i}{\partial x_j} = \frac{\partial b_j u_i}{\partial x_j} + B_0 \frac{\partial u_i}{\partial x_3} + \frac{\partial}{\partial x_3} \langle u_3 b_i - u_i b_3 \rangle + \eta \nabla^2 b_i, \end{gather}$$

(2.3c)$$\begin{gather}\frac{\partial c}{\partial t} + \frac{\partial u_j c}{\partial x_j} ={-} u_3 \frac{\partial }{\partial x_3} \langle C \rangle + \frac{\partial}{\partial x_3} \langle u_3 c \rangle + \kappa \nabla^2 c. \end{gather}$$

In these equations, summation upon repeated indices is used. Both mean velocity and magnetic fields are negligible in this configuration. Here, $\varPi = p + (b_j b_j - \langle b_j b_j \rangle )/2$ is the generalized pressure, with $p=P-\langle P \rangle$ the fluctuating reduced pressure. Equations for the two-point correlations are given in Appendix A.2. In addition to the usual nonlinear and dissipative terms, the above set of (2.3a)–(2.3c) contains several linear terms: the production of velocity and magnetic fluctuations through the mean magnetic field, namely $B_0 \partial _3 \boldsymbol u$ and $B_0 \partial _3 \boldsymbol b$, and also the production of scalar fluctuations through the enlargement of the turbulent mixing zone via $-u_3 \partial _3 \langle C \rangle$. Note that the inverse mean concentration gradient can be used to evaluate the mixing zone width, with $1/L \simeq \partial _3 \langle C \rangle$. In practice, here we use the integral definition (Andrews & Spalding Reference Andrews and Spalding1990)

(2.4)\begin{equation} L(t) = 6 \int \langle C \rangle ( 1 - \langle C \rangle ) \,\mathrm{d}\kern 0.06em x_3.\end{equation}

Equations (2.3) indicate that for an initially perturbed interface between the heavy and light fluids, vertical velocity fluctuations $u_3$ are created through the buoyancy term $- 2 \mathcal {A} g c$, which then produces additional concentration fluctuations through $-u_3 \partial _3 \langle C \rangle$ and generates magnetic fluctuations as well through the term $-B_0 \partial _3 u_i$. In the following, we will investigate how the induced magnetic field can alter the development of the turbulent mixing zone.

2.2. Numerical set-up

We have performed direct numerical simulations (DNS) of the incompressible Navier–Stokes equations (2.1a)–(2.1e) under the Boussinesq and MHD approximations, with the pseudo-spectral code STRATOSPEC (Briard, Gostiaux & Gréa Reference Briard, Gostiaux and Gréa2020). The domain is a triply periodic box of dimensions $(L_{x_1}, L_{x_2}, L_{x_3})$, with always $L_{x_1} = L_{x_2} = 2 {\rm \pi}$.

The present simulations have resolution from $1024^2 \times 2048$ points to $4096^3$ points, and are performed with a pencil decomposition, using from 2048 up to 65 536 cores. A classical spectral Fourier collocation method is used with two-third rule dealiasing, so that the maximum wavenumber is $k_{max}=N_1/3$, with $N_1$ the number of points in the $x_1$ direction. The P3DFFT algorithm was implemented to perform massively parallel fast Fourier transforms (Pekurovsky Reference Pekurovsky2012). The time increment is determined using a third-order low-storage strong-stability-preserving Runge–Kutta scheme, with implicit treatment of diffusive terms.

The initial conditions are detailed in Appendix B.1. The main features are that the velocity and fluctuating magnetic fields are initially zero. The initial perturbation of the interface between the heavy and light fluids is driven by a three-dimensional concentration variance spectrum $E_{cc}$, whose peak wavenumber $k_p$ and initial variance $\overline {c^2}(t=0)$ are given in table 1.

Table 1.

Parameters of the DNS: Number of points $N$, Vertical extent $L_{x_3}$ (with always $L_{x_1}=L_{x_2}=2{\rm \pi}$), Buoyancy strength $\mathcal {A}g$, with $\mathcal {A}=0.05$ the Atwood number and $g$ the acceleration, Mean vertical magnetic field intensity $B_0$ (scaled as a Alfvén velocity), Kinematic viscosity $\nu = \kappa$, Magnetic Prandtl number $P_M = \nu /\eta$, and Peak wavenumber $k_p$ of the initial concentration spectrum (all runs have initially $\overline {c^2} (t=0) = 0.022$ and $\varTheta (t=0)=0.867$), Self-similar value of the magnetic to kinetic energy ratio $R_{ss}$, Ratio of the maximum and Kolmogorov wavenumbers $k_{max}/k_\eta$ at the final time $t_{end}$ of the simulation.

The parameters of the different DNS are gathered in table 1. The simulations starting with label R1 refer to resolution $1024^2 \times 2048$, with a computational domain vertically elongated, with $L_{x_1}=L_{x_2} = L_{x_3}/2=2{\rm \pi}$. Some of these simulations were duplicated in the cubic domain $L_{x_1}=L_{x_2} = L_{x_3}=2{\rm \pi}$, with resolution $1024^3$, but are not presented in table 1 for the sake of conciseness. They only serve to estimate vertical confinement effects, to address the spatial resolution of the simulations in Appendix B.2 and for illustration purposes, like in figure 2. Hence, only the R1 simulations in the elongated domain are presented in table 1 since they permit to reach a fully developed turbulent asymptotic state.

Figure 2.

Concentration field $C$ from the simulations R1B0G10, R1B020G10 and R1B050G10 ($1024^3$) for different mean magnetic field intensities – from left to right, $B_0=0$, $B_0=0.20$ and $B_0=0.50$ – at four different instants – from top to bottom, $t=0$, $t=2$, $t=5$ and $t=8$. The colourmap is chosen so that the pure fluids ($C=0$ and $C=1$) are transparent, with a sharp colour gradient around $C=1/2$.

Simulations R2 with resolution $2048^3$ and cubic domain $(2{\rm \pi} )^3$ are used to address the spatial resolution in Appendix B.2, along with the effects of changing either the peak wavenumber of the initial condition or the diffusion coefficients in Appendix D. The only simulation with resolution $2048^2 \times 4096$ serves to reach a more turbulent self-similar state of the MRTI. Finally, simulations R4 have the finest resolution with $4096^3$ points and are used to address spectral scalings and scale-by-scale energy budget.

The subsequent letters in the labels are B for the mean magnetic field intensity $B_0$, and G for the value of the gravitational acceleration $g$ (note that $\mathcal {A}$ is not modified in this study). For some cases, an additional letter a or b is used to discriminate between different diffusion coefficients. Note that the peak wavenumber $k_p$ is different for the R1, R2 and R4 simulations.

3.1. Observations

Before performing a thorough analysis of the flow generated by the magnetic Rayleigh–Taylor instability, we first show in figure 2 the 3-D concentration field $C$, extracted from the R1 simulations at different instants, with varying mean magnetic field intensity $B_0$. The colourmap is chosen so that the pure fluids ($C=0$ and $C=1$) are transparent, with a sharp colour gradient around $C=1/2$. From left to right, $B_0$ is increased, and from top to bottom, time goes by. A supplementary movie is available at https://doi.org/10.1017/jfm.2023.1053 of the developing instability.

Starting from the same initial condition (the same perturbed interface), it is clear that at $t=2$, the onset of the Rayleigh–Taylor instability is already greatly modified by the presence of the mean magnetic field. In the purely hydrodynamic case, typical bubbles of different sizes appear. For $B_0=0.2$, the initial perturbations are vertically stretched. For the even stronger Alfvén velocity $B_0=0.5$, the instability is delayed, with the growth rate of the perturbations damped by the mean magnetic field (Chandrasekhar Reference Chandrasekhar1961).

Later on, at $t=5$, turbulent structures are visible for the non-magnetic case, with multiple colours on the 3-D concentration field, illustrating mixing at different levels. For $B_0=0.2$, bubbles have also formed at the edges of the developing elongated structures, which shows some resemblance with the previous case $B_0=0$ at $t=2$, but with a wider mixing layer. For $B_0=0.5$, the initial perturbations have only been stretched, revealing smooth elongated fingers, with apparently no mixing.

At $t=8$, the turbulent mixing zone has grown in size for both $B_0=0$ and $B_0=0.2$, but interestingly, it is wider for the latter case. Despite the initial inhibition of the shear instabilities, turbulent mixing has been enhanced afterwards. Finally, for $B_0=0.5$, the smooth elongated fingers are still being stretched, with bubbles eventually forming at their edges, with no mixing between ascending and descending structures.

The phenomenology depicted here is interesting and reveals several noteworthy features. Compared with the purely hydrodynamic case, the mean vertical magnetic field initially delays the instability, but later enhances turbulent mixing, as seen with $B_0=0$ and $B_0=0.2$ at $t=5$. This was first evidenced in the early 2-D simulations of Jun et al. (Reference Jun, Norman and Stone1995) and is confirmed here in 3-D high-resolution simulations. In the regime where magnetic effects dominate, structures are elongated and smooth fingers appear, rather than turbulent bubbles; small-scale shear instabilities are inhibited and there is no turbulent mixing, yielding a rapid growth of the mixing zone $L$.

In the two following parts, in a simplified and idealized framework, a criterion for the fingers to be destabilized is determined, and some explanations are proposed for the appearance of elongated structures.

3.2. Shear instability between ascending and descending structures

In this part, we perform a simplified linear stability analysis to provide a criterion for which the vertically elongated fingers observed previously are destabilized by small-scale shear instabilities.

We consider a 2-D inviscid framework and choose to model the two elongated fingers by two semi-infinite planes, separated at some initial time $t_0$ by a flat interface. The falling pure heavy fluid of density $\rho _1$ and the ascending light one of density $\rho _2$, within the Boussinesq approximation, have opposite constant velocities $\pm V$. The mean magnetic field $B_0$ is oriented upward and the gravitational acceleration $g$ downward. The configuration is sketched in figure 3 with, for this section only, $x_1=0$ indicating the interface position ‘at rest’ between ascending and descending structures. It is worth noting that the framework just described corresponds in fact, after a $90^\circ$ clockwise rotation, to an almost classical Kelvin–Helmholtz configuration (Chandrasekhar Reference Chandrasekhar1961), with a mean magnetic field parallel to the interface, but in which the gravity would be horizontal, opposite to $B_0$.

Figure 3.

Sketch of the 2-D configuration for the present linear stability analysis.

The details of the linear stability analysis are gathered in Appendix C for the sake of conciseness and only the result is given here: the interface between the ascending and descending fluids is unstable if

(3.1)\begin{equation} V^2 > B_0^2,\end{equation}

meaning that the typical velocity of the fingers should exceed the Alfvén velocity imposed by the vertical mean magnetic field. This simple criterion, also proposed by Jun et al. (Reference Jun, Norman and Stone1995), is compared with DNS in § 4. It provides only a lower bound for transition to turbulence to occur. Indeed, the real elongated structures are not semi-infinite planes, have finite thickness and are separated by a diffuse interface. These features participate in the reduction of the instability growth rate (Miura & Pritchett Reference Miura and Pritchett1982). Hence we expect transition to turbulence to occur for velocities larger than $B_0$.

Finally, note that the stability of a shear flow with an aligned mean magnetic field was addressed by Hughes & Tobias (Reference Hughes and Tobias2001) for arbitrary background velocity and magnetic fields. In the simple case of a uniform $B_0$ and constant piecewise velocity field, it reduces to (3.1).

3.3. Homogeneous turbulence approach

In the previous section, the destabilization of the elongated fingers due to the shear between ascending and descending structures has been investigated. When these smooth fingers break and the flow becomes turbulent, it is expected that due to magnetic tension, structures remain aligned with magnetic field lines, namely stretched in the vertical direction. To analyse this phenomenology in more detail, we use the spectral formalism, which is a corner stone of the developments and predictions made later in § 6.4. For simplicity, we assume that turbulent fluctuations within the core of the mixing layer are homogeneous in the three directions and that they ‘see’ only the mean concentration gradient $\partial _3 \langle C \rangle =1/L$: in the non-magnetic case, this framework is known as unstably stratified homogeneous turbulence (Griffond, Gréa & Soulard Reference Griffond, Gréa and Soulard2014; Soulard, Griffond & Gréa Reference Soulard, Griffond and Gréa2014; Burlot et al. Reference Burlot, Gréa, Godeferd, Cambon and Griffond2015; Griffond, Gréa & Soulard Reference Griffond, Gréa and Soulard2015; Gréa et al. Reference Gréa, Burlot, Godeferd, Griffond, Soulard and Cambon2016; Briard, Iyer & Gomez Reference Briard, Iyer and Gomez2017). Then, we discard the nonlinearities along with diffusion terms, which do not act at the large scales considered, keeping only interactions with the mean flow, which is reminiscent of the rapid distortion theory (Hunt & Carruthers Reference Hunt and Carruthers1990). The resulting equations for the spectral vertical fluctuations are

(3.2a)$$\begin{gather} \frac{\partial \hat{u}_3}{\partial t} ={-} 2 \mathcal{A}g \hat{c}(\boldsymbol k) P_{33}(\boldsymbol k) + \mathrm{i} k_3 B_0 \hat{b}_3(\boldsymbol k), \end{gather}$$

(3.2b)$$\begin{gather}\frac{\partial \hat{b}_3}{\partial t} = \mathrm{i} k_3 B_0 \hat{u}_3(\boldsymbol k), \end{gather}$$

(3.2c)$$\begin{gather}\frac{\partial \hat{c}}{\partial t} ={-} \frac{1}{L} \hat{u}_3(\boldsymbol k), \end{gather}$$

where $\hat {{\cdot }}$ represents the Fourier transform and $P_{ij}=\delta _{ij}-k_i k_j/k^2$ the projector onto the plane perpendicular to the wavevector $\boldsymbol k$. The complete spectral equations for the fluctuating fields are given in Appendix A.3. Combining the previous equations readily gives

(3.3)\begin{equation} \frac{\partial^2 \hat{u}_3}{\partial t^2} + ( k^2 B_0^2 \cos^2 \theta - N^2 \sin^2 \theta ) \hat{u}_3 = 0, \end{equation}

where $\theta$ is the angle between $\boldsymbol k$ and the vertical ${\boldsymbol n}_3$, and $N=\sqrt {2 \mathcal {A}g/L}$ the buoyancy frequency (similar equations for $\hat {c}$ and $\hat {b}_3$ can be obtained). In the absence of a magnetic field, it is recovered that the modes with wavevector $\boldsymbol k$ pertaining to the equatorial plane ($\theta ={\rm \pi} /2$) are the most amplified (Soulard et al. Reference Soulard, Griffond and Gréa2014). The striking difference in the magnetic case is that for all other orientations, $\theta \neq {\rm \pi}/2$, the modes are damped by the mean magnetic field, which forces the structures to remain aligned with the field lines.

4.1. Overall dynamics, anisotropy and mixing

In this part, we focus on one-point statistics, which provide quantitative information about the dynamics, transition to turbulence and mixing. Additional details are given in Appendix D, notably regarding the Reynolds number, which can be defined in multiple ways.

The criterion (3.1), derived from the previous inviscid linear analysis, is now compared with DNS: we choose $V = (\overline {u_3^2})^{1/2}$ for the typical velocity of the ascending and descending fingers. It is revealed in figure 4(a) that for $V>B_0$, there is indeed a change of regime that occurs later for larger $B_0$. Elongated smooth fingers are destabilized by shear instabilities, mixing occurs, causing a net reduction of the vertical velocity: this is transition to turbulence. Afterwards, the growth rate of the vertical velocity is decreased compared to the initial phase. For the case $B_0=0.5$, $V$ hardly reaches the Alfvén velocity, consistent with the observations of figure 2. In contrast, for the weak field $B_0=0.1$, the initial growth is quite comparable to the RTI case so that $V$ reaches large values before destabilization. As expected from the discussion at the end of § 3.2, transition to turbulence happens after the criterion (3.1) is fulfilled. The ratio $M_A=V/B_0$ can be interpreted as an Alfvén Mach number and shows that the turbulent state of the MRTI corresponds to $M_A>1$.

Figure 4.

Effects of the mean magnetic field intensity $B_0$ in the R1 simulations. (a) Normalized vertical velocity with the Alfvén velocity $\sqrt {\overline {u_3^2}}/B_0$, see criterion (3.1). (b) Mixing zone size $L$ (2.4).

We pursue the analysis by investigating the mixing layer growth in figure 4(b). As observed in figure 2, the presence of a vertical mean magnetic field slows down the onset of the instability: a stronger $B_0$ results in a slower initial growth. This is qualitatively in agreement with the theory since the smallest wavelengths have their growth rate damped (Chandrasekhar Reference Chandrasekhar1961). Later on, with the magnetic field acting as a tension that inhibits small-scale shear instabilities, the elongated fingers grow more rapidly than the merging bubbles. With larger $B_0$, transition to turbulence is delayed and the width of the mixing layer can be significantly larger than in the hydrodynamic case. Then, turbulent mixing occurs and the growth of $L(t)$ is slowed down. Dimensional analysis indicates that two regimes are possible for the mixing layer: either $L \sim B_0 t$ or $L \sim \mathcal {A}g t^2$. Consistent with the previous discussion, the former corresponds to the transient regime of rapid growth with $M_A<1$: structures are rapidly stretched without mixing and grow at a constant velocity (Gréa & Briard Reference Gréa and Briard2023), whereas the latter, which is the standard self-similar regime of the RTI (Youngs Reference Youngs1994), is recovered after transition to turbulence and corresponds to $M_A>1$. The turbulent mixing layer growth rate, modified by the presence of induced magnetic energy, is thoroughly analysed later in § 6.

It has been observed in figure 2 that a larger $B_0$ means that mixing occurs later, since the mean magnetic field inhibits shear instabilities at small scales. To quantify mixing, we use the mixing parameter $\varTheta$ (Youngs Reference Youngs1994) defined as

(4.1)\begin{equation} \varTheta = \frac{\displaystyle \int \langle C (1-C) \rangle \,\mathrm{d}\kern 0.06em x_3}{\displaystyle \int \langle C \rangle \langle (1-C) \rangle \mathrm{d}\kern 0.06em x_3} = 1- 6 \overline{c^2}, \end{equation}

and presented in figure 5(a). Indeed, during the initial phase of elongated fingers, the mixing parameter is small $\varTheta \simeq 0.4$. When the vertical turbulent velocity exceeds the Alfvén velocity and transition to turbulence occurs, the mixing parameter increases, a phenomenology consistent with the observations of Cook, Cabot & Miller (Reference Cook, Cabot and Miller2004). The final values of the mixing parameter for $B_0 \leq 0.2$ are quite close to the non-magnetic case, with $0.76 \leq \varTheta \leq 0.80$, which is a classical range (Gréa Reference Gréa2013). Interestingly, for $B_0=0.3$, the mixing parameter is larger than in the non-magnetic case, with $\varTheta \simeq 0.83$: despite the long transient regime, mixing is eventually more efficient. This can be interpreted as follows: due to the vertical stretching of the structures, the surface available for mixing and mass exchanges at the time of transition to turbulence is greater.

Figure 5.

Effects of the mean magnetic field intensity $B_0$ in the R1 simulations: (a) mixing parameter $\varTheta$ (4.1); (b) global anisotropy indicator $\sin ^2 \gamma$ (4.2).

The directional anisotropy of the concentration field can be analysed through the dimensionality parameter $\sin ^2 \gamma$, expressing the ratio between the amplitudes of the two first spherical harmonics, defined as

(4.2)\begin{equation} \sin^2 \gamma = \frac{\displaystyle \int \mathcal{E}_{cc}(\boldsymbol k) \sin^2 \theta \,\mathrm{d}V_k}{\displaystyle \int \mathcal{E}_{cc}(\boldsymbol k) \,\mathrm{d}V_k}, \end{equation}

where $\mathcal {E}_{cc}$ is the spectral concentration variance such that $L_{x_3} \int \mathcal {E}_{cc} \, \mathrm {d}V_k = L\overline {c^2}$, with $\mathrm {d}V_k$ the infinitesimal spherical volume in spectral space. This parameter takes the value 2/3 for an isotropic field, such that $2/3 \leq \sin ^2 \gamma \leq 1$ for vertically stretched structures and such that $0 \leq \sin ^2 \gamma \leq 2/3$ for flat ones.

This anisotropy indicator is presented in figure 5(b) for various $B_0$. When there is an effective transition to turbulence, i.e. when $\varTheta$ increases in figure 5(a), $\sin ^2 \gamma$ is significantly reduced due to the partial return to isotropy of the smallest scales of the flow. In the asymptotic regime – which is not reached for $B=0.5$ – a stronger $B_0$ results in a larger $\sin ^2 \gamma$, which indeed corresponds to vertically stretched structures. This shows that there is a persistent effect of the mean magnetic field on the Rayleigh–Taylor instability, even after transition to turbulence.

Another way to observe the delayed transition to turbulent mixing in the MRTI is to analyse the time evolution of the concentration probability density function $f(C)$. To make the visualization easier, we restrict our attention to the middle plane $x_3=0$ and investigate how the distribution $f(C)$ evolves with time for $B_0=0$ and $B_0=0.2$ in figures 6(a) and 6(b), respectively. Initially, due to the perturbed interface, there are mostly values around $C \simeq 1/2$. However, when the instability develops, pure heavy fluid falls down while light fluid goes up, showing the strong values of $f(C)$ for $C=0$ and $C=1$ (the maximum value is set at $3$ for readability). Then, around $2 \leq t \leq 3$, turbulent mixing starts for $B_0=0$ and a wide range of intermediate $C$-values appears. For the magnetic case, it takes much longer for turbulent mixing to occur.

Figure 6.

Probability density function $f(C)$ of the concentration field at the plane $x_3=0$ as a function of time: (a) for $B_0=0$ (R1B0G10); (b) for $B_0=0.2$ (R1B020G10).

In figures 7(a) and 7(b), the time evolution of the kinetic energy $\mathcal {K}_{uu}=\overline {u^2}/2$ and magnetic energy $\mathcal {K}_{bb}=\overline {b^2}/2$, normalized by $\dot {L}^2$, are presented. After a transient growth, both ratios strongly decrease: this corresponds to the regime of rapid increase of $L(t)$, associated with the stretching of vertical structures with no mixing. Then, after transition to turbulence, $\mathcal {K}_{uu}/\dot {L}^2$ and $\mathcal {K}_{bb}/\dot {L}^2$ tend towards constant values at late times, except for $B_0=0.5$ where the flow is still in a transient phase. As discussed before, since $L(t)$ grows like $t^2$ in the self-regime of the MRTI, and thus $\dot {L} \sim t$, this shows that both $\mathcal {K}_{uu}$ and $\mathcal {K}_{bb}$ grow also like $t^2$.

Figure 7.

Effects of the mean magnetic field intensity $B_0$ in the R1 simulations: (a) normalized kinetic energy $\mathcal {K}_{uu}$; (b) normalized magnetic energy $\mathcal {K}_{bb}$.

Moreover, while the asymptotic value of $\mathcal {K}_{uu}/\dot {L}^2$ seems quite independent of the vertical mean magnetic field for $0.1 \leq B_0 \leq 0.3$, the level of induced turbulent magnetic energy increases. It will be shown later in § 6 that the ratio $\mathcal {K}_{bb}/\mathcal {K}_{uu}$ is of primary importance to determine the growth rate of the mixing layer in the self-similar state of the MRTI.

Now that the global dynamics has been described, and the effects of a mean magnetic field upon mixing and anisotropy exhibited, we quantify transition to turbulence in the MRTI by investigating the turbulent Reynolds number

(4.3)\begin{equation} Re_T = \frac{\mathcal{K}_{uu}^2}{\nu \varepsilon_{uu}},\end{equation}

where $\mathcal {K}_{uu} = \overline {u_i u_i}/2$ is the kinetic energy and $\varepsilon _{uu} = \nu \overline {\partial _j u_i \partial _j u_i}$ its dissipation rate. In figure 8(a), the transition to turbulence is evidenced in the magnetic configurations by a strong decrease of the Reynolds number, corresponding to an important dissipation of energy. This is correlated to the onset of small-scale shear instabilities that makes $\varTheta$ increase and $\sin ^2\gamma$ decrease. This is quite different from the hydrodynamic case, where $Re_T$ always increases. In addition, transition to turbulence happens for larger $Re_T$ when $B_0$ increases, basically around $Re_T \simeq 135$ for $B_0=0$, up to almost $Re_T = 10^3$ for $B_0=0.3$. Other definitions of the Reynolds number are discussed in Appendix D.2.

Figure 8.

(a) Effects of the mean magnetic field intensity $B_0$ on the turbulent Reynolds number $Re_T$, defined in (4.3), for the R1 simulations. The cross indicates $t_{tr}$ for $B_0=0.3$. (b) Normalized transition time to turbulence $t_{tr}/T_{ref}$ as a function of the magnetic Froude number $F_M$ for all simulations of table 1. Colours refer to the diffusion coefficients: the green circle indicates the R1 simulation with $B_0=0.2$ and $P_M=1$ to be compared with the two simulations R1B020G10P1 and R1B020G10P2 with $P_M<1$.

We now turn our attention to the time $t_{tr}$ at which transition to turbulence occurs. Elongated structures are observed to break after the sudden decrease of $Re_T$ in figure 8(a). Therefore, we choose $t_{tr}$ as the time at which $Re_T$ reaches its minimum value, after the first maximum (for example, around $t\simeq 8$ for $B_0=0.3$). In addition, we define characteristic time, velocity and length scales as

(4.4a)$$\begin{gather} T_{ref} = \left( \frac{\nu}{(\mathcal{A}g)^2} \right)^{1/3} = \left( \frac{\eta}{(\mathcal{A}g)^2} \right)^{1/3} P_M^{1/3}, \end{gather}$$

(4.4b)$$\begin{gather}V_{ref} = ( \mathcal{A}g \eta )^{1/3} = ( \mathcal{A}g \nu )^{1/3} P_M^{{-}1/3}, \end{gather}$$

(4.4c)$$\begin{gather}L_{ref} = V_{ref} T_{ref} = \left( \frac{\eta^2}{\mathcal{A}g} \right)^{1/3} P_M^{1/3}, \end{gather}$$

where $P_M=\nu /\eta$ is the magnetic Prandtl number. These reference scales are chosen independent of the peak wavenumber $k_p$ at which perturbations are injected in the simulations, because we assume that $k_p$ is small enough for the asymptotic state to not depend significantly on it (the effect of $k_p$ is discussed later in figures 13a and 13b). Finally, we define global magnetic Froude and Reynolds numbers, which are related as

(4.5a,b)\begin{equation} F_M = \frac{B_0}{V_{ref}}, \quad R_M = \frac{L_{ref} B_0}{\eta} = P_M^{1/3} F_M.\end{equation}

In figure 8(b), when considering all simulations of table 1, the normalized transition time to turbulence $t_{tr}/T_{ref}$ shows two scalings as function of the magnetic Froude number $F_M$. For $F_M \leq 2$, $t_{tr}$ is relatively independent of $F_M$ because induced magnetic fields are not strong enough to modify significantly the dynamics. While for $F_M \geq 2$, $t_{tr}$ evolves linearly with $F_M$, and this illustrates once again the delayed transition to turbulence when increasing the mean magnetic field intensity. This scaling with $F_M$ is used later on to characterize the self-similar regime of the MRTI.

In figure 8(b), scaling with $R_M$ rather than $F_M$ would have slightly shifted towards the left the two simulations R1B020G10P1 and R1B020G10P2 with $P_M < 1$ (green squares). The black square circled in green indicates simulation R1B020G10 where $P_M=1$ for comparison.

4.2. Energy budget and irreversible mixing

In this section, the emphasis is put on the MRTI energy budget. In the Rayleigh–Taylor instability framework, the variation of potential energy with respect to its initial value corresponds to the energy injected into the system, and is defined as (Cook et al. Reference Cook, Cabot and Miller2004)

(4.6)\begin{equation} E_P=g \overline{x_3(\rho-\rho(0))} = \frac{g}{L} \int (\langle \rho \rangle - \rho(0) ) x_3 \, \mathrm{d}\kern 0.06em x_3 \simeq{-}\frac{\mathcal{A}g \rho_0 L}{12}. \end{equation}

The last equality stands for a linear mean density profile in the mixing zone, which is well verified in the present simulations. The total energy of the system is $E = E_K + E_P$, where $E_K = \rho _0 \mathcal {K}$ is the turbulent energy. The equation for $\mathcal {K}$ (which is either $\mathcal {K}_{uu}$ in the hydrodynamic case or $\mathcal {K}_{uu} + \mathcal {K}_{bb}$ in the MHD case) is given in (A3) and is derived in Appendix A.2: turbulent energy is produced by the mass flux $\mathcal {F}=\overline {u_3 c}$, participates to the growth of the mixing zone ($\mathcal {K}\dot {L}/L$ term) and is finally dissipated at the rate $\varepsilon$. In addition, the mass flux is related to potential energy: indeed, after time derivation of (2.4) and using (2.1c), one ends up with $\dot {L} = -12 \overline {u_3 c}$ (Gréa Reference Gréa2013; Soulard et al. Reference Soulard, Griffond and Gréa2014). Hence, the equation for the total energy reads

(4.7)\begin{equation} \left( \frac{\mathrm{d} }{\mathrm{d} t} + \frac{\dot{L}}{L} \right) (\rho_0 \mathcal{K} + E_P ) ={-} \rho_0 \varepsilon, \end{equation}

which can also be written in a manner comparable to that of Cabot & Cook (Reference Cabot and Cook2006), after time-integration,

(4.8)\begin{equation} LE = \rho_0 L\mathcal{K} + L E_P ={-} \int_0^t \rho_0 L \varepsilon \, \mathrm{d}t' ={-} \varPsi.\end{equation}

We recall at this point that the imposed mean magnetic field does not bring any additional energy to the system: $\rho _0 B_0^2/2$ remains constant throughout the simulation, it only creates a new reservoir of energy through $\mathcal {K}_{bb}$.

The energy budget is displayed in figure 9(a) for both the $B_0=0$ and $B_0=0.3$ cases. Turbulent energy $E_K$ first increases and then decreases at transition to turbulence when dissipation increases significantly due to small-scale shear instabilities. Transition occurs later for $B_0=0.3$, and the proportion of dissipated energy increases by roughly $15\,\%$ compared to the case $B_0=0$. This is consistent with the analysis presented later in § 6, where it will be shown that the global drag of the flow is increased by the presence of turbulent magnetic energy $\mathcal {K}_{bb}$. Note that the sum $(E_K+\varPsi )/|E_P|$ is different from unity at early times because the expression at the right-hand side of (4.6) is not verified: indeed, this relation is true only for a developed mixing layer in which $\langle C \rangle$ is linear (see in figure 11 the profiles of $\langle C \rangle$).

Figure 9.

Total energy budget for the R1 simulations: (a) energy budget (4.8) for $B_0=0.3$ (solid lines) and $B_0=0$ (thin dashed lines); (b) detailed energy budget for the case $B_0=0.3$, with the kinetic (dash-dotted) and magnetic (dashed) contributions of turbulent energy (blue) and dissipation (red). The thin solid lines indicate the sum of the energy and dissipation contributions, as in figure 9(a).

The detailed energy budget for $B_0=0.3$ is presented in figure 9(b): kinetic energy $\mathcal {K}_{uu}$ increases more rapidly than its magnetic counterpart $\mathcal {K}_{bb}$ and represents at late time a more significant part of the injected potential energy. Similar conclusions can be drawn for the dissipation. The magnetic to kinetic energy ratio itself is addressed later on in § 6 where its constancy in the self-similar regime will be used for modelling.

There is a striking feature regarding the conversion of potential energy into turbulent energy: in figure 10(a), the amount of injected potential energy converted into vertical turbulent energy $\rho _0 (\overline {u_3^2} + \overline {b_3^2})$ seems independent of the mean magnetic field intensity $B_0$ (except for $B_0=0.5$ where the asymptotic state is not reached) and is roughly $60\,\%$. Note that the ratio of total turbulent energy on injected potential energy is much more scattered, as seen in figure 9(a), with a decrease from $50\,\%$ ($B_0=0$) to $35\,\%$ ($B_0=0.3$). Therefore, since similar asymptotic states are reached in terms of mixing (see $\varTheta$ in figure 5a), the previous outcome shows that it is dominantly the vertical part of total turbulent energy that drives mixing, the main effect of induced magnetic fields being to decrease the transverse turbulent energy.

Figure 10.

Effects of the mean magnetic field intensity $B_0$ in the R1 simulations. (a) Ratio of total vertical turbulent energy to potential energy $\rho _0 (\overline {u_3^2} + \overline {b_3^2})/|E_p|$. (b) Ratio of the mixing zone length computed from the sorted concentration field $\tilde {L}$ on the usual one $L$. Time is normalized by its final value. Inset shows mixing efficiency $\eta _{ME}$ defined in (4.9) at the final time.

Eventually, only part of the injected potential energy produces irreversible mixing. This can be evaluated by computing the mixing length $\tilde {L}$ based on the sorted concentration field (Winters et al. Reference Winters, Lombard, Riley and D'Asaro1995): the definition of $\tilde {L}$ is given in (D2) along with some details on sorted fields in Appendix D.1. The ratio $\tilde {L}/L$ is shown in figure 10(b) and is highly correlated to $\varTheta$ in figure 5(a). It is remarkable that more irreversible mixing can be obtained in the presence of a strong enough mean magnetic field, for instance, $B_0=0.3$. Nevertheless, it takes longer to reach the asymptotic value of $\tilde {L}/L$ with $B_0=0.3$ than with $B_0=0$. Hence, what is the most efficient configuration to produce irreversible mixing? There is a compromise between the amount of irreversible mixing obtained and the time it takes.

To answer this question, the mixing efficiency is introduced, defined as (Peltier & Caufield Reference Peltier and Caufield2003; Davies Wykes, Hughes & Dalziel Reference Davies Wykes, Hughes and Dalziel2015)

(4.9)\begin{equation} \eta_{ME} = \frac{\displaystyle \int \varepsilon_{bpe} \,\mathrm{d}t}{\displaystyle \int (\varepsilon_{bpe} + \varepsilon) \,\mathrm{d}t}, \end{equation}

where $\varepsilon _{bpe}$ is the dissipation rate of background potential energy, evaluated with the sorted concentration field. This quantity accounts for all the time history of the flow, and its final value is indicated in the inset of figure 10(b) for various $B_0$. Consistent with the discussion of Davies Wykes & Dalziel (Reference Davies Wykes and Dalziel2014), one has $\eta _{ME} \leq 0.5$ whatever the value of $B_0$. Furthermore, it appears that the hydrodynamic case remains the most efficient. Despite having less irreversible mixing, the self-similar state, in which $\tilde {L}/L$ is constant, is reached more rapidly. In comparison, the case $B_0=0.2$ has the same final value for $\tilde {L}/L$, but $\eta _{ME}$ is smaller, because transition to turbulence has been delayed, and so one has to wait longer to reach the same state of irreversible mixing. The case $B_0=0.3$ balances this delay by reaching a larger amount of irreversible mixing. This observation that magnetic configurations are less effective than $B_0=0$ is consistent with the kinetic energy budget of figure 9(a), where dissipation was significantly increased in the presence of a mean magnetic field.

4.3. Vertical profiles

Finally, we investigate the vertical profiles of horizontally averaged quantities. First, we focus on the mean concentration $\langle C \rangle$ in figure 11(a,b), whose extent corresponds roughly to the mixing zone size. At $t=2$, in agreement with figure 4(b), the non-magnetic case has the largest developed mixing zone, whereas the strongest Alfvén velocity $B_0=0.5$ yields the smallest $L$. However, at later time $t=10$, this is for $B_0=0.3$ that the mixing layer is the largest, consistently with what was observed before. Furthermore, it is remarkable that $\langle C \rangle$ remains roughly linear in the turbulent mixing layer despite the presence of a mean magnetic field. A cubic fit would nicely describe $\langle C \rangle$ in the mixing region, provided the virtual origin is increased when $B_0 \neq 0$ since transition to turbulence is delayed (Boffetta, De Lillo & Musacchio Reference Boffetta, De Lillo and Musacchio2010). For $B_0=0.5$ however, the mixing layer is large as well, but the details of the flow are quite different: indeed, $\langle C \rangle$ is more curved because transition to turbulence has not occurred yet at $t=10$.

Figure 11.

Effects of the mean magnetic field intensity $B_0$ on horizontally averaged profiles for the R1 simulations. At (a,c,e,g) $t=2$ and (b,d,f,h) $t=10$, for (a,b) the mean concentration $\langle C \rangle$, (c,d and e,f) the vertical kinetic $\langle u_3^2 \rangle /|E_p|$ and magnetic $\langle b_3^2 \rangle /|E_p|$ energies, respectively, normalized by the injected potential energy, and (g,h) the concentration variance $\langle c^2 \rangle$.

Moving on with the vertical kinetic energy normalized by the released potential energy $\langle u_3^2 \rangle /|E_p|$, with $E_p$ defined in (4.6), it is clear in figure 11(c,d) that at $t=2$, turbulence is damped by the presence of the mean magnetic field, except for $B_0=0.1$: indeed, in figure 4(a), the typical vertical velocity already exceeds the Alfvén one at $t=2$. Later at $t=10$, the profiles of the vertical kinetic energy are much wider. Notably, the extent of $\langle u_3^2 \rangle$ is wider for $B_0=0.3$ and $B_0=0.2$ than for the non-magnetic case, with a flatter profile and a maxima slightly less intense. However, for $B_0=0.5$, kinetic energy grows strongly without mixing and the profile remains smooth.

In figure 11(e,f), the vertical magnetic energy $\langle b_3^2 \rangle$ is quite correlated to $\langle u_3^2 \rangle$ at $t=2$: a stronger $B_0$ results in less fluctuations. The cases $B_0=0.1$ and $B_0=0.2$ are comparable because the former is approaching transition to turbulence, whereas for the latter, $L(t)$ starts growing rapidly. Later at $t=10$, larger mean magnetic field intensities yield stronger induced magnetic energy. More details about the shape and time evolution of the magnetic energy profiles are given in Appendix D.4. Note that for all cases, either at $t=2$ or $t=10$, the magnitude of the vertical magnetic energy is smaller than the vertical kinetic energy.

The discussion is different when it comes to the concentration variance $\langle c^2 \rangle$ in figure 11(g,h). At $t=2$, there are two distinct effects: the peak of $\langle c^2 \rangle$ is greater for $0.1 \leq B_0 \leq 0.3$ than $B_0=0$ since there is less mixing due to the early damping of the RTI, but still a growth of the perturbations. However, for $B_0=0.5$, the interface has hardly evolved, as shown in figure 2. Then, at $t=10$, the concentration variance is significantly reduced. In particular, in the centre of the mixing zone, the concentration variance is smaller for $B_0=0.3$ than $B_0=0$, once again showing that this configuration is more efficient at mixing the two fluids, as already revealed in figure 5(a). However, the maximum concentration variance remains large for $B_0=0.5$ because the strong magnetic field still inhibits small-scale shear instabilities and thus mixing.

5.1. Dynamics in physical space

First, vertical slices of the turbulent flow of run R4B0G10b are presented in figure 12 at times $t=3.5$ and $t=10$. A wide range of structures of various scales and shapes is observed. One can see, for example, small patches of light fluid surrounded by heavy fluid, sharp gradients (in black) between light and heavy fluids, and large patches of mixed fluids that penetrate the reservoirs of pure fluids (dark blue and dark red).

Figure 12.

Non-magnetic case R4B0G10b. Vertical slice of the concentration field at (a) $t=3.5$ and (b) $t=10$, with the same colourmap as figure 2, except for the pure fluids which are not transparent any longer. The isotropic and vertical Taylor Reynolds numbers at $t=10$ are $Re_\lambda = 210$ and $Re_\lambda ^{(z)} = 314$.

The ‘isotropic’ Taylor Reynolds number $Re_\lambda = \sqrt {20 \mathcal {K}_{uu}^2/(3 \nu \varepsilon _{uu})}$ is $Re_\lambda = 210$. The Reynolds number based solely on the vertical velocity is $Re_\lambda ^{(z)} = 314$. The ratio $Re_\lambda ^{(z)}/Re_\lambda ^{(x)}$ is of order $\simeq 1.6$, somehow smaller than the values reported by Cook & Dimotakis (Reference Cook and Dimotakis2001) for the RTI at larger Atwood numbers.

We are now interested in the main integrated quantities useful to characterize anisotropy and mixing, namely the mixing parameter $\varTheta$ defined in (4.1), the global anisotropy parameter $\sin ^2 \gamma$ defined in (4.2), the ratio of the vertical kinetic energy $R_{33} = \overline {u_3^2}$ to the transverse one $0.5(R_{11}+R_{22})$ and the effective Atwood number $\mathcal {A}_e$ defined by Cook et al. (Reference Cook, Cabot and Miller2004):

(5.1)\begin{equation} \mathcal{A}_e = \frac{\sqrt{\langle {\rho'}^2 \rangle}|_0}{\langle \rho \rangle|_0} = \frac{2\mathcal{A} \sqrt{\langle {c}^2 \rangle}|_0}{1+2\mathcal{A}(\langle C \rangle |_0 -1/2)},\end{equation}

where $\rho '$ is the fluctuating density and $\langle {\cdot } \rangle |_0$ refers to the mid plane ($x_3=0$) value of the horizontally averaged profile.

All these quantities are presented in figure 13(a) for the two non-magnetic cases R4B0G10a and R2B0G10b, which have different peak wavenumber $k_p$ but the same diffusion coefficients $\nu =\kappa =10^{-4}$. For all quantities, the transient regime is slightly shorter for R4B0G10a, which has the largest initial peak wavenumber ($k_p=50$).

Figure 13.

Non-magnetic cases R4B0G10a ($4096^3$, $\nu =10^{-4}$, $k_p=50$, solid lines) and R2B0G10b ($2048^3$, $\nu =10^{-4}$, $k_p=40$, dash-dotted lines). (a) Global anisotropy $\sin ^2 \gamma$ (black), mixing parameter $\varTheta$ (blue), horizontal to vertical kinetic energy ratio (red) and normalized effective Atwood number $\mathcal {A}_e/\mathcal {A}$ (green). (b) Growth rate $\alpha$ of the turbulent mixing zone (5.2) for simulation R4B0G10a, computed from DNS (black, $\dot {L}^2/(8 \mathcal {A}g L)$) and from theoretical prediction (5.3) (blue). The inset shows for comparison the $\alpha$ from simulations R2B0G10b (black dash-dotted line) and R4B0G10b (red solid line).

The final value of the mixing parameter for R4B0G10a (and R4B0G10b) is $\varTheta =0.78$, slightly lower than cases R2B0G10b and R1B0G10 (see figure 5a) which both have $k_p=40$ and reach $\varTheta =0.8$. In addition, the effective Atwood number $\mathcal {A}_e$ is larger for R4B0G10a than R2B0G10b, showing that there has been less mixing for $k_p=50$, which is consistent with the smaller value of $\varTheta$ for R4B0G10a.

The global anisotropy parameters are quite similar, $\sin ^2 \gamma \simeq 0.725$, corresponding to slightly vertically elongated structures. Anisotropy can also be observed through the ratio $0.5(R_{11}+R_{22})/R_{33}$ which tends to $0.30$, a value consistent with those reported by Dimonte et al. (Reference Dimonte2004).

In the self-similar regime, the mixing zone size evolves as (Youngs Reference Youngs1994)

(5.2)\begin{equation} L(t) = 2 \alpha \mathcal{A} g t^2,\end{equation}

where $\alpha$ is the growth rate of the ‘bubbles’. Since we are within the Boussinesq approximation, the growth of the turbulent mixing zone is symmetric between the light and heavy fluids sides. The growth rate is evaluated using the relation $\alpha = \dot {L}^2/(8 \mathcal {A}g L)$ and displayed in figure 13(b); its value at $t=11$ is $\alpha =0.021$ for case R4B0G10a, which is comparable to other measured values in numerical simulations (Dimonte et al. Reference Dimonte2004). The theoretical prediction (Gréa Reference Gréa2013; Soulard, Griffond & Gréa Reference Soulard, Griffond and Gréa2016)

(5.3)\begin{equation} \alpha_{hyd} = \frac{\sin^4 \gamma (1-\varTheta)^2}{1+ \sin^2 \gamma (1-\varTheta)}, \end{equation}

which relates the growth rate to the mixing and anisotropy parameters, is also presented in the same figure. This is in excellent agreement with the numerical result. Parametric investigations of the transition to turbulence in the RTI later confirmed this scaling for a wide variety of initial conditions (Morgan & Black Reference Morgan and Black2020).

Still in figure 13(b), the growth rate of the mixing zone for R2B0G10b is displayed in the inset (dash-dotted line). The growth rate is slightly smaller, $\alpha =0.017$, due to the stronger molecular mixing that slows down the enlargement of the turbulent mixing zone. This is compatible with the smaller effective density contrast observed in figure 13(a) for R2B0G10b (Cook et al. Reference Cook, Cabot and Miller2004). The growth rate for R4B0G10b is presented as well (red solid line) and reaches $\alpha =0.022$ at $t=11$, which shows that the asymptotic turbulent state of the RTI is quite independent of diffusion coefficients, provided they are small enough, at least for this initial condition.

5.2. Spectral scale-by-scale analysis

We continue by analysing the main turbulent spectra for the RTI using simulation R4B0G10b. In the following, we note with $\hat {{\cdot }}$ the Fourier transforms of the fluctuating fields. We define the kinetic energy $E_{uu}$, concentration variance $E_{cc}$ and vertical flux $E_{uc}$ spectra as averages of the spectral densities $\mathcal {E}_{xx}$ on spherical shells $S_k$ of radius $k$:

(5.4a)$$\begin{gather} E_{uu}(k) = \int_{S_k} \mathcal{E}_{uu}(\boldsymbol k)\, {\rm d}^2\boldsymbol k \quad \text{with } \mathcal{E}_{uu}(\boldsymbol k) = \hat{u}_i(\boldsymbol k) \hat{u}_i(-\boldsymbol k)/2, \end{gather}$$

(5.4b)$$\begin{gather}E_{cc}(k) = \int_{S_k} \mathcal{E}_{cc}(\boldsymbol k)\, {\rm d}^2\boldsymbol k \quad \text{with } \mathcal{E}_{cc}(\boldsymbol k) = \hat{c}(\boldsymbol k) \hat{c}(-\boldsymbol k), \end{gather}$$

(5.4c)$$\begin{gather}E_{uc}(k) = \int_{S_k} \mathcal{E}_{uc}(\boldsymbol k)\, {\rm d}^2\boldsymbol k \quad \text{with } \mathcal{E}_{uc}(\boldsymbol k) = \hat{c}(\boldsymbol k) \hat{u}_3(-\boldsymbol k). \end{gather}$$

These three spectra are presented in figure 14(a) at the latest time for R4B0G10b, revealing a wide inertial range spanning more than one decade. For $E_{uu}$ and $E_{cc}$, the inertial scaling is close to $k^{-5/3}$, whereas it is steeper for the flux $E_{uc}$. The Kolmogorov wavenumber $k_\eta =(\varepsilon _{uu}/\nu ^3)^{1/4}$ is indicated as well, showing that the simulation is well resolved. The values of $k_{max}/k_\eta$ at the final time for all simulations are given in table 1. The time evolution of $E_{uu}$ is shown in figure 14(b) from $t=0.5$ up to $t=11$, showing the dynamical development of the inertial range, along with the enlargement of the integral scale. The infrared slope at large scales remains at its initial value $k^4$.

Figure 14.

Main spectra at $t=11$ ($Re_\lambda = 235$) for the non-magnetic case R4B0G10b. (a) Kinetic energy $E_{uu}$, scalar variance $E_{cc}$ and vertical flux $E_{uc}$ spectra; the vertical dashed line indicates the Kolmogorov wavenumber $k_\eta$. (b) Time evolution of $E_{uu}$ from $t=0.1$ to $t=11$. (c) Compensated spectra with Kolmogorov (5.5a), Corrsin–Obukhov (5.5b) and Lumley scalings (5.5c): note that for the latter, the intensity is divided by four for clarity.

We now consider the compensated spectra to assess the inertial scalings

(5.5a)$$\begin{gather} E_{uu}(k) = C_K \varepsilon_{uu}^{2/3} k^{{-}5/3}, \end{gather}$$

(5.5b)$$\begin{gather}E_{cc}(k) = C_{CO} \varepsilon_{cc} \varepsilon_{uu}^{{-}1/3} k^{{-}5/3}, \end{gather}$$

(5.5c)$$\begin{gather}E_{uc}(k) = C_F \varepsilon_{uu}^{1/3} \varLambda k^{{-}7/3}, \end{gather}$$

where $\varepsilon _{uu}=2\nu \int k^2 E_{uu}\, {\rm d}k$ is the kinetic energy dissipation rate, $\varepsilon _{cc}=2\kappa \int k^2 E_{cc} \,{\rm d}k$ is the scalar variance dissipation rate and $\varLambda =1/L$ is the mean concentration gradient. The two first scalings (5.5a) and (5.5b) constitute the classical Kolmogorov–Obukhov theory derived for isotropic turbulence, whereas (5.5c) was derived by Lumley (Reference Lumley1967) for the anisotropic flux spectrum.

In figure 14(c), the Kolmogorov–Obukhov scalings for $E_{uu}$ and $E_{cc}$ are reasonably well assessed, with a plateau for the compensated spectra spanning roughly one decade. These inertial scalings are in agreement with results obtained by Boffetta et al. (Reference Boffetta, Mazzino, Musacchio and Vozella2009), Soulard & Griffond (Reference Soulard and Griffond2012) and with the theoretical predictions derived from Soulard (Reference Soulard2012), based on a Monin–Yaglom relation adapted to the framework of the turbulent RTI. Regarding the spectral flux, the inertial range slope is shallower than $k^{-7/3}$, closer to $k^{-2}$, but larger Reynolds numbers are required to conclude on this matter (Briard, Gomez & Cambon Reference Briard, Gomez and Cambon2016; Soulard & Gréa Reference Soulard and Gréa2017).

Despite a significant bottleneck, the Kolmogorov constant $C_K \simeq 1.5$ is recovered, along with the Corrsin–Obukhov one, $C_{CO} \simeq 0.65$ (Sreenivasan Reference Sreenivasan1996). The constant related to the vertical flux is $C_F \simeq 4$, compatible with previous measurements in isotropic turbulence (O'Gorman & Pullin Reference O'Gorman and Pullin2005; Briard et al. Reference Briard, Gomez and Cambon2016).

We pursue the spectral analysis by evaluating the different terms in the equations for the kinetic energy and scalar variance spectra:

(5.6a)$$\begin{gather} \left( \frac{\partial}{\partial t} + 2 \nu k^2 \right) E_{uu}(k) = T_{uu}^{(u)}(k) + P_{uu}(k), \end{gather}$$

(5.6b)$$\begin{gather}\left( \frac{\partial}{\partial t} + 2 \kappa k^2 \right) E_{cc}(k) = T_{cc}(k) + P_{cc}(k), \end{gather}$$

where $T_{uu}^{(u)}$ and $T_{cc}$ are the spherically averaged nonlinear terms for kinetic energy and concentration variance, respectively, whose expressions are given in Appendix A.3. The production terms are $P_{uu} = - 2 \mathcal {A}g E_{uc}$ and

(5.7)\begin{equation} P_{cc}= \int_{S_k} [\hat{c}(\boldsymbol k) \mathcal{P}_c(-\boldsymbol k) + \hat{c}(-\boldsymbol k) \mathcal{P}_c(\boldsymbol k) ]\, \mathrm{d}^2 \boldsymbol k, \quad \text{with } \mathcal{P}_c(\boldsymbol k) ={-} (\hat{u}_3 * \widehat{\partial_{x_3} \langle C \rangle})(\boldsymbol k), \end{equation}

with $*$ the convolution product. Later on in § 6.4, this production term will be approximated by $P_{cc} \simeq -2 E_{uc}/L$. For conciseness hereafter, the dissipative terms are written by $D_{uu}=2\nu k^2 E_{uu}$ and $D_{cc}=2\kappa k^2 E_{cc}$.

These contributions are presented in figures 15(a) and 15(b). We observe a direct transfer of kinetic energy and concentration variance from large to small scales, a balance between dissipation and nonlinear transfers at the smallest scales, whereas there is an excess of production through the mass flux at the largest scales, allowing the kinetic energy and scalar variance to grow with time. The fluxes $\varPi _{xx} = - \int _0^k T_{xx}(p) \,\mathrm {d}p$ normalized by the corresponding dissipation rates are presented as well as insets, showing the conservative nature of the nonlinear transfers.

Figure 15.

Nonlinear transfers, production and dissipative terms at $t=11$ for the non-magnetic case R4B0G10b (a) for the kinetic energy and (b) for the concentration variance. The insets show the normalized fluxes.

6.1. Elongated structures and anisotropy

First, a vertical slice of the turbulent flow from simulation R4B02G10b at two different times is presented in figure 16, showing once again the striking difference in terms of structures compared with the non-magnetic case in figure 12.

Figure 16.

Case R4B02G10b. Vertical slice of the concentration field, with the same colourmap as figure 12, at two different times: (a) $t=3.5$ and (b) $t=10$, where the isotropic, vertical and total Taylor Reynolds numbers are respectively $Re_\lambda = 177$, $Re_\lambda ^{(z)} = 245$ and $Re_\lambda ^{(t)} = 189$.

The smooth elongated fingers are quite narrow and eventually destabilize under the effects of shear instabilities, giving birth to a wide variety of scales. Nevertheless, compared with the non-magnetic case, one can still see an imprint of the initial elongated fingers in the fully turbulent regime, with vertically stretched structures. The Taylor Reynolds numbers, for the second image at $t=10$, are $Re_\lambda = 177$ for the isotropic version and $Re_\lambda ^{(z)} = 245$ for the vertical component. One can also use a total Taylor Reynolds number, based on the total energy $\mathcal {K}= \mathcal {K}_{uu}+\mathcal {K}_{bb}$ and its dissipation $\varepsilon$ rather than the kinetic energy only: this yields $Re_\lambda ^{(t)} = 189$; see Appendix D for the effects of $B_0$ on the various Reynolds numbers.

By a simple visual comparison of figures 12 and 16 at $t=10$, it follows that the mixing zone is much larger for the magnetic case, but the overall flow is less turbulent, as confirmed by the lower $Re_\lambda$.

We now turn our attention to global statistics, like for the non-magnetic case, to quantify anisotropy and mixing. The two simulations R4B02G10a and R4B02G10b are presented in figure 17(a) to show the effects of decreasing the diffusion coefficients from $\nu = \kappa = \eta = 10^{-4}$ (R4B02G10a) to $5.10^{-5}$ (R4B02G10b). Note that the peak wavenumber in the R4 simulations ($k_p=50$) is well below the diffusivity cutoff $k_c=(\mathcal {A}g/(2\nu ^2))^{1/3} \simeq 184\unicode{x2013}292$.

Figure 17.

Magnetic cases R4B02G10b ($\nu =5.10^{-5}$, solid lines) and R4B02G10a ($\nu =1.10^{-4}$, dash-dotted lines). (a) Global anisotropy $\sin ^2 \gamma$ (black), mixing parameter $\varTheta$ (blue), horizontal to vertical kinetic energy ratio (red), and normalized effective Atwood number $\mathcal {A}_e/\mathcal {A}$. (b) Growth rate $\alpha$ of the turbulent mixing zone for simulation R4B02G10b, computed from DNS (black, $\dot {L}^2/(8 \mathcal {A}g L)$) and from theoretical prediction (5.3) (blue). Simulation R4B02G10a is shown as well (dash-dotted line).

There is not much difference between the two cases for the directional anisotropy of the scalar field $\sin ^2 \gamma$ on the one hand, and for the anisotropy of the velocity field $0.5(R_{11}+R_{22})/R_{33}$ on the other one. Still, one may interpret the slightly delayed transition to turbulence for the case with larger diffusion coefficients as a consequence of secondary shearing instabilities being quite smoothed. Regarding quantities more sensitive to mixing like $\varTheta$ and $\mathcal {A}_e$, there are some differences in the transient regime. Similar observations were made for the effect of viscosity on the purely hydrodynamic case in Appendix D.3.

If we now evaluate the global effect of a vertical mean magnetic field by comparing the anisotropy and mixing indicators of R4B0G10b and R4B02G10b, there is a remarkable outcome. While mixing is slightly affected by the mean magnetic field – $\varTheta \simeq 0.78$ for $B_0=0$ and $\simeq 0.80$ for $B_0=0.2$ – anisotropy is quite different: the more elongated structures in figure 16 compared with figure 12 are directly reflected by a larger value of $\sin ^2 \gamma$, namely $\simeq 0.81$ for R4B02G10b against $\simeq 0.71$ for R4B0G10b. Consistently, the vertical component of the kinetic energy is roughly three times stronger for $B_0=0.2$ than in the non-magnetic case.

In figure 17(b), we compute the growth rate of the mixing zone using $\alpha = \dot {L}^2/(8 \mathcal {A}g L)$, assuming that the self-similar regime is still defined by (5.2). There is now a strong difference with the prediction (5.3) both in the transient regime and in the asymptotic case: the reasons for the discrepancy with $\alpha _{hyd}$ are addressed later in § 6.4.

Compared with the non-magnetic case in figure 13(b), the growth rate remains much longer at its maximal value, which means that the mixing layer grows much faster in the transient regime with an imposed vertical mean magnetic field. Then, the growth rate settles around a smaller value, $\alpha =0.0126$, confirming that the enhanced mixing due to the breaking of stretched structures slows down the overall growth.

6.2. Spectral scale-by-scale analysis

In the following, we pursue the analysis of the MRTI by a spectral scale-by-scale investigation of transfer and production terms. We now define the magnetic energy spectrum

(6.1)\begin{equation} E_{bb}(k) = \int_{S_k} \mathcal{E}_{bb}(\boldsymbol k)\, {\rm d}^2\boldsymbol k \quad \text{with } \mathcal{E}_{bb}(\boldsymbol k) = \hat{b}_i(\boldsymbol k) \hat{b}_i(-\boldsymbol k)/2, \end{equation}

along with the total energy spectrum $E = E_{uu}+E_{bb}$. The main spectra in figure 18(a), extracted from simulation R4B02G10b at $t=11$, are quite different from the spectra of the hydrodynamic case in figure 18(b). Considering for instance $E_{uu}$, the inertial range is narrower (consistent with a smaller Taylor Reynolds number), and the integral scale is much smaller, which is a consequence of the elongated structures observed in figure 16 which have smaller horizontal extent. The same remark applies for both $E_{cc}$ and $E_{uc}$. In addition, note that the mass flux $E_{uc}$ is much closer to $E_{uu}$ in magnitude than in the non-magnetic case, due to the strong vertical velocity.

Figure 18.

Main spectra at $t=10$ ($Re_\lambda = 191$) for case R4B02G10b. (a) Kinetic energy $E_{uu}$, concentration variance $E_{cc}$, vertical flux $E_{uc}$ and magnetic energy $E_{bb}$ spectra; the vertical dashed lines indicate the Kolmogorov wavenumber $k_\eta$. The inset shows the time evolution of $E_{bb}$. (b) Comparison of the three spectra $E_{uu}$, $E_{cc}$ and $E_{uc}$ at $t=11$ between R4B02G10b (solid lines) and R4B0G10b (dash-dotted lines). (c) Compensated spectra for R4B02G10b: $E_{uu}$ with scaling (6.2a), $E_{bb}$ with scaling (6.2b) and $E_{cc}$ with scaling (6.4). The thin dash-dotted line indicates scaling (5.5b) for $E_{cc}$.

Also of interest is the relative magnitude of the kinetic and magnetic energy spectra. In the literature of strong homogeneous MHD turbulence (SHMHD), sustained with a vertical mean magnetic field, one has in general a residual energy spectrum such that $E_r = E_{uu} - E_{bb} \leq 0$, meaning an excess of magnetic energy at all scales (Muller & Grappin Reference Muller and Grappin2005; Boldyrev et al. Reference Boldyrev, Perez, Borovsky and Podesta2011; Beresnyak Reference Beresnyak2014). Here, in the MRTI, the phenomenology is different since there is an additional large-scale production term, which is the conversion of potential energy into kinetic energy through the mass flux $E_{uc}$. As a result, kinetic energy exceeds magnetic energy at large scales in figure 18(a), and then a negative residual energy spectrum $E_r$ is recovered for smaller scales.

In SHMHD, there exists a controversy regarding the inertial scaling of the total energy spectrum (see discussion by Briard & Gomez Reference Briard and Gomez2018): either $E \sim k^{-3/2}$ (Muller & Grappin Reference Muller and Grappin2005; Boldyrev et al. Reference Boldyrev, Perez, Borovsky and Podesta2011; Mason et al. Reference Mason, Perez, Boldyrev and Cattaneo2012; Perez et al. Reference Perez, Mason, Boldyrev and Cattaneo2012), the so-called Iroshnikov–Kraichnan scaling (Iroshnikov Reference Iroshnikov1964; Kraichnan Reference Kraichnan1965); or $E \sim k^{-5/3}$ (Beresnyak Reference Beresnyak2014) following the theory of Goldreich & Sridhar (Reference Goldreich and Sridhar1995) for anisotropic turbulence. Both scaling were obtained numerically over the years, with the common feature that the residual spectrum $E_r$ is negative, no matter its scaling.

Since it might be difficult to discriminate between two scalings like $k^{-5/3}$ and $k^{-3/2}$ at the limited Reynolds number reached here, we propose in figure 18(c) to analyse fully compensated spectra. Thus, a very large or small value of the plateau would indicate that the scaling is erroneous. Following the observations made in figure 18(a), where $E_{bb}$ is shallower than $E_{uu}$ and exceeds it at small scales, we propose the following inertial scalings:

(6.2a)$$\begin{gather} E_{uu}(k) = C_K {\varepsilon}^{2/3} k^{{-}5/3}, \end{gather}$$

(6.2b)$$\begin{gather}E_{bb}(k) = C_{IK} \sqrt{\varepsilon B_0} k^{{-}3/2}, \end{gather}$$

the latter being the Iroshnikov–Kraichnan scaling applied only to $E_{bb}$. Scaling (6.2a) differs only from scaling (5.5a) by the use of the total dissipation rate $\varepsilon = \varepsilon _{uu} + \varepsilon _{bb}$ rather than $\varepsilon _{uu}$ only. These two scalings seem appropriate in figure 18(c): the Kolmogorov constant in MRTI $C_K \simeq 1.7$ is slightly larger than in the non-magnetic case ($C_K \simeq 1.5$) and the Iroshnikov–Kraichnan constant is $C_{IK} \simeq 1$. It is not possible to discriminate whether the total energy spectrum $E$ scales like $k^{-3/2}$ or $k^{-5/3}$, since $E_{uu}$ dominates at the largest scales of the inertial range and then $E_{bb}$ at smaller ones.

The scalings (6.2a) and (6.2b) differ from numerical simulations of SHMHD by Alexakis (Reference Alexakis2013), where $E_{uu} \sim k^{-3/2}$ and $E_{bb} \sim k^{-5/3}$: these latter scalings are compatible with solar wind data (Boldyrev et al. Reference Boldyrev, Perez, Borovsky and Podesta2011). The differences are not surprising since in the present case, turbulence is unstably stratified, which creates a strong production term at large scales, which is quite different from the SHMHD framework.

If we now consider the concentration variance spectrum $E_{cc}$, its inertial scaling is shallower than $k^{-5/3}$ and deviates from the Corrsin–Obukhov scaling (5.5b), as revealed in figure 18(c) by the dash-dotted line. An alternative scaling can be proposed using classical arguments. Assuming a constant transfer rate of scalar variance at scale $1/k$, one may write

(6.3)\begin{equation} \varepsilon_{cc} = \frac{k E_{cc}(k)}{\tau_{tr}(k)}, \end{equation}

where $\tau _{tr}$ is the transfer time. In hydrodynamics, one has $\tau _{tr} = \tau _V = (k^2 \varepsilon _{uu})^{-1/3}$ and this would yield (5.5b). However, in MHD, the transfer time is reduced due to the propagation and interaction of Alfvén waves, so that $\tau _{tr} = \tau _V^2/\tau _A$, with $\tau _A = (k B_0)^{-1}$ (Galtier, Politano & Pouquet Reference Galtier, Politano and Pouquet1997; Zhou, Matthaeus & Dmitruk Reference Zhou, Matthaeus and Dmitruk2004). Hence, the inertial scaling of the scalar variance spectrum reads

(6.4)\begin{equation} E_{cc}(k) = \tilde{C}_{CO} \varepsilon_{cc} \varepsilon_{uu}^{{-}2/3} B_0 k^{{-}4/3}, \end{equation}

where $\tilde {C}_{CO}$ is a new Corrsin–Obukhov constant. This scaling is assessed in figure 18(c) and the new constant obtained is very close to the previous one in the non-magnetic case, $\tilde {C}_{CO} \simeq C_{CO}$. Whether this scaling would hold at larger Reynolds numbers is an open question. Still, the $k^{-4/3}$ inertial slope reflects that the transfer towards small scales is less effective in MRTI, probably because structures are collimated along the vertical magnetic field lines.

We pursue the spectral analysis by considering the equations for the kinetic, magnetic and scalar spectra,

(6.5a)$$\begin{gather} (\partial_t + 2 \nu k^2 ) E_{uu}(k) = T_{uu}^{(u)}(k) + T_{uu}^{(\textit{b})}(k) + P_{uu}(k) + P_{uu}^{(\textit{b})}(k), \end{gather}$$

(6.5b)$$\begin{gather}(\partial_t + 2 \eta k^2) E_{bb}(k) = T_{bb}(k) + P_{bb}(k), \end{gather}$$

(6.5c)$$\begin{gather}(\partial_t + 2 \kappa k^2 ) E_{cc}(k) = T_{cc}(k) + P_{cc}(k), \end{gather}$$

where the new production term related to $B_0$ is

(6.6)\begin{equation} P_{bb}(k) ={-}P_{uu}^{(\textit{b})}(k) = B_0 \int_{S_k} k \cos\theta \, {\rm Im}[\mathcal{E}_{ub}(\boldsymbol k)] \,\mathrm{d}^2 \boldsymbol k, \end{equation}

with ${\rm Im} [{\cdot }]$ the imaginary part and $\mathcal {E}_{ub} = \hat {u}_i(-\boldsymbol k) \hat {b}_i(\boldsymbol k)$ the spectral density of cross-helicity. The new nonlinear transfer $T_{uu}^{(\textit {b})}$ reflects the effects of the Lorentz force on the velocity field, and is such that $T_{uu}^{(\textit {b})} + T_{bb}$ is a conservative term.

All these terms are presented in figures 19(a)–19(d). For the concentration variance spectrum equation in figure 19(b), the scale-by-scale budget is quite similar to the non-magnetic case in figure 15(b), except that due to the lower Reynolds number, there is a much wider region of the wavenumber space where production ($P_{cc}$) and dissipation ($D_{cc}$) coexist.

Figure 19.

Nonlinear transfers, production and dissipative terms at $t=10$ for simulation R4B02G10b. (a) Kinetic energy. (b) Concentration variance. (c) Magnetic energy. (d) Total energy.

Comparing now the redistribution terms $P_{uu}^{(\textit {b})}$ and $P_{bb}$, they have opposite roles for $E_{uu}$ and $E_{bb}$ in figures 19(a) and 19(c), respectively. Indeed, they show that kinetic energy is converted to magnetic energy at large scales, whereas the opposite occurs at small scales. Regarding nonlinear terms, $T_{uu}^{(u)}$ transfers kinetic energy from large to small scales, as usual. Whereas the effect of the Lorentz force $T_{uu}^{(b)}$ is to extract kinetic energy on a wide range of scales. This can be understood as follows: at the onset of the instability, the magnetic field lines remain naturally along $B_0$. Then, when turbulence becomes strong enough, the magnetic tension acts as a damping of kinetic energy to prevent the stretching of field lines. Regarding the magnetic energy, it is transferred from large to small scales by $T_{bb}$.

Considering now the total energy in figure 19(d), we write $T=T_{uu}^{(u)} + T_{uu}^{(\textit {b})} + T_{bb}$ and observe that it cascades from large to small scales, in a fashion much like kinetic energy in the non-magnetic case (see figure 15a), and the only production term remaining is $P=P_{uu}$, namely production through the mass flux.

6.3. Scale-by-scale anisotropy and self-similarity

After presenting the energy spectra and scale-by-scale energy budgets in § 5.2 for the hydrodynamic case and in § 6.2 for the magnetic one, we now briefly investigate anisotropy and self-similarity of large scales.

The spatial scalar anisotropy indicator $\sin ^2\gamma (k)$ (corresponding to $\sin ^2\gamma$ defined in (4.2) without integration over $k$) is shown in figures 20(a) and 20(b) at $t=10$ for $B_0=0$ and $B_0=0.2$, respectively. For both configurations, $\sin ^2\gamma (k)$ is the strongest at the largest scales ($k< k_{uu}$), indicating vertically elongated structures. However, in the inertial range ($k_{uu} < k < k_\eta$), there is a significant small-scale anisotropy in the magnetic case, since $\sin ^2\gamma (k)$ is much larger than the isotropic value $2/3$ reached in the hydrodynamic case. The transverse to vertical kinetic energy ratio $E_\perp /E_{33}$, with $E_\perp = (E_{11}+E_{22})/2$, illustrates a similar feature, with $E_\perp /E_{33}$ tending towards the isotropic value $1$ for $B_0=0$ at small scales, while remaining very low for $B_0=0.2$. These two indicators show that the imposed vertical mean magnetic field induces a strong persistent small-scale anisotropy.

Figure 20.

Scale-by-scale anisotropy at $t=10$ for simulations (a) R4B0G10b with $B_0=0$ and (b) R4B02G10b with $B_0=0.2$. Scalar anisotropy parameter $\sin ^2\gamma (k)$, magnetic to kinetic energy ratio $E_{bb}/E_{uu}$ and transverse to vertical kinetic energy ratio $E_\perp /E_{33}$. The Kolmogorov wavenumber $k_\eta$ and the wavenumber $k_{uu}$ of the maximum of $E_{uu}$ are indicated as thin vertical dashed lines.

Then, the magnetic to kinetic energy ratio $E_{bb}/E_{uu}$ appears to be roughly constant at large scales in figure 20(b), a feature that will be used in the next part to characterize the asymptotic self-similar state. For the smallest scales, magnetic energy exceeds kinetic energy, as already revealed in figure 18(a).

We conclude this part with considerations about the self-similarity of the largest scales of the flow. We use the following self-similar scaling $E_{uu}/(V_N^2 L^3)$ (Soulard, Griffond & Gréa Reference Soulard, Griffond and Gréa2015), where $V_N=LN$ is the reference velocity with $N=\sqrt {2Ag/L}$ the buoyancy frequency. With this scaling, the normalized kinetic energy spectra at various times collapse at the largest scales of the flow, both for $B_0=0$ and $B_0=0.2$, in figures 21(a) and 21(b), respectively. Furthermore, we recover two features previously mentioned. First, the self-similar scaling is obtained much more rapidly in the hydrodynamic configuration, because the presence of a mean magnetic field delays transition to turbulence. Second, $E_{uu}$ is less intense and peaked at smaller scales for $B_0=0.2$ compared to $B_0=0$, as already seen in figure 18(b).

Figure 21.

Large-scale self-similar scaling for the kinetic energy spectra $E_{uu}/(V_N^2 L^3)$ at $t=0 \ldots 10$ for (a) R4B0G10b with $B_0=0$ and (b) R4B02G10b with $B_0=0.2$.

6.4. The asymptotic state of the magnetic Rayleigh–Taylor instability

One of the main observations from the two previous parts is that the growth rate of the mixing zone can be significantly damped by the presence of the mean magnetic field: it decreases from $\alpha = 0.022$ in the hydrodynamic case to $\alpha =0.012$ for $B_0=0.20$, as revealed in figure 17(b) for R4B02G10b. As a consequence, the prediction (5.3) fails and overestimates the asymptotic growth rate value. This is not intuitive, since the initial more rapid growth of the mixing layer – due to the inhibition of small-scale shear instabilities – yields larger values for $L(t)$ in the asymptotic regime, as observed in figure 4(b). The objective of what follows is to understand and predict this damping of the growth rate in the MRTI.

6.4.1. Background: buoyancy-drag equation

The outcome that the mean vertical magnetic field has a persistent imprint on the asymptotic fully turbulent regime of the Rayleigh–Taylor instability is an important result, and needs now to be interpreted and modelled. First, we can consider the buoyancy-drag equation (Davies & Taylor Reference Davies and Taylor1950; Ramshaw Reference Ramshaw1998; Dimonte Reference Dimonte2000)

(6.7)\begin{equation} \ddot{L} ={-} C_d \frac{\dot{L}^2}{L} + C_b \mathcal{A}g,\end{equation}

where $C_d$ and $C_b$ are respectively the drag and buoyancy coefficients. Assuming a self-similar growth of the mixing zone width (5.2), one gets for the growth rate

(6.8)\begin{equation} \alpha_{bd} = \frac{C_b}{4 + 8 C_d}.\end{equation}

Hence, one could argue that to reduce $\alpha _{bd}$, where ‘$bd$’ stands for buoyancy-drag, the mean magnetic field either increases the drag coefficient or decreases the buoyancy one. To better understand how the large-scale dynamics is modified in the MRTI and how the presence of a vertical mean magnetic field impacts the growth rate, we consider the Rapid Acceleration (RA) model, established by Gréa (Reference Gréa2013) for the hydrodynamic Rayleigh–Taylor instability. For the sake of clarity, we first recall the main ingredients that will be useful to extend the approach to the MRTI.

Within the RA framework, nonlinearities are discarded along with dissipative terms, and only interactions between the fluctuations and the mean field are kept. After some algebra, for the classical Rayleigh–Taylor instability, the buoyancy-drag equation (6.7) takes the form

(6.9)\begin{equation} \ddot{L} ={-} \underbrace{\frac{1}{2}}_{C_d} \frac{\dot{L}^2}{L} + \mathcal{A}g \underbrace{4 \sin^2 \gamma (1-\varTheta)}_{C_b}. \end{equation}

It is remarkable that this equation derives from a Lagrangian formulation, and it was noted by Gréa (Reference Gréa2013), following Ramshaw (Reference Ramshaw1998), that additional dissipative effects could modify the drag coefficient as $C_d = \mathcal {D}+1/2$. This additional turbulent dissipation $\mathcal {D}$ that occurs within the RT instability was modelled by $C_d C_b =2$, as demonstrated by Poujade & Peybernes (Reference Poujade and Peybernes2010). Combining (6.8) and (6.9) along with $C_d C_b =2$ recovers the hydrodynamic prediction (5.3).

6.4.2. Simplified equations for the MRTI at large scales

We wish to establish a buoyancy-drag equation for the MRTI within the Rapid Acceleration model. The first step is to derive the simplified evolution equations for the various spectral densities of interest, without nonlinear and dissipative terms, which are the vertical kinetic energy $\mathcal {E}_{uu}^{(33)}(\boldsymbol k) = \hat {u}_3(\boldsymbol k) \hat {u}_3(-\boldsymbol k)$, the vertical magnetic energy $\mathcal {E}_{bb}^{(33)}(\boldsymbol k) = \hat {b}_3(\boldsymbol k) \hat {b}_3(-\boldsymbol k)$, the vertical mass flux $\mathcal {E}_{uc}(\boldsymbol k) = \hat {c}(\boldsymbol k) \hat {u}_3(-\boldsymbol k)$, the vertical magnetic flux $\mathcal {E}_{bc}(\boldsymbol k) = \hat {c}(\boldsymbol k) \hat {b}_3(-\boldsymbol k)$, the concentration variance $\mathcal {E}_{cc}(\boldsymbol k) = \hat {c}(\boldsymbol k) \hat {c}(-\boldsymbol k)$ and the cross-helicity $\mathcal {E}_{ub}(\boldsymbol k) = \hat {b}_i(\boldsymbol k) \hat {u}_i(-\boldsymbol k)$. One eventually gets

(6.10a)$$\begin{gather} \frac{\partial \mathcal{E}_{uu}^{(33)}}{\partial t} ={-} 4 \mathcal{A}g {\rm Re} [ \mathcal{E}_{uc} ] \sin^2 \theta - 2 k \cos \theta B_0 {\rm Im} [\mathcal{E}_{ub}], \end{gather}$$

(6.10b)$$\begin{gather}\frac{\partial \mathcal{E}_{uc}}{\partial t} ={-}2 \mathcal{A}g \mathcal{E}_{cc} \sin^2 \theta - \frac{1}{L} \mathcal{E}_{uu}^{(33)} - \mathrm{i} k \cos \theta B_0 \mathcal{E}_{bc},\end{gather}$$

(6.10c)$$\begin{gather}\frac{\partial \mathcal{E}_{cc}}{\partial t} ={-}\frac{2}{L} \mathcal{E}_{uc}, \end{gather}$$

(6.10d)$$\begin{gather}\frac{\partial \mathcal{E}_{bb}^{(33)}}{\partial t} ={+} 2 k \cos \theta B_0 {\rm Im} [\mathcal{E}_{ub} ], \end{gather}$$

(6.10e)$$\begin{gather}\frac{\partial \mathcal{E}_{ub}}{\partial t} ={-}2 \mathcal{A}g \mathcal{E}_{bc}^* \sin^2 \theta + \mathrm{i} k \cos \theta B_0 (\mathcal{E}_{uu}^{(33)} - \mathcal{E}_{bb}^{(33)} ), \end{gather}$$

(6.10f)$$\begin{gather}\frac{\partial \mathcal{E}_{bc}}{\partial t} ={-} \frac{1}{L} \mathcal{E}_{ub}^* - \mathrm{i} k \cos \theta B_0 \mathcal{E}_{uc}, \end{gather}$$

where $()^*$ reflects the complex conjugate. Here, we have assumed homogeneous statistics for the production terms in (6.10b) and (6.10c): the full equations are given in Appendix A.3. All these spectral densities are related to the spherically averaged spectra previously addressed, see (5.4) and (6.1).

In contrast to the hydrodynamic case, the presence of a mean magnetic field introduces an explicit dependence to the wavevector $\boldsymbol k$ and adds multiple interactions between the different fields. An important outcome of these equations is that for the particular angle $\theta ={\rm \pi} /2$, one recovers the RA hydrodynamic equations written for $\theta ={\rm \pi} /2$. However, for this particular angle for which perturbations grow the fastest (see § 3.3), there is no feedback of $B_0$, meaning that the subtle effects of the mean magnetic field are contained in nonlinear interactions. A way to take these effects into account, which differs from the RA approach, is to spherically integrate these equations so that all $\theta$-contributions are accounted for, yielding

(6.11a)$$\begin{gather} \frac{\partial}{\partial t} (E_{uu}^{(33)} + E_{bb}^{(33)}) ={-} 4 \mathcal{A}g \int \sin^2 \theta {\rm Re} [ \mathcal{E}_{uc} ] \,\mathrm{d}^2 S_k , \end{gather}$$

(6.11b)$$\begin{gather}\frac{\partial E_{uc}}{\partial t} ={-}2 \mathcal{A}g \int \sin^2 \theta \mathcal{E}_{cc} \mathrm{d}^2 S_k - \frac{E_{uu}^{(33)}}{L} - \mathrm{i} B_0 \int k \cos \theta \mathcal{E}_{bc} \mathrm{d}^2 S_k, \end{gather}$$

(6.11c)$$\begin{gather}\frac{\partial E_{cc}}{\partial t} ={-}\frac{2}{L} E_{uc}, \end{gather}$$

with $\mathrm {d}^2 S_k = k^2 \sin \theta \,\mathrm {d}\theta \,\mathrm {d}\phi$ the surface integration at constant $k$ in spherical coordinates. Note that summing the kinetic and magnetic energy spectra allows to not consider the redistribution of energy between these two reservoirs.

To close this system, two additional assumptions must be made: (i) we assume a constant magnetic to kinetic energy ratio $R = E_{bb}^{(33)}/E_{uu}^{(33)}$ at large scales; (ii) at the leading order $\theta \to {\rm \pi}/2$, the last right-hand side term in (6.11b), which scales like $\cos \theta$, can be neglected.

Assumption (i) is justified by the fact that in figure 22(a), the ratio of magnetic to kinetic energy $\mathcal {K}_{bb}/\mathcal {K}_{uu}$ becomes roughly constant in the asymptotic state. This is further corroborated by the observations made previously in figure 18(a), namely that the kinetic and magnetic spectra have similar integral scales, and that the ratio $E_{bb}/E_{uu}$ is indeed roughly constant at large scales in figure 20(b). Note that $R = \overline {b^2}/\overline {u^2}$ can also be interpreted as the ratio between the Lorentz term $(\boldsymbol b \boldsymbol {\cdot } \boldsymbol \nabla ) \boldsymbol b \sim b^2/L$ and the advection term $(\boldsymbol u \boldsymbol {\cdot } \boldsymbol \nabla ) \boldsymbol u \sim u^2/L$. It follows that induction is never dominant during the mixing layer development. Nevertheless, the induced magnetic fields strongly modify the large-scale dynamics. Assumption (ii) can also be justified a posteriori since DNS results indicate that the correlation $\overline {b_3 c}$ remains orders of magnitude smaller than $\overline {u^2}$, $\overline {c^2}$ and $\overline {u_3c}$.

Figure 22.

(a) Time evolution of the ratio $R=\mathcal {K}_{bb}/\mathcal {K}_{uu}$ for various $B_0$ and diffusion coefficients; time is normalized by its final value $t_{end}$, which is larger for vertically elongated domains. (b) Evolution of $C_b C_d$ as a function of $B_0$ for various runs of table 1. Symbols correspond to different resolutions ($\times$, R1; $+$, R2; $\square$, R4) and colours to different diffusion coefficients.

Applying these two assumptions, the previous system of equations yields, at leading order in $\theta$:

(6.12a)$$\begin{gather} (1+R) \frac{\partial E_{uu}^{(33)}}{\partial t} ={-} 4 \mathcal{A}g E_{uc}, \end{gather}$$

(6.12b)$$\begin{gather}\frac{\partial E_{uc}}{\partial t} ={-}2 \mathcal{A}g E_{cc} - \frac{1}{L} E_{uu}^{(33)}, \end{gather}$$

(6.12c)$$\begin{gather}\frac{\partial E_{cc}}{\partial t} ={-}\frac{2}{L} E_{uc}. \end{gather}$$

This new system of equations describes the large-scale dynamics of the MRTI, and is very close to the system derived by Gréa (Reference Gréa2013) for the hydrodynamic case, except that here, it is spherically integrated. The present system (6.12) allows to relate the magnetic to kinetic energy ratio $R$ to the growth rate $\alpha$ and the buoyancy and drag coefficients through (6.8).

The following step is to determine which of $C_b$ and $C_d$ is modified by the presence of turbulent magnetic energy. To do so, we come back to the system (6.11) and expand the spectra in Legendre polynomials, as was done by Gréa (Reference Gréa2013), using in addition the constant magnetic to kinetic energy ratio assumption. After some algebra, one gets the MRTI buoyancy-drag equation

(6.13)\begin{equation} \ddot{L} ={-} \left( \frac{1}{2} + \frac{R}{2R+3} \right) \frac{\dot{L}^2}{L} + 4 \mathcal{A}g \sin^2 \gamma (1-\varTheta), \end{equation}

where the drag coefficient is now $C_d = 1/2 + R/(2R+3)$, increased by a factor $R/(2R+3) > 0$ with respect to the non-magnetic case, for which $R=0$ and (6.9) is recovered. Equation (6.13) is an important result since it proves analytically, under the assumption of constant magnetic to kinetic energy ratio, that the growth rate $\alpha$ of the mixing zone is indeed reduced in the presence of a mean magnetic field, due to an increase of the drag coefficient. However, this is not a predictive result since turbulent diffusion is not taken into account: hence, we have to propose a new closure between the drag coefficient $C_d$ and the buoyancy coefficient $C_b$ in the MRTI framework.

6.4.3. The damped growth rate of the MRTI

Before choosing a new closure for $C_d$ and $C_b$ that accounts for the effects of a vertical mean magnetic field in the asymptotic state, we investigate how the quantity $C_b C_d = C_\alpha$ varies for various $B_0$, $g$ and diffusion coefficients in figure 22(b). Using for the buoyancy coefficient $C_b = 4 \sin ^2 \gamma (1-\varTheta )$, $C_d$ is determined so that $\alpha _{bd}$ in (6.8) recovers the numerical value $\alpha = \dot {L}^2/(8 \mathcal {A}gL)$ in the asymptotic state with a standard least square fit.

The result is rather clear in figure 22(b): there is an increase of $C_b C_d = C_\alpha$ for larger values of $B_0$, and then saturation. Moreover, unlike the purely hydrodynamic case where asymptotically $\alpha = \alpha _{hyd}$ irrespective of the kinematic viscosity, $C_\alpha$ depends here on the diffusion coefficients $\nu = \kappa = \eta$ in the presence of a vertical mean magnetic field. This is expected, since we have observed in figure 22(a) that the magnetic to kinetic energy ratio, and hence the drag coefficient, depends on the turbulence intensity.

Furthermore, we also observe a reduction of the drag coefficient when the magnetic Prandtl number is decreased (green crosses): this is also expected since magnetic energy is more rapidly dissipated with larger $\eta$ so that the magnetic to kinetic energy ratio becomes smaller.

This analysis also shows that $C_\alpha$ does not depend on $\mathcal {A}g$: indeed, for simulations R1B01G20 and R1B02G20, where $g$ is increased from 10 to 20, the same $C_\alpha$ is obtained despite a slightly larger $\alpha$.

To propose a new closure for $C_b$ and $C_d$, it is quite natural to extend the methodology proposed by Poujade & Peybernes (Reference Poujade and Peybernes2010), but starting now from the system (6.12), where we keep our assumption that the magnetic to kinetic energy ratio is constant in the asymptotic state. Then, we assume that the large scales of the various spectra evolve self-similarly, which is shown for $E_{uu}$ in figure 21(b), so that one may write $E_{xy} = E_{xy}^0 k^s t^{n_{xy}}$, where $xy$ stands for either $uu$, $uc$ or $cc$, and with the amplitudes $E_{xy}^0$ and the exponents $n_{xy}$ independent of space and time. The infrared slope $s$ is the same for all spectra, which is also confirmed in figure 18(a). Finally, using (6.13), which states that the additional effects of the vertical mean magnetic field are only contained in $C_d$, one ends up with a new closure for the drag coefficient

(6.14)\begin{equation} C_d = \frac{2}{C_b} (1+R) + \frac{R}{2}.\end{equation}

Hence, using (6.8), the new prediction for the bubble growth rate in the magnetic Rayleigh–Taylor instability is

(6.15)\begin{equation} \alpha_{MRTI} = \frac{\sin^4 \gamma (1-\varTheta)^2}{(1+R) [1+ \sin^2 \gamma (1-\varTheta)]} = \frac{\alpha_{hyd}}{1+R}. \end{equation}

This expression illustrates the direct impact of the vertical mean magnetic field on the growth rate of the turbulent mixing zone: potential energy is not only converted into kinetic energy, but also into magnetic energy, which enhances drag by creating elongated structures. A numerical application of (6.15) shows that for a growth rate of order $\alpha =0.022$ in the hydrodynamic case, it could be reduced to $\simeq 0.0183$ or $\simeq 0.0142$ for a magnetic to kinetic energy ratio $0.2 \leq R \leq 0.55$, which represents a $35\,\%$ damping of the growth rate in the most turbulent cases considered here.

Similar to the anisotropy indicator $\sin ^2 \gamma$ and the mixing parameter $\varTheta$, the magnetic to kinetic energy ratio $R$ is flow dependent. It appears from figure 22(a) that the asymptotic value of $R$ obviously depends on parameters like $B_0$ and the diffusion coefficients. In the following, we define $R_{ss}$ as the self-similar value of $R$, obtained by averaging over the two last times of the simulations (over the range $t_{end}-2 \leq t \leq t_{end}$): these values are reported in table 1. It is clear in figure 23(a) that $R_{ss}$ scales linearly with the magnetic Froude number $F_M = B_0/V_{ref}$, defined in (4.5a,b), and then reaches a sort of plateau for the largest values of $F_M$. Although additional simulations at larger Reynolds numbers would be required to conclude this matter, it seems that the induced turbulent magnetic energy is bounded in the MRTI with an imposed vertical mean magnetic field. When the magnetic Prandtl number $P_M$ is decreased (green squares in figure 23b), $R_{ss}$ is consistently smaller than the $P_M=1$ case (green circle) since magnetic energy is more dissipated.

Figure 23.

(a) Self-similar value of the magnetic to kinetic energy ratio $R_{ss}$ as a function of the magnetic Froude number $F_M$: colours as in figure 8(b). (b) Ratio of the measured growth rate $\alpha$ on the hydrodynamic prediction $\alpha _{hyd}$ defined in (5.3) (blue circles) and on the MHD prediction $\alpha _{MRTI}$ defined in (6.15) (red squares). The thin dashed lines are visual guides.

The new prediction (6.15) for the growth rate in the MRTI is tested against the simulations in figure 23(b) for all cases of table 1. The growth rates are evaluated from simulations using $\alpha = \dot {L}^2/(8 \mathcal {A}g L)$ and are averaged over the two last times.

There is roughly a $\pm 5\,\%$ uncertainty when evaluating the growth rates, for example, coming from the dependence upon a delayed time $t_0$ if fits are used. In the present case, we have verified that the asymptotic values for $\alpha$ were consistent between an evaluation with $L$ and $\dot {L}$, and for instance by fitting $\sqrt {L/(\mathcal {A}g)}$ with a first-order polynomial.

It is clear that the ratio $\alpha /\alpha _{hyd}$ of measured growth rate on the hydrodynamic prediction (5.3) is well below unity, meaning that the prediction always overestimates the growth rate when the intensity of the mean magnetic field increases. The horizontal guide line shows that $\alpha \simeq 0.67 \alpha _{hyd}$ for large $B_0$. In contrast, the ratio $\alpha /\alpha _{MRTI}$ of measured growth rate on the MHD prediction (6.15) is close to unity, with excellent agreement for the largest Alfvén velocities. For smaller ones, namely $B_0 \leq 0.1$, the magnetic energy spectra are not developed enough so that the assumption of a constant ratio $R$ is questionable; still, it provides better results than the non-magnetic prediction.