%20-%20https://sketch.com%20--%3e%3ctitle%3e3ECBBF50-FD08-412A-943D-2CFC5E025A16%3c/title%3e%3cdesc%3eCreated%20with%20sketchtool.%3c/desc%3e%3cg%20id='Stylesheet'%20stroke='none'%20stroke-width='1'%20fill='none'%20fill-rule='evenodd'%3e%3cg%20transform='translate(-926.000000,%20-668.000000)'%20id='BUTTONS'%3e%3cg%20transform='translate(100.000000,%20597.000000)'%3e%3cg%20id='BTN-/-RIGHTS'%20transform='translate(816.000000,%2065.000000)'%3e%3cg%20id='ICON-/-Rights'%20transform='translate(10.000000,%206.000000)'%3e%3cpath%20d='M13.5752841,0%20L1.34943182,0%20C0.59375,0%200,0.586666667%200,1.33333333%20L0,22.6666667%20C0,23.4133333%200.620738636,24%201.34943182,24%20L17.6505682,24%20C18.40625,24%2019,23.3866667%2019,22.6666667%20L19,5.36%20L13.5752841,0%20Z%20M13.5752841,1.89333333%20L17.0838068,5.36%20L13.5752841,5.36%20L13.5752841,1.89333333%20Z%20M1.34943182,22.6666667%20L1.34943182,1.33333333%20L12.1988636,1.33333333%20L12.1988636,5.36%20C12.1988636,6.10666667%2012.7926136,6.69333333%2013.5482955,6.69333333%20L17.6235795,6.69333333%20L17.6505682,22.6666667%20L1.34943182,22.6666667%20Z'%20id='Shape'%20fill='%23595959'%20fill-rule='nonzero'%3e%3c/path%3e%3ccircle%20id='Oval'%20stroke='%23595959'%20stroke-width='1.5'%20fill='%23595959'%20cx='9.5'%20cy='13.5355339'%20r='5.5'%3e%3c/circle%3e%3cpolyline%20id='Path-2'%20stroke='%23FFFFFF'%20stroke-width='1.5'%20transform='translate(9.500000,%2013.035534)%20rotate(45.000000)%20translate(-9.500000,%20-13.035534)%20'%20points='8%2016.0355339%2011%2016.0355339%2011%2010.0355339'%3e%3c/polyline%3e%3c/g%3e%3c/g%3e%3c/g%3e%3c/g%3e%3c/g%3e%3c/svg%3e)
1. Introduction
The magnetic field plays a vital role in structuring the atmosphere of stars like the Sun and, in general, in the dynamics of astrophysical plasma. The interplay of non-propagating (NP) turbulent motions and propagating plasma waves is a fundamental process in astrophysical magnetohydrodynamics (MHDs) (see, e.g. Marino & Sorriso-Valvo Reference Marino and Sorriso-Valvo2023). In addition to the presence of Alfvén waves, acoustic waves and MHD Rossby waves in the solar atmosphere (see, e.g. Straus et al. Reference Straus, Fleck, Jefferies, Cauzzi, McIntosh, Reardon, Severino and Steffen2008; Wiśniewska et al. Reference Wiśniewska, Musielak, Staiger and Roth2016; Grant et al. Reference Grant, Jess, Zaqarashvili, Beck, Socas-Navarro, Aschwanden, Keys, Christian, Houston and Hewitt2018; Zaqarashvili et al. Reference Zaqarashvili, Lomineishvili, Leitner, Hanslmeier, Gömöry and Roth2021), several observations provided strong evidence for the presence also of internal gravity waves (IGW) (Straus et al. Reference Straus, Fleck, Jefferies, Cauzzi, McIntosh, Reardon, Severino and Steffen2008; Kneer & González Reference Kneer and González2011; Nagashima et al. Reference Nagashima, Löptien, Gizon, Birch, Cameron, Couvidat, Danilovic, Fleck and Stein2014; Vigeesh & Roth Reference Vigeesh and Roth2020). The IGW propagating through the lower solar atmosphere, where buoyancy plays the role of the dominant resorting force, may be generated by turbulent motions close to the visible solar surface (see Hague & Erdélyi Reference Hague and Erdélyi2016). Because the magnetic activity is ubiquitous throughout the solar atmosphere, so it is expected that the behaviour of IGW is to be affected (see Vigeesh, Jackiewicz & Steiner Reference Vigeesh, Jackiewicz and Steiner2017). In this work, we study the dynamics of Alfvén waves, magnetogravity waves and NP modes occurring in stably stratified MHD turbulence subject to a uniform mean magnetic field for a electrically conductive, Boussinesq fluid. In the frame of the Boussinesq approximation (see Spiegel & Veronis Reference Spiegel and Veronis1960), higher-frequency acoustic waves are filtered out and, hence, the Boussinesq MHD equations (Davidson Reference Davidson2013) can be used to study the turbulent motions in stably stratified regions of stars and gas giants below their surfaces (see, e.g. Skoutnev Reference Skoutnev2023).
The introduction of a uniform background magnetic field 
$\boldsymbol {B}_0$ or the buoyancy force (under the Boussinesq approximation) reduces the number of inviscid/ideal invariants in three-dimensional (3-D) incompressible MHD. In the first case, the magnetic helicity is no longer conserved while in the second it is the cross-correlation between the velocity and the magnetic field. For the Boussinesq MHD equations, the Ertel potential vorticity (PV), denoted here by 
$\tilde {\varPi }_\kappa,$ which is the fundamental correlation between vorticity and stratification (e.g. Pedlosky Reference Pedlosky2013), is not a Lagrangian invariant in MHD because it removes the baroclinic torque in the extended vorticity equation, but not the counterpart of the Lorentz force. In contrast, the so-called ‘magnetic induction potential scalar’ (MIPS) (i.e. the scalar product of the magnetic field vector and the buoyancy scalar gradient, denoted here by 
$\tilde {\varPi })$ is a Lagrangian invariant for a non-diffusive, electrically conducting, Boussinesq fluid as shown by Salhi et al. (Reference Salhi, Lehner, Godeferd and Cambon2012) (see also Salhi et al. Reference Salhi, Baklouti, Godeferd, Lehner and Cambon2017; Salhi & Cambon Reference Salhi and Cambon2023).
Consider the two systems: stratified MHD turbulence subject to a uniform mean magnetic field (
$\boldsymbol {B}_0=B_0\hat {\boldsymbol {z}})$ and (non-magnetized) stratified turbulence with system rotation (
$\boldsymbol {\varOmega } =(f/2) \hat {\boldsymbol {z}}).$ Here, 
$f>0$ is the Coriolis parameter and 
$B_0>0$ is the Alfvén velocity scaled from the mean magnetic field. For both systems, the mean density gradient aligns with the unit vertical vector 
$\hat{\boldsymbol{z}},$ with constant strength 
$N$ (i.e. the Brunt–Väisälä frequency). Accordingly, 
$\tilde {\varPi }_\kappa$ and 
$\tilde {\varPi }$ take the form (see § 2)
$$\begin{gather} \tilde{\varPi}_\kappa=fN+N\omega_z +f\partial_z\vartheta+\boldsymbol{\omega}\boldsymbol{\cdot} \boldsymbol{\nabla} \vartheta, \end{gather}$$
$$\begin{gather}\tilde{\varPi}=B_0N+Nb_z +B_0\partial_z\vartheta+\boldsymbol{b}\boldsymbol{\cdot} \boldsymbol{\nabla}\vartheta. \end{gather}$$
Here, 
$\omega _z=\boldsymbol {\omega }\boldsymbol {\cdot }\hat{\boldsymbol{z}}$ and 
$b_z=\boldsymbol {b}\boldsymbol {\cdot }\hat{\boldsymbol {z}}$ are the vertical components of the vorticity vector, 
$\boldsymbol {\omega }=\boldsymbol {\nabla }\times \boldsymbol {u},$ and the fluctuating magnetic field, respectively, and 
$\vartheta$ denotes the fluctuating buoyancy scalar. Note that the constant parts 
$fN$ and 
$B_0N$ do not participate in the dynamics and therefore can be ignored. The question that may arise from (1.1) is whether the similarity between 
$\tilde {\varPi }_\kappa$ and 
$\tilde {\varPi }$ is or not limited to a formal analogy.
It appears useful to briefly recall some previous studies characterizing the role played by 
$\varPi _\kappa$ in the dynamics of geophysical flows and stratified turbulence. The evolution of weak, homogeneous turbulence with strong rotation and stratification was recently revisited by Scott & Cambon (Reference Scott and Cambon2024). In the linear limit, the flow consists of oscillatory spectral modes which represent inertia-gravity waves with frequency 
$\omega _{ig}$,
\begin{equation} \omega^2_{ig}=f^2k^{{-}2}k_\parallel^2+N^2k^{{-}2}k_\perp^2, \end{equation}
and time-independent ones that express a NP component of the flow. Here 
$k_{\perp }=\Vert \boldsymbol {k}\times \hat {\boldsymbol {z}}\Vert$ and 
$k_\parallel =\boldsymbol {k}\boldsymbol {\cdot }\hat {\boldsymbol {z}}$ are, respectively, the horizontal and vertical components of the wavevector 
$\boldsymbol {k}.$ The evolution under weak nonlinearity, which is the purpose of the wave turbulence theory, was obtained without stratification (e.g. Galtier Reference Galtier2003; Bellet et al. Reference Bellet, Godeferd, Scott and Cambon2006), whereas the presence of the NP mode renders more problematic the case with coupled stratification. whereas the presence of the NP mode makes it more problematic the case with coupled stratification. The NP mode corresponds to PV and its dynamics, when it is decoupled from waves. It is very close to what one has in the quasigeostrophic (QG) theory (see Pedlosky Reference Pedlosky2013), which is one of the cornerstones in the study of atmospheric and oceanic flows since its development by Charney (Reference Charney1948, Reference Charney1971). The new results with direct numerical simulations (DNS) without forcing for the NP component are consistent with the prediction of Charney (Reference Charney1971), both in terms of power laws and the lack of a cascade (see also Leith Reference Leith1980; Embid & Majda Reference Embid and Majda1998; Kurien, Wingate & Taylor Reference Kurien, Wingate and Taylor2008). In the case of the evolution of both the NP component and the inertia-gravity waves one, via wave turbulence theory, the emphasis is put on angle-dependent spectra for a wide range of the ratio 
$f/N,$ as the most detailed description of anisotropy.
On the other hand, an inverse cascade is suggested by previous quasilinear studies as follows. The inviscid statistical tendency at low Rossby 
$R_o=U/(fL)$ and Froude 
$F_r=U/(NL)$ has been explored theoretically and numerically by, e.g. Bartello (Reference Bartello1995), Kurien, Smith & Wingate (Reference Kurien, Smith and Wingate2006), Kurien et al. (Reference Kurien, Wingate and Taylor2008), Herbert, Pouquet & Marino (Reference Herbert, Pouquet and Marino2014) and Herbert et al. (Reference Herbert, Marino, Rosenberg and Pouquet2016). Here, 
$U$ and 
$L$ are characteristic velocity and length scales, respectively. Bartello (Reference Bartello1995) employed statistical mechanics to consider the effects of combined conservation of total (
${\rm kinetic}+{\rm potential}$) energy and potential enstrophy (i.e. the 
$L_2$ norm of PV). Indeed, in the limit 
$F_r\rightarrow 0$ and/or 
$R_o\rightarrow 0,$ the potential enstrophy can be approximated by its quadratic part, and hence, the two inviscid invariants can be expressed in terms of the abovementioned normal modes. In the study by Herbert et al. (Reference Herbert, Pouquet and Marino2014), it was shown that restricting the equilibrium probability distribution to the slow manifold produces an inverse cascade. In contrast, taking into account the whole phase space, including the waves, results in a direct cascade (see also Lucarini et al. Reference Lucarini, Blender, Herbert, Ragone, Pascale and Wouters2014).
Moderate and strong nonlinearity is taken into account in several pseudospectral DNS studies in triperiodic boxes: the presence an inverse cascade of energy is evidenced, but most of these computations are forced and are eventually affected by finite-box effects, sometime also using elongated boxes. It is shown that turbulent QG flow produces an inverse energy cascade (see, e.g. Vallgren & Lindborg Reference Vallgren and Lindborg2010) resulting in a total (
${\rm kinetic}+{\rm potential}$) energy spectrum scaling as 
$E(k)\propto k^{-5/3}$ for the low wavenumbers, and a forward cascade of potential enstrophy, analogous to the enstrophy cascade in two-dimensional turbulence, with 
$E(k)\propto k^{-3}$ at high wavenumbers (see also Nastrom & Gage Reference Nastrom and Gage1983; Lilly Reference Lilly1989; Maltrud & Vallis Reference Maltrud and Vallis1991). Note that the QG enstrophy cascade from the large-scale circulation is one source of PV at intermediate and small scales in the atmosphere and ocean (see, e.g. Waite & Richardson Reference Waite and Richardson2023). Recently, van Kan & Alexakis (Reference van Kan and Alexakis2022), who studied forced rapidly rotating and stably stratified turbulence performing DNS with an elongated domain, showed that there is an inverse cascade even for weak stratification (large 
$Fr$) but for sufficiently large 
$\lambda =(hRo)^{-1}$ where 
$h$ is the ratio of the domain height to the energy injection scale. Other previous DNS studies of rapidly rotating stratified turbulence report the observation of a split energy cascade where the ratio 
$N/f$ has been identified as a control parameter: rotating and stratified flows at moderate values 
$N/f$ develop inverse cascade (see, e.g. Smith & Waleffe Reference Smith and Waleffe2002; Marino et al. Reference Marino, Rosenberg, Herbert and Pouquet2015; Herbert et al. Reference Herbert, Marino, Rosenberg and Pouquet2016).
In the stratified MHD flows with mean magnetic fields 
$(B_0\ne 0),$ the normal mode analysis in the inviscid linear limit indicates that, in the case of small 
$F_r$ and small Alfvén–Mach number 
$M=U/B_0,$ there are three decoupling normal modes: a NP mode, fast Alfvén and fast magnetogravity waves whose dispersion relations are
\begin{equation} \omega_a^2=B_0^2k_\parallel^2,\quad \omega_{ag}^2=\omega_a^2+\omega_g^2= B_0^2k^2_\parallel+N^2k^{{-}2}k_\perp^2, \end{equation}
respectively. In counterpart, in the case of small 
$F_r$ and finite 
$M,$ there are fast gravity waves with frequency 
$\omega _g$ and slow modes: a NP mode and slow Alfvén waves with frequency 
$\omega _a.$
For the assessment of the effect of stratification and of a background magnetic field, separately, we refer to the cases of purely stratified turbulence (PST) and of purely MHD turbulence with a mean magnetic field. Here we briefly report some results from the literature relative to these two configurations.
As it was shown in previous studies, the buoyancy Reynolds number 
$R{eb}=U_h^3/(\nu L_h N^2)$ can be considered a primary control parameter in characterizing PST, as far as is concerned energy transfer, dissipation, mixing and dispersion properties (Smyth & Moum Reference Smyth and Moum2000; Billant & Chomaz Reference Billant and Chomaz2001; Brethouwer et al. Reference Brethouwer, Billant, Lindborg and Chomaz2007; Bartello & Tobias Reference Bartello and Tobias2013; Waite Reference Waite2013; Maffioli & Davidson Reference Maffioli and Davidson2016). Here, 
$U_h$ and 
$L_h$ denote characteristic horizontal velocity and length 
$L_h$ scales, respectively, and 
$\nu$ being the kinematic viscosity. The potential relevance of 
$R_{eb}$ is based on the shearing dynamics observed to occur when strong stratification leads to the tendency for vertical decoupling of horizontal motions (Riley & de Bruyn Kops Reference Riley and de Bruyn Kops2003; de Bruyn Kops & Riley Reference de Bruyn Kops and Riley2019). In the lower-
$R_{eb}$ regime (called the viscosity-affected stratified flow regime (Watanabe, Zheng & Nagata Reference Watanabe, Zheng and Nagata2022)), the flow is characterized by thin layers of horizontal flow, glued together by the viscosity acting along the vertical direction Waite (Reference Waite2013). In previous studies, using and/or comparing anisotropic multimodal eddy damped quasinormal Markovian (EDQNM) and DNS, the layering was quantified by angle-dependent spectra, showing the concentration of energy towards vertical wavevectors (Godeferd & Cambon Reference Godeferd and Cambon1994); it was shown in Liechtenstein, Godeferd & Cambon (Reference Liechtenstein, Godeferd and Cambon2005) that this anisotropic nonlinear effect affected almost exclusively the NP mode. In that case, the inviscid potential enstrophy can be approximated by its quadratic part which dominates the higher-order terms allowing for the use of equilibrium statistical mechanics which indicates the lack of an inverse cascade of energy (see Waite & Bartello Reference Waite and Bartello2004; Herbert et al. Reference Herbert, Pouquet and Marino2014). In the high-
$R_{eb}$ regime (also called strongly stratified turbulence regime), stratified turbulence is characterized by layerwise ‘pancake’ flows with associated Kelvin–Helmoltz instabilities, consistent with the breakdown of the quadratic approximation of PV (see Waite Reference Waite2013). Note that the scaling analysis of Billant & Chomaz (Reference Billant and Chomaz2001) in strongly stratified turbulence regime gives 
$F_{rv}=U_h/(NL_v)\sim 1,$ where 
$L_h$ is the horizontal length scale. We also note that the numerical simulations by Kimura & Herring (Reference Kimura and Herring2012) for sufficiently strong stratification (
$R_{eb}>1)$ indicate that the wave spectrum is a steeper than 
$k_\perp ^{-2}$, while that for the NP mode is consistent with 
$k_\perp ^{-3}$. Our numerical simulations for the purely stratified case are characterized by 
$R_{eb}<1,$ 
$F_{rh}<0.1$ and 
$F_{rv}< 0.1$ meaning that the flows we considered in those cases do correspond to a stratified weak wave MHD turbulence regime that is rather affected by viscosity (see § 4).
The MHD turbulence in the presence of a mean magnetic field has been the subject of many studies (Iroshnikov Reference Iroshnikov1963; Kraichnan Reference Kraichnan1965; Shebalin, Matthaeus & Montgomery Reference Shebalin, Matthaeus and Montgomery1983; Galtier et al. Reference Galtier, Nazarenko, Newell and Pouquet2000; Nazarenko Reference Nazarenko2011; Alexakis & Biferale Reference Alexakis and Biferale2018). Indeed, the presence of a mean magnetic field (here, 
$\boldsymbol {B}_0=B_0\hat {\boldsymbol {z}})$ in (non-stratified) MHD turbulence supports the development of a high degree of anisotropy with most of the energy cascading perpendicular to 
$\boldsymbol {B}_0$. The turbulent fluctuations are elongated along the mean magnetic field and the anisotropy is scale-dependent (Müller, Biskamp & Grappin Reference Müller, Biskamp and Grappin2003). Thus, two time scales can be introduced: the nonlinear eddy turnover time the nonlinear eddy turnover time 
$\tau _{nl}^{-1}=k_\perp b_\lambda$ where 
$b_\lambda$ the root mean square (r.m.s.) magnetic fluctuations at the scale 
$\lambda =k_\perp ^{-1},$ and the Alfvén time 
$\tau _a^{-1}=k_\parallel B_0.$ The latter time scale can be interpreted as the interaction time between two counter-propagating Alfvén wave packets, in the 
$\pm \boldsymbol {B}_0$ directions. Depending on the ratio 
$\tau _{nl}/\tau _a$ different turbulent regimes can be distinguished (Zhou Reference Zhou2010). In the case where 
$\tau _a^{-1}=k_\parallel B_0\gg \tau _{nl}^{-1}=k_\perp b_\lambda,$ or equivalently, when the linear terms dominate, there is a weak turbulence regime. Otherwise, turbulence is called strong. A direct evidence transition from weak to strong MHD turbulence has been presented by Meyrand, Galtier & Kiyani (Reference Meyrand, Galtier and Kiyani2016) that showed how the change of regime is characterized, among other things, by a variation of the slope of the energy spectrum going from approximately 
$-2$ to 
$-3/2$ and by an increase of the ratio 
$\tau _a/\tau _{nl}.$ We note that Galtier et al. (Reference Galtier, Nazarenko, Newell and Pouquet2000) had analysed the weak turbulence limit of MHD turbulence for large 
$B_0$ and showed that the anisotropic spectrum scales as 
$E(k_\perp )\sim k_\perp ^{-2}.$ On the other hand, in the strong MHD turbulence regime, phenomenological models taking into account anisotropy do predict that the anisotropic energy spectrum would scale as 
$k_\perp ^{-5/3}$ (i.e. Kolmogrov's exponent (Goldreich & Sridhar Reference Goldreich and Sridhar1995) or as 
$k_\perp ^{-3/2}$ (Boldyrev Reference Boldyrev2005)).
The outline of the paper is structured as follows. In § 2, we present the mathematical formulation used in our study, including the Boussinesq-MHD equations as well as the equations for conserved quantities of an inviscid system. The normal mode decomposition in the inviscid linear limit is provided in § 3. The main DNS results are then presented and discussed in § 4. In § 5, we report our concluding remarks.
2. Mathematical formulation
2.1. The Boussinesq-MHD equations
The dynamics of a stably stratified fluid subjected to a mean uniform magnetic field, 
$\boldsymbol {B}_0$ can be described, under the Boussinesq approximation, by the momentum equation for the velocity, 
$\boldsymbol {u},$ the induction equation for the magnetic field, 
$\boldsymbol {b},$ and an equation for the buoyancy scalar, 
$\vartheta$ (e.g. Davidson Reference Davidson2013),
$$\begin{gather} \partial_t{\boldsymbol{u}}+\left (\boldsymbol{u}\boldsymbol{\cdot}\boldsymbol{\nabla}\right)\boldsymbol{u}={-}\boldsymbol{\nabla}p+\left (\left ({\boldsymbol{B}}_0+\boldsymbol{b}\right)\boldsymbol{\cdot} \boldsymbol{\nabla}\right)\boldsymbol{b}+N\vartheta \hat{\boldsymbol{z}} +\nu \nabla^2 {\boldsymbol{u}}, \end{gather}$$
$$\begin{gather}\partial_t{\boldsymbol{b}}+\left (\boldsymbol{u}\boldsymbol{\cdot}\boldsymbol{\nabla}\right)\boldsymbol{b}= \left ({\boldsymbol{B}}_0\boldsymbol{\cdot}\boldsymbol{\nabla}\right){\boldsymbol{u}} + \left ({\boldsymbol{b}}\boldsymbol{\cdot}\boldsymbol{\nabla}\right){\boldsymbol{u}}+\eta\nabla^2{\boldsymbol{b}}, \end{gather}$$
$$\begin{gather}\partial_t{\vartheta}+\left (\boldsymbol{u}\boldsymbol{\cdot}\boldsymbol{\nabla}\right)\vartheta={-}N\hat{\boldsymbol{z}}\boldsymbol{\cdot}\boldsymbol{u}+\kappa\nabla^2\vartheta, \end{gather}$$together with the incompressibility condition and the divergence-free condition for the magnetic field,
\begin{equation} \boldsymbol{\nabla}\boldsymbol{\cdot} {\boldsymbol{u}}=0,\quad \boldsymbol{\nabla}\boldsymbol{\cdot} {\boldsymbol{b}}=0. \end{equation}
Here, the magnetic field is scaled using Alfvén velocity units, i.e. it is divided by 
$\sqrt {\varrho _0\mu _0}$ where 
$\varrho _0$ and 
$\mu _0$ are a constant reference density and the magnetic permeability of the fluid, and the fluctuation of the buoyancy scalar 
$\vartheta =-(g/(N\varrho _0))\rho$ has the dimension of a velocity where 
$\rho$ is the density fluctuation and 
$N,$ such that
\begin{equation} N^2={-}\frac{g}{\rho_0}\frac{{\rm d}\varrho}{{\rm d}z}\geqslant 0 \end{equation}denotes the Brunt–Väisälä frequency which is a positive constant since we assume a linear background stable stratification,
\begin{equation} {\varrho=\varrho_0\left [1-\left (N^2/g\right )z\right]},\quad \boldsymbol{g}={-}g\hat{\boldsymbol{z}}, \end{equation}
with 
$g$ being the acceleration of gravity, and 
$\hat {\boldsymbol {z}}$ denotes a vertical upward unit vector such that 
$(\hat {\boldsymbol {x}},\hat {\boldsymbol {y}},\hat {\boldsymbol {z}})$ is the canonical basis of 
$\mathbb {R}^3.$ In the present paper, we consider that the magnetic and thermal Prandtl numbers, 
$P_m=\nu /\mu$ and 
$P_r=\nu /\kappa$ are equal to unity.
In (2.1), 
$p$ denotes the total pressure fluctuations (including the magnetic pressure) divided by 
$\varrho _0,$ and 
$\nu,$ 
$\eta$ and 
$\kappa$ denote the kinematic viscosity and the magnetic and density diffusivities, respectively. The second and third terms on the right-hand side of the momentum equation (2.1a) represent the Lorentz and the buoyancy forces, respectively. For simplicity, in the present study, we consider that the unperturbed magnetic field 
$\boldsymbol {B}_0$ aligns with the 
$\hat {\boldsymbol {z}}$-direction,
\begin{equation} \boldsymbol{B}_0=B_0\hat{\boldsymbol{z}}, \end{equation}
where 
$B_0$ is a positive constant. Note that the pressure 
$p$ can be eliminated by computing the divergence of the momentum equation in (2.1) together with 
$\boldsymbol {\nabla }\boldsymbol {\cdot }\boldsymbol {u}=0$ and 
$\boldsymbol {\nabla }\boldsymbol {\cdot }\boldsymbol {b}=0,$ and then solving the resulting elliptic equation for the pressure 
$p$ to get
\begin{equation} p(\boldsymbol{x},t)={-}\nabla^{{-}2}\left (\boldsymbol{\nabla}\boldsymbol{\cdot}\left (\boldsymbol{u}\boldsymbol{\cdot}\boldsymbol{\nabla}\boldsymbol{u}- \boldsymbol{b}\boldsymbol{\cdot}\boldsymbol{\nabla}\boldsymbol{b}\right)-N\boldsymbol{\nabla}\vartheta\boldsymbol{\cdot}\hat{\boldsymbol{z}}\right). \end{equation}From (2.1) and (2.2a,b) we deduce that the sum of the kinetic, magnetic and potential energies,
\begin{equation} { E}=\tfrac{1}{2}\int\left ( \Vert \boldsymbol{u}(\boldsymbol{x})\Vert^2+\Vert \boldsymbol{b}(\boldsymbol{x})\Vert^2+\vartheta (\boldsymbol{x})^2\right){\rm d}^3\boldsymbol{x} \end{equation}is conserved for a non-diffusive Boussinesq fluid,
\begin{equation} \partial_t{E}={-}\int\left (\nu\Vert \boldsymbol{\omega}\Vert^2+\eta\Vert\boldsymbol{\,j}\Vert^2+\kappa\Vert \boldsymbol{\nabla} \vartheta\Vert^2\right){\rm d}^3\boldsymbol{x}, \end{equation}
where 
$\boldsymbol {\omega }=\boldsymbol {\nabla }\times \boldsymbol {u}$ is the vorticity vector and 
$\boldsymbol {j}$ denotes the normalized current density. Likewise, from (2.1) and (2.2a,b) we derive the equation for the so-called MIPS (Salhi et al. (Reference Salhi, Lehner, Godeferd and Cambon2012); see also Salhi et al. (Reference Salhi, Baklouti, Godeferd, Lehner and Cambon2017); Salhi & Cambon (Reference Salhi and Cambon2023))
$$\begin{gather} \tilde{\varPi}=\tilde{\boldsymbol{b}}\boldsymbol{\cdot} \boldsymbol{\nabla} \tilde{\vartheta}= NB_0+N\hat{\boldsymbol{z}}\boldsymbol{\cdot}\boldsymbol{b} +B_0\hat{\boldsymbol{z}}\boldsymbol{\cdot}\boldsymbol{\nabla}\vartheta+\boldsymbol{b}\boldsymbol{\cdot}\boldsymbol{\nabla}\vartheta, \end{gather}$$
$$\begin{gather}\partial_t\tilde{\varPi}+\left (\boldsymbol{u}\boldsymbol{\cdot}\boldsymbol{\nabla}\right)\tilde{\varPi}=\eta\left (\boldsymbol{\nabla} \tilde{\vartheta}\right)\boldsymbol{\cdot}\nabla^2\tilde{\boldsymbol{b}}+\kappa \tilde{\boldsymbol{b}}\boldsymbol{\cdot}\nabla^2\left (\boldsymbol{\nabla} \tilde{\vartheta}\right), \end{gather}$$
which makes it clear that 
$\tilde {\varPi }$ is conserved following the flow in the absence of magnetic and thermal diffusion. Note that the constant part 
$\varPi _0=NB_0$ does not participate in the dynamics and can, therefore, be neglected.
2.2. Non-dimensionalized version of the Boussinesq MHD equations
In order to non-dimensionalize the above Boussinesq MHD equations we introduce characteristic scales for the physical variables: 
$L$ is the length scale in both the horizontal and the vertical directions and 
$U$ is the velocity scale. As in the study of the case of non-magnetized rotating and stratified system by Embid & Majda (Reference Embid and Majda1998), a non-dimensionalized version of the Boussinesq MHD equations can be written in abstract form as
\begin{equation} \partial_{t'}\boldsymbol{Y}+{\boldsymbol{\mathcal{L}}}(\boldsymbol{Y})+{\boldsymbol{\mathcal{B}}}(\boldsymbol{Y},\boldsymbol{Y})-{\boldsymbol{\mathcal{D}}}(\boldsymbol{Y})= \boldsymbol{{0}}, \end{equation}
where 
$t'=tU/L$ denotes the dimensionless time unit, 
$\boldsymbol {Y}$ is given by
\begin{equation} \boldsymbol{Y}=\left(\begin{array}{c} \boldsymbol{u}'\\ \boldsymbol{b}'\\ \vartheta'\end{array}\right)=U^{{-}1} \left (\begin{array}{c} \boldsymbol{u}\\ \boldsymbol{b}\\ \vartheta\end{array}\right), \end{equation}the linear operator is given by
\begin{equation} {\boldsymbol{\mathcal{L}}}(\boldsymbol{Y})=\left (\begin{array}{c} -F_r^{{-}1}\vartheta'\hat{\boldsymbol{z}}-M^{{-}1}\partial_{x'_3}\boldsymbol{b}'- F_r^{{-}1}\boldsymbol{\nabla}'\left (\nabla'^{{-}2}\partial_{x'_3}\vartheta'\right)\\ -M^{{-}1}\partial_{x'_3}\boldsymbol{u}'\\ F_r^{{-}1}\boldsymbol{u}'\boldsymbol{\cdot}\hat{\boldsymbol{z}} \end{array}\right ), \end{equation}the quadratic operator is given by
\begin{equation} {\boldsymbol{\mathcal{B}}}\left(\boldsymbol{Y},\boldsymbol{Y}\right)=\left (\begin{array}{c} \boldsymbol{u}'\boldsymbol{\cdot} \boldsymbol{\nabla}'\boldsymbol{u}'-\boldsymbol{b}'\boldsymbol{\cdot}\boldsymbol{\nabla}'\boldsymbol{b}' -\boldsymbol{\nabla}'\left (\nabla'^{{-}2}\left (\boldsymbol{\nabla}'\boldsymbol{\cdot} \left (\boldsymbol{u}'\boldsymbol{\cdot} \boldsymbol{\nabla}'\boldsymbol{u}'-\boldsymbol{b}'\boldsymbol{\cdot} \boldsymbol{\nabla}'\boldsymbol{b}'\right)\right)\right)\\ \boldsymbol{u}'\boldsymbol{\cdot} \boldsymbol{\nabla}'\boldsymbol{b}'-\boldsymbol{b}'\boldsymbol{\cdot} \boldsymbol{\nabla}'\boldsymbol{u}'\\ \boldsymbol{u}'\boldsymbol{\cdot}\boldsymbol{\nabla}'\vartheta' \end{array}\right ) \end{equation}and the diffusion operator is given by
\begin{equation} {\boldsymbol{\mathcal{D}}}\left(\boldsymbol{Y}\right)=\left(\begin{array}{c} R_e^{{-}1}\nabla'^{2}\boldsymbol{u}'\\ R_e^{{-}1}P_{m}^{{-}1}\nabla'^{2}\boldsymbol{b}'\\ R_e^{{-}1}P_{r}^{{-}1}\nabla'^{2}\vartheta' \end{array}\right ). \end{equation}
Here, 
$\boldsymbol {\nabla }'=L^{-1}\boldsymbol {\nabla }$ and the non-dimensional numbers
\begin{equation} F_r=\frac{U}{NL},\quad M^{{-}1}=\frac{B_0}{U}, \quad R_e=\frac{U}{\nu L} \end{equation}are the Froude, the inverse Alfvén–Mach and Reynolds numbers, respectively. In a non-dimensional form, the expression of MIPS reduces to
\begin{equation} \varPi'=b'_3+F_rM^{{-}1}\partial_{x'_3}\vartheta'+F_r\boldsymbol{b}'\boldsymbol{\cdot}\boldsymbol{\nabla}'\vartheta'. \end{equation}First, we consider the asymptotic limit of the small Froude number and the finite Alfvén–Mach number (hereinafter, it is referred to as case (II)). In that case, we consider the scaling
\begin{equation} F_r=\varepsilon\ll 1,\quad M={O}(1), \quad R_e={O}(1),\quad P_m=P_r={O}(1) \quad \text{as}\ \varepsilon \rightarrow 0. \end{equation}Accordingly, (2.10) can be rewritten as
$$\begin{gather} \partial_{t'}\boldsymbol{Y}+\varepsilon^{{-}1}{\boldsymbol{\mathcal{L}}}_F(\boldsymbol{Y})+{\boldsymbol{\mathcal{L}}}_S(\boldsymbol{Y})+ {\boldsymbol{\mathcal{B}}}(\boldsymbol{Y},\boldsymbol{Y})-{\boldsymbol{\mathcal{D}}}(\boldsymbol{Y})=\boldsymbol{{0}}, \end{gather}$$
$$\begin{gather}\boldsymbol{Y}(t=0)=\boldsymbol{Y}(\boldsymbol{x}), \end{gather}$$
where the linear operator 
${\boldsymbol {\mathcal {L}}}=\varepsilon ^{-1}{\boldsymbol {\mathcal {L}}}_F+{\boldsymbol {\mathcal {L}}}_S$ in (2.10) decomposes into a fast component 
${\boldsymbol {\mathcal {L}}}_F$ and a slow component 
${\boldsymbol {\mathcal {L}}}_S$ such that
\begin{equation} {\boldsymbol{\mathcal{L}}}_F(\boldsymbol{Y})=\left (\begin{array}{c} -\vartheta'\hat{\boldsymbol{z}}-\boldsymbol{\nabla}'\left (\nabla'^{{-}2}\partial_{x'_3}\vartheta'\right)\\ 0\\ \boldsymbol{u}'\boldsymbol{\cdot}\hat{\boldsymbol{z}} \end{array}\right ),\quad {\boldsymbol{\mathcal{L}}}_S(\boldsymbol{Y})=\left (\begin{array}{c} -M^{{-}1}\partial_{x'_3}\boldsymbol{b}'\\ -M^{{-}1}\partial_{x'_3}\boldsymbol{u}'\\ 0 \end{array}\right), \end{equation}
while the expression of 
$\varPi '$ takes the form
\begin{equation} \varPi'=b'_3+\varepsilon \left (M^{{-}1}\partial_{x'_3}\vartheta'+\boldsymbol{b}'\boldsymbol{\cdot}\boldsymbol{\nabla}'\vartheta'\right). \end{equation} In the case of small 
$F_r$ and small 
$M$ (hereinafter, it is referred to as case (I)), we use the following scaling:
\begin{equation} F_r=\varepsilon,\quad M=C^{{-}1}\varepsilon,\quad R_e={O}(1),\quad P_m=P_r={O}(1) \quad \text{as}\ \varepsilon \rightarrow 0, \end{equation}
where 
$C$ is the proportionality coefficient. With this scaling the equation (2.10) can be rewritten as
$$\begin{gather} \partial_{t'}\boldsymbol{Y}+\varepsilon^{{-}1}{\boldsymbol{\mathcal{L}}}_F(\boldsymbol{Y})+{\boldsymbol{\mathcal{B}}}(\boldsymbol{Y},\boldsymbol{Y})- {\boldsymbol{\mathcal{D}}}(\boldsymbol{Y})=\boldsymbol{{0}}, \end{gather}$$
$$\begin{gather}\boldsymbol{Y}(t=0)=\boldsymbol{Y}(\boldsymbol{x}), \end{gather}$$
where the linear operator 
$\boldsymbol {\mathcal {L}}=\varepsilon ^{-1}{\boldsymbol {\mathcal {L}}}_F$ has only a fast component 
${\boldsymbol {\mathcal {L}}}_F$ (while 
${\boldsymbol {\mathcal {L}}}_S$ is zero),
\begin{equation} {\boldsymbol{\mathcal{L}}}_F(\boldsymbol{Y})=\left (\begin{array}{c} -\vartheta'\hat{\boldsymbol{z}}-C\partial_{x'_3}\boldsymbol{b}'- \boldsymbol{\nabla}'\left (\nabla'^{{-}2}\partial_{x'_3}\vartheta'\right)\\ -C\partial_{x'_3}\boldsymbol{u}'\\ \boldsymbol{u}'\boldsymbol{\cdot}\hat{\boldsymbol{z}} \end{array}\right ), \end{equation}
while the expression of 
$\varPi '$ is of the form
\begin{equation} \varPi'=b'_3+C\partial_{x'_3}\vartheta'+\varepsilon \boldsymbol{b}'\boldsymbol{\cdot}\boldsymbol{\nabla}'\vartheta'. \end{equation}
More importantly, the non-dimensional form of the total energy equation (i.e. (2.8)) is free of the singular parameter 
$\varepsilon$ in (2.19a,b) or in (2.23), and therefore it is valid uniformly in 
$\varepsilon.$
It clearly appears that in both the asymptotic limit of cases (I) and (II), the scaling results in the introduction of the singular term 
$\varepsilon ^{-1}{\boldsymbol {\mathcal {L}}}_F,$ and in the appearance of separated fast and slow scales of motion. This can be shown clearly in the linear analysis of the Boussinesq MHD equations without diffusion. Accordingly, we consider the linear equation associated with either (2.19a,b) or (2.23),
\begin{equation} \partial_{t'}\boldsymbol{Y}+\varepsilon^{{-}1}{\mathcal{L}}_F(\boldsymbol{Y})=\boldsymbol{{0}}. \end{equation}
We assume periodic boundary conditions in the spatial domain and utilize the Fourier eigenfunctions of the operator 
${\boldsymbol {\mathcal {L}}}_F$ (see the next section).
For case (II), the analysis of (2.25), which will be presented in the next section, shows that there are NP modes with frequency 
$\omega _{np} (\boldsymbol {k})=0$ and fast (gravity) waves with frequency 
$\pm \omega _g(\boldsymbol {k})$,
\begin{equation} \omega_g(\boldsymbol{k})=\frac{Nk_\perp}{k}=\varepsilon^{{-}1}\frac{U}{L}\frac{k_\perp}{k} \end{equation}
where 
$\boldsymbol {k}$ is the wavevector, 
$k=\Vert \boldsymbol {k}\Vert$ and 
$k_\perp =\Vert \boldsymbol {k}\times \hat{\boldsymbol{z}}\Vert$ denotes the horizontal wavenumber. Therefore, we may conclude that, in case (II), the equations admit fast (gravity) waves moving on (dimensionless) time scales of order 
${O}(\varepsilon ^{-1})$ and NP modes moving on (dimensionless) time scales of order 
${O}(1).$
For case (I) the analysis of (2.25), which will be presented in the next section, shows that there are NP modes with frequency 
$\omega _{np} (\boldsymbol {k})=0$ and fast waves composed of Alfvén waves of frequency 
$\pm \omega _a$ and magnetogravity waves of frequency 
$\pm \omega _{ag}$,
\begin{equation} \omega_a=B_0\vert k_\parallel\vert=\varepsilon^{{-}1}CU\vert k_\parallel\vert,\quad \omega_{ag}=\sqrt{\omega_a^2+\omega_g^2}= \varepsilon^{{-}1}\sqrt{C^2U^2k_\parallel^2+\frac{U^2}{L^2}\frac{k_\perp^2}{k^2}} \end{equation}
where 
$k_\parallel \equiv k_3=\boldsymbol {k}\boldsymbol {\cdot }\hat{\boldsymbol {z}}$ is the vertical wavenumber. In case (I), the equations admit fast (Alfvén and magnetogravity) waves moving on (dimensionless) time scales of order 
${O}(\varepsilon ^{-1})$ and NP modes moving on (dimensionless) time scales of order 
${O}(1)$.
3. Normal mode decomposition
We consider the linear equation associated with either (2.19a,b) or (2.23),
\begin{equation} \partial_t\boldsymbol{Y}+\varepsilon^{{-}1}{\mathcal{L}}_F(\boldsymbol{Y})=\boldsymbol{{0}}, \end{equation}and introduce the Fourier decomposition
\begin{equation} \boldsymbol{Y}(\boldsymbol{x},t)= \sum_{\boldsymbol{k} \in \hat{\mathcal{T}}}U\hat{\boldsymbol{Y}}(\boldsymbol{k},t)\text{e}^{\mathrm{i} \boldsymbol{x}\ \boldsymbol{\cdot}\ \boldsymbol{k}} \end{equation}
where 
$\hat {\mathcal {T}}=\mathbb {Z}^3$ is the set of wavevectors and 
$\mathrm {i}^2=-1.$ In Fourier space, the incompressibility constraint and the divergence-free condition for the magnetic field read 
$\boldsymbol {k}\boldsymbol {\cdot }\hat {\boldsymbol {u}}=0$ and 
$\boldsymbol {k}\boldsymbol {\cdot }\hat {\boldsymbol {b}}=0,$ respectively, which signifies that both 
$\hat {\boldsymbol {u}}$ and 
$\hat {\boldsymbol {b}}$ belong to the plane perpendicular to the wavevector 
$\boldsymbol {k}$.
In physical space, we even start from eight components, including the pressure, that is to say: three for the velocity vector; three for the magnetic fluid; one for the buoyancy scalar; and one for the pressure. The effect of the divergence-free condition is to reduce to two independent components both the velocity field and the magnetic field in the 3-D Fourier space, in the plane normal to 
$\boldsymbol {k}$. The pressure is removed since the Poisson equation (see (2.6)) is equivalent to an algebraic relation along 
$\boldsymbol {k}.$ Accordingly, in the 3-D Fourier space, only five components remain: two for the velocity; two for the magnetic field; and one for the buoyancy scalar.
3.1. Case of small 
$F_r$ and small $M$ (case (I))
First, we perform the Fourier analysis of the fast operator 
${\boldsymbol {\mathcal {L}}}_F$ in the case of small 
$F_r$ and small 
$M$ (case (I)). Given (2.23) the substitution of (3.2) into (3.1) leads to the linear differential system
\begin{equation} d_t{\hat{\boldsymbol{Y}}}={\hat{\boldsymbol{\mathcal{L}}}}\boldsymbol{\cdot}\hat{\boldsymbol{Y}} \end{equation}
where the non-zero elements of the matrix 
$\hat {\boldsymbol {\mathcal {L}}}(\boldsymbol {k})$ are
$$\begin{gather} \hat{\mathcal{L}}_{14}=\hat{\mathcal{L}}_{25}=\hat{\mathcal{L}}_{36}=\hat{\mathcal{L}}_{41}=\hat{\mathcal{L}}_{52}= \hat{\mathcal{L}}_{63}=\mathrm{i} k_3B_0, \end{gather}$$
$$\begin{gather}\hat{\mathcal{L}}_{17}=\frac{k_1k_3}{k^2} N,\quad \hat{\mathcal{L}}_{27}=\frac{k_2k_3}{k^2}N,\quad \hat{\mathcal{L}}_{37}={-}\frac{k_\perp^2}{k^2}N,\quad \hat{\mathcal{L}}_{73}=N. \end{gather}$$
The spectrum of the matrix 
$\hat {\boldsymbol {\mathcal {L}}}(\boldsymbol {k})$ is found as
\begin{equation} \text{Sp}\,\hat{\boldsymbol{\mathcal{L}}}(\boldsymbol{k})=\{\lambda_1={-}\mathrm{i}\omega_a, \lambda_2=\mathrm{i}\omega_a, \lambda_3=0, \lambda_4={-}\mathrm{i}\omega_{ag}, \lambda_5=\mathrm{i}\omega_{ag}\}, \end{equation}
where the eigenvalues 
$\lambda _1$ and 
$\lambda _2$ are of multiplicity 2, and
\begin{equation} \omega_a=\vert k_3\vert B_0,\quad \omega_g=k^{{-}1}{k_\perp} N,\quad \omega_{ag}=\sqrt{\omega_a^2+\omega_g^2}=\sqrt{k_3^2B_0^2+(k_\perp{/}k)^2N^2} \end{equation}denote the frequencies of the Alfvén, gravity and magnetogravity waves, respectively, as indicated at the end of the previous section.
The associated eigenvectors, say 
$\boldsymbol {Z}_\ell$ (
$\ell =1,2,\ldots,5),$ such that 
$\hat {\boldsymbol {\mathcal {L}}}\boldsymbol {\cdot }{\boldsymbol {Z}}_\ell =\lambda _\ell \boldsymbol {Z}_\ell,$ are, respectively, the columns of the following rectangular matrix:
\begin{equation} \boldsymbol{M}(\boldsymbol{k})= \frac{1}{\sqrt{2}kk_\perp} \left (\begin{array}{ccccc} {kk_2} & {kk_2} & 0 & k_1k_3 & -k_1k_3\\ {-kk_1} & {-kk_1} & 0 & k_2k_3 & -k_2k_3\\ 0 & 0 & 0 & -{k_\perp^2} & {k_\perp^2}\\ {-s kk_2} & {s kk_2} & \sqrt{2}k_1k_3{\omega_{g}^+} & -k_1k_3{s\omega_{a}^+} & -k_1k_3{s\omega_{a}^+}\\ {s kk_1} & {-s kk_1} & \sqrt{2}k_2k_3{\omega_{g}^+} & -k_2k_3{s\omega_{a}^+} & -k_2k_3{s\omega_{a}^+}\\ 0 & 0 & -\sqrt{2}{k_\perp^2}\omega_g^+ & {k_\perp^2}s\omega_a^+ & {k_\perp^2}s\omega_a^+\\ 0 & 0 & -\mathrm{i}\sqrt{2}{kk_\perp}s\omega_a^+ & -\mathrm{i}{kk_\perp}\omega_g^+ & -\mathrm{i}{kk_\perp}\omega_g^+ \end{array}\right) \end{equation}
if 
$k_\perp \ne 0$ and
\begin{equation} \boldsymbol{M}(\boldsymbol{k})= \frac{1}{{2}} \left(\begin{array}{ccccc} 1 & 1 & 0 & 1 & -1\\ -1 & -1 & 0 & 1 & -1\\ 0 & 0 & 0 & 0 & 0\\ -s & s & 0 & -s & -s\\ s & -s & 0 & -s & -s\\ 0 & 0 & 0 & 0 & 0\\ 0 & 0 & -2\mathrm{i} s & 0 & 0 \end{array}\right) \end{equation}
if 
$k_\perp =0.$ Here,
\begin{equation} s=\frac{\vert k_3\vert}{k_3}={\pm} 1,\quad \omega_a^+=\frac{\omega_a}{\omega_{ag}},\quad \omega_g^+=\frac{\omega_g}{\omega_{ag}},\quad {\text{so that}\ \omega_a^{+2}+\omega_g^{+2}=1.} \end{equation}All the eigenvectors above have been normalized so as to give an orthonormal basis,
\begin{equation} \boldsymbol{Z}_\ell =\sum_{j=1}^7M_{j\ell}\boldsymbol{e}_j,\quad \boldsymbol{Z}_\ell^{\dagger}\boldsymbol{\cdot} \boldsymbol{Z}_m=\delta_{\ell m},\quad (\ell,m=1,2,\ldots,5) \end{equation}
where 
$(\boldsymbol {e}_1,\boldsymbol {e}_2,\ldots,\boldsymbol {e}_7)$ is the canonical basis of 
$\mathbb {C}^7,$ 
${\dagger}$ denotes the transpose and complex conjugation and 
$\delta _{\alpha \beta }$ denotes the Kronecker delta. The corresponding Fourier eigenfunctions are an orthogonal family. Therefore, we can decompose the vector 
$\hat {\boldsymbol {Y}}$ as
\begin{equation} \hat{\boldsymbol{Y}}=A_1^-\boldsymbol{Z}_1+A_1^+\boldsymbol{Z}_2+A_0\boldsymbol{Z}_3+A_2^-\boldsymbol{Z}_4+A_2^+\boldsymbol{Z}_5 \end{equation}where
$$\begin{gather} A_1^\mp=\frac{1}{\sqrt{2}k_\perp}\left (k_2\hat{u}_1-k_1\hat{u}_2\mp s\left (k_2\hat{b}_1-k_1\hat{b}_2\right)\right), \end{gather}$$
$$\begin{gather}A_0={-} \left (\mathrm{i} s\omega_a^+\hat{\vartheta}+\frac{k}{k_\perp}\omega_g^+\hat{b}_3\right), \end{gather}$$
$$\begin{gather}A^\mp_2={-}\frac{1}{\sqrt{2}}\left (\frac{k}{k_\perp}\hat{u}_3\mp\left (\frac{k}{k_\perp}s\omega_a^+\hat{b}_3-\mathrm{i}\omega_g^+\hat{\vartheta}\right)\right), \end{gather}$$
if 
$k_\perp \ne 0,$ and
$$\begin{gather} A_1^\mp=\tfrac{1}{{2}}\left (\hat{u}_1-\hat{u}_2\mp s\left (\hat{b}_1-\hat{b}_2\right)\right), \end{gather}$$
$$\begin{gather}A_0={-}\mathrm{i} s\hat{\vartheta}, \end{gather}$$
$$\begin{gather}A^\mp_2=\tfrac{1}{{2}}\left (\hat{u}_1+\hat{u}_2\mp s\left (\hat{b}_1+\hat{b}_2\right)\right), \end{gather}$$
if 
$k_\perp =0.$ The mode 
$A_0$ is associated with zero eigenvalue, while 
$A_1^\pm$ are associated with the eigenvalues 
$\pm \mathrm {i} \omega _a$ characterizing the Alfvén waves and 
$A_2^\pm$ are associated with the eigenvalues 
$\pm \mathrm {i} \omega _{ag}$ characterizing the magnetogravity waves. In addition, the modes 
$A_1^\pm$ characterizing Alfvén waves are expressed in terms of the horizontal components of the velocity and magnetic field, while the modes 
$A_2^\pm$ characterizing magnetogravity waves are expressed in terms of the buoyancy scalar and the vertical component of the velocity.
The quadratic part of the 
$L_2$ norm of MIPS (see (2.24))
\begin{equation} \varGamma_q(t)=\frac{1}{2}\int_{\mathcal{D}}\left (B_0\partial_z\vartheta+Nb_3\right)^2{\rm d}^3\boldsymbol{x}= \frac{1}{2}\sum_{\boldsymbol{k} \in {\mathcal{T}}} \left \vert \mathrm{i} \omega_a\hat{\vartheta}+ \omega_g\frac{k}{k_\perp}\hat{b}_3\right \vert^2, \end{equation}
can be expressed in terms of the spectral density of energy of the NP mode, denoted by 
${\mathcal {E}}^{(np)}(\boldsymbol {k},t),$ as
\begin{equation} \varGamma_q(t)=\frac{1}{2}\sum_{\boldsymbol{k} \in {\mathcal{T}}}\omega_{ag}^2\left \vert A_0(\boldsymbol{k},t)\right\vert^2=\sum_{\boldsymbol{k} \in {\mathcal{T}}}\omega_{ag}^2{\mathcal{E}}^{(np)}(\boldsymbol{k},t). \end{equation}
In the limit 
$F_r\rightarrow 0$ and 
$M\rightarrow 0,$ the 
$L_2$ norm of MIPS can be approximated by its quadratic part 
$\varGamma _q.$
In this case, one can consider the two inviscid invariants that are 
$E$ (total energy) and 
$\varGamma _q$ and study the system as a problem of equilibrium statistical mechanics. Indeed, statistical theory is particularly informative in the case of (non-magnetized) rotating stratified turbulence, since no inverse cascade appears in PST, whereas it is predicted in the presence of rotation, with the occurrence of a negative temperature states. This is possible because the NP mode including PV is affected by system rotation. On the other hand, the MIPS, that is the NP mode in the present case, is not affected by rotation, at least in the linear limit. Accordingly, application of the statistical theory, with and without rotation, should give no presumption of inverse cascade.
On the other hand, we note that, for unit magnetic and thermal Prandtl numbers, the transport equation for 
$\mathcal {E}^{(np)}(\boldsymbol {k},t),$ for example, reads
\begin{equation} \left (\partial_t+2\nu k^2\right){\mathcal{E}}^{(np)}(\boldsymbol{k},t)={\mathcal{N}}_{nl}^{(np)}(\boldsymbol{k},t), \end{equation}
where 
${\mathcal {N}}_{nl}^{(np)}$ denotes the nonlinear transfer term. Therefore, in the linear inviscid limit, the energy of the NP mode, as well as the energy of Alfvén waves and the energy of magnetogravity waves, do not evolve in time.
3.2. Case of small $F_r$ and finite $M$ (case (II))
For case (II), the non-zero elements of the matrix 
$\hat {\boldsymbol {\mathcal {L}}}(\boldsymbol {k})$ are
\begin{equation} \hat{\mathcal{L}}_{17}=\frac{k_1k_3}{k^2} N,\quad \hat{\mathcal{L}}_{27}=\frac{k_2k_3}{k^2}N,\quad \hat{\mathcal{L}}_{37}={-}\frac{k_\perp^2}{k^2}N,\quad \hat{\mathcal{L}}_{73}=N. \end{equation}
In this case, the spectrum of the matrix 
$\hat {\boldsymbol {\mathcal {L}}}(\boldsymbol {k})$ is found as
\begin{equation} \text{Sp}\,\hat{\boldsymbol{\mathcal{L}}}(\boldsymbol{k})=\{\lambda_1=\lambda_2=\lambda_3=0, \lambda_4={-}\mathrm{i}\omega_{g}, \lambda_5=\mathrm{i}\omega_{g}\}. \end{equation}
The associated eigenvectors, say 
$\boldsymbol {Z}_\ell$ (
$\ell =1,2,\ldots,5),$ such that 
$\hat {\boldsymbol {\mathcal {L}}}\boldsymbol {\cdot }{\boldsymbol {Z}}_\ell =\lambda _\ell \boldsymbol {Z}_\ell,$ are, respectively, the columns of the rectangular matrix given by (3.7) provided that 
$\omega _g^+$ is replaced by one and 
$\omega _a^+$ by zero. In this case, the normal modes 
$A_0$ and 
$A_1^\pm$ are associated with zero eigenvalues, whereas the modes 
$A_2^\pm$ are associated with the eigenvalues 
$\pm \mathrm {i} \omega _g$ characterizing the gravity waves,
$$\begin{gather} A_1^\mp= \frac{1}{\sqrt{2}k_\perp}\left (k_2\hat{u}_1-k_1\hat{u}_2\mp s\left (k_2\hat{b}_1-k_1\hat{b}_2\right)\right), \end{gather}$$
$$\begin{gather}A_0={-}\frac{k}{k_\perp}\hat{b}_3, \end{gather}$$
$$\begin{gather}A^\mp_2={-}\frac{1}{\sqrt{2}}\left (\frac{k}{k_\perp}\hat{u}_3\pm\mathrm{i}\hat{\vartheta}\right), \end{gather}$$
if 
$k_\perp \ne 0,$ and
$$\begin{gather} A_1^\mp= \tfrac{1}{{2}}\left (\hat{u}_1-\hat{u}_2\mp s\left (\hat{b}_1-\hat{b}_2\right)\right), \end{gather}$$
$$\begin{gather}A_0=0, \end{gather}$$
$$\begin{gather}A^\mp_2= \tfrac{1}{{2}}\left (\hat{u}_1+\hat{u}_2\mp s\left (\hat{b}_1+\hat{b}_2\right)\right), \end{gather}$$
if 
$k_\perp =0.$ It clearly appears that, for the case of small 
$F_r$ and finite 
$M$ (case (II)) for which the slow component of the linear operator is not zero (see (2.18)), the modes 
$A_2^\pm$ are associated with fast (gravity) waves, while 
$A_0$ and 
$A_1^\pm$ characterize slow modes.
Because in both the cases (I) and (II) the mode 
$A_0$ is associated with a zero eigenvalue, we will call it a NP mode. Likewise, the modes 
$A_1^\pm,$ which are associated with zero eigenvalues only in case (II) and not in case (I), we will call them slow Alfvén waves.
The present DNS results, which will be analysed in 
$\S \,4$, correspond more to case (II) than to case (I) (see figure 1). They characterize a weak MHD turbulence regime (with or without stratification) because the ratio of the nonlinear eddy turnover time to the Alfvén time, 
$\tau _{nl}/\tau _{a},$ exceeds unity (see figure 2a).
(a) Time evolution of horizontal Froude number 
$F_{rh}=U_h/(NL_h)$ for (non-magnetized) PST with 
$N=5$ and stratified MHD turbulence with 
$(B_0=0.4, N=5)$ and 
$(B_0=0.4, N=10)$ (runs B04N5 and B04N10, see table 2). (b) Time evolution of Alfvén–Mach number, 
$M=U_h/B_0,$ for MHD turbulence with or without stratification: 
$(B_0=0.4, N=2)$ 
$(B_0=0.4, N=5)$ and 
$(B_0=0.4, N=10)$ (runs B04N2, B04N5 and B04N10, see table 2) and 
$(B_0=0.4, N=0)$.
(a) Time evolution of the ratio of the nonlinear eddy turnover time to the Alfvén time, 
$\tau _{nl}/\tau _{a}=(B_0/L_v)/(U_b/L_h)$ for (non-magnetized) PST with 
$N=5$ and stratified MHD turbulence with 
$(B_0=0.4, N=5)$ and 
$(B_0=0.4, N=10)$ (runs B04N5 and B04N10, see table 2). (b) Variation of the ratio of gravity wave frequency to Alfvén wave frequency, 
$\omega _g/\omega _a=(k_\perp /k)N/(B_0\vert k_\parallel \vert )$ versus 
$k_\perp /k_\parallel$ for 
$(B_0=0.4,N=5)$ and several values of 
$k_\parallel$.
3.3. Asymptotic behaviour of energy ratios from linear theory
In the inviscid linear limit, an analytical solution for 
$\hat {\boldsymbol {Y}}(\boldsymbol {k},t)$ can be found (see Salhi et al. Reference Salhi, Baklouti, Godeferd, Lehner and Cambon2017). From such a solution, we deduce that for initial isotropic conditions with zero initial magnetic and buoyancy scalar fluctuations and a unit value for the magnetic and thermal Prandtl numbers (
$\nu =\eta =\kappa ),$ the spectral density of energies (kinetic, magnetic and potential) takes the form
$$\begin{gather} \mathcal{E}_\kappa(\boldsymbol{k},t)=\frac{1}{2}\left \vert \hat{\boldsymbol{u}}\right\vert^2=\frac{E_\kappa(k,0)}{8{\rm \pi} k^2} \left (\cos^2(\omega_at)+\cos^2(\omega_{ag}t)\right)\exp({-2\nu k^2t}), \end{gather}$$
$$\begin{gather}\mathcal{E}_m(\boldsymbol{k},t)=\frac{1}{2}\left \vert \hat{\boldsymbol{b}}\right\vert^2= \frac{E_\kappa(k,0)}{8{\rm \pi} k^2}\left (\sin^2(\omega_a t)+\omega_a^{+2}\sin^2(\omega_{ag}t)\right)\exp({-2\nu k^2t}), \end{gather}$$
$$\begin{gather}\mathcal{E}_p(\boldsymbol{k},t)=\frac{1}{2}\left \vert \hat{\vartheta}\right\vert^2=\frac{E_\kappa(k,0)}{8{\rm \pi} k^2}\left (\omega_g^{+2} \sin^2(\omega_{ag}t)\right)\exp({-2\nu k^2t}), \end{gather}$$
$$\begin{gather}E_\kappa(t)=\int_{\mathbb{R}^3}\mathcal{E}_\kappa \,{\rm d}^3\boldsymbol{k},\quad E_m(t)=\int_{\mathbb{R}^3}\mathcal{E}_m \,{\rm d}^3\boldsymbol{k},\quad E_p(t)=\int_{\mathbb{R}^3}\mathcal{E}_p \,{\rm d}^3\boldsymbol{k}, \end{gather}$$
where 
$E_\kappa (k,t)= 2{\rm \pi} \int _0^{\rm \pi} \mathcal {E}_\kappa (\boldsymbol {k},t)\sin \theta \,{\rm d}\theta,$ denotes the radial spectrum of kinetic energy and 
$\theta$ being the angle between the wavevector 
$\boldsymbol {k}$ and the vertical axis (of unit vector 
$\hat {\boldsymbol {z}}$, the same notation is used for the other fields).
It follows that, when 
$\omega _g/\omega _a\ll 1,$ which is the case when 
$k_\perp /\vert k_\parallel \vert \ll 1$ or when 
$\vert k_\parallel \vert$ is large (see figure 2b), the kinetic and magnetic energies behave similarly to their counterparts in non-stratified MHD turbulence,
\begin{equation} \mathcal{E}_\kappa(\boldsymbol{k},t)\propto (1+\cos (2\omega_at))\exp({-2\nu k^2t}),\quad {\mathcal{E}_m(\boldsymbol{k},t)\propto (1-\cos (2\omega_at))\exp({-2\nu k^2t})}, \end{equation}while the potential energy is relatively small,
\begin{equation} {\mathcal{E}_p(\boldsymbol{k},t)\propto (\omega_a^2/\omega_g^2) (1-\cos (2\omega_at))\exp({-2\nu k^2t})}. \end{equation}
When 
$\omega _g/\omega _a\gg 1,$ which is the case when 
$\vert k_\parallel \vert$ is small and 
$k_\perp /\vert k_\parallel \vert >1$ (see figure 2b), the energies exhibit a decaying oscillatory behaviour with frequency 
$\omega _g,$
\begin{equation} \left.\begin{gathered} \mathcal{E}_\kappa(\boldsymbol{k},t)\propto (3+\cos (2\omega_gt))\exp({-2\nu k^2t}),\\ \mathcal{E}_m(\boldsymbol{k},t),\ \mathcal{E}_p(\boldsymbol{k},t)\propto (1-\cos (2\omega_gt))\exp({-2\nu k^2t}). \end{gathered}\right\} \end{equation}
Additionally, theoretical insights can be obtained by analysing the ratios between the energy radial spectra. Indeed, it is found that, at small scales, 
$B_0k/N\gg 1,$ the local Alfvén ratio 
$E_\kappa (k,t)/E_m(k,t)$ approaches one indicating an equipartition of energy between the magnetic and kinetic components, whereas the ratios 
$E_p(k,t)/E_\kappa (k,t)$ and 
$E_p(k,t)/E_m(k,t)$ behave like 
$k^{-1}$ for long times and independently of the form of the initial spectrum 
$E_\kappa (k,0)$ (see Salhi et al. (Reference Salhi, Baklouti, Godeferd, Lehner and Cambon2017) and their figure 4). In counterpart, at large scales, 
$B_0k/N\ll 1,$ the ratio 
$E_p(k,t)/E_m(k,t)$ approaches one signifying that there is equipartition of energy between magnetic and potential components, while the local Alfvén ratio 
$E_\kappa (k,t)/E_m(k,t)$ approaches 
$2$ for long times.
By using the initial spectrum 
$E_\kappa (k,0)\propto k^n\exp (-k^2/k_p^2)$ with 
$n=2$ or 
$n=4$ and 
$\nu =\eta =\kappa =0.0025$ (as in the present DNS), we integrate numerically the solution (3.21) for several values of the ratio 
$B_0/(L_i N),$ where 
$L_i$ is the (isotropic) integral length scale (see table 1) and 
$k_p=6$ is the peak wavenumber. We found that the Alfvén ratio 
$E_\kappa (t)/E_m(t)$ approaches 
$2$ for small 
$B_0/(NL_i)$ and 1 for large 
$B_0/(NL_i).$ As for the ratio 
$E_\kappa (t)/E_p(t),$ it also approaches 
$2$ for small 
$B_0/(NL_i),$ while for large 
$B_0/(NL_i)$ it asymptotes to a large value as 
$t\gg 1.$ This is illustrated by figure 3 obtained for 
$B/(NL_i)=0.05, 2.9, 29.$
Initial isotropic conditions for all the runs (see table 2). Here, 
${L_i=2{\rm \pi} \int _0^{+\infty } (E_\kappa (k)/k)\,\textrm {d}k/E_\kappa }$, 
${E_\kappa =(1/2)\langle \boldsymbol {u}^2\rangle =\int _0^{+\infty }E_\kappa (k)\,\textrm {d}k},$ 
${u_{rms}=\sqrt {2E_\kappa }}$, 
${\tau _{nl}=L_i/u_{rms}}$, 
${R_e=u_{rms}L_i/\nu },$ 
${\varepsilon _\kappa =\nu \langle \partial _{x_j}u_i\partial _{x_j}u_i\rangle },$ 
${R_{e\lambda }=(2/3)E_\kappa \sqrt {15/(\nu \varepsilon _\kappa )}},$ 
$\nu$ and 
${\eta =\nu ^{3/4}/\varepsilon ^{1/4}}$ denote the (isotropic) integral length scale, kinetic energy, r.m.s. velocity, initial isotropic eddy turnover time, kinetic Reynolds number, kinetic dissipation rate, microscale Reynolds number, kinematic viscosity and Kolmogorov scale, respectively, and 
$\langle \boldsymbol {\cdot }\rangle$ denotes an ensemble average.
Linear theory predictions for energy ratios of (a) kinetic to magnetic, 
$E_\kappa (Nt)/E_m(Nt),$ (b) kinetic to potential, 
$E_\kappa (Nt) /E_p (Nt)$ for several values of 
$B_0/(L_iN)=0.05, 2.9, 29$ and 
$\nu =\eta =\kappa =0.0025.$ Initially, the turbulence is isotropic with zero magnetic and potential energies. The initial radial spectrum of the kinetic energy used is of the form 
$E(k,0)\propto k^n\exp (-k^2/k_p^2)$ with 
$n=2$ or 
$n=4$ and 
$k_p=6.$
As a result, we can conclude that the fact that the decay rates of the kinetic and magnetic energies are the same as for pure Alfvenic decay in the linear limit (see Moffatt Reference Moffatt1967) is at least due to the fact that the magnetic and thermal Prandtl numbers are equal to one. Indeed, the study by Sreenivasan & Maurya (Reference Sreenivasan and Maurya2021) obtains the decay of kinetic (
$E_\kappa (t)\sim t^{-1/2})$ and magnetic (
$E_m(t)\sim t^{-5/4})$ energy for frequencies 
$\omega _g \gg \omega _a$ with 
$\nu =\kappa =0$ and 
$\eta \ne 0$ through some approximations in this limit, so that, 
$E_\kappa (t)/E_m(t)\sim t^{3/4}$ for long times.
4. Direct numerical simulations for stratified MHD turbulence
The numerical simulations performed were started from isotropic initial conditions with zero magnetic and buoyancy fluctuations (see table 1). Eleven runs were performed, including two runs for non-stratified MHD cases for which 
$B_0 = 0.2$ and 
$B_0 = 0.4,$ three runs for purely stratified cases (i.e. PST) for which 
$N= 2,$ 
$N = 5$ and 
$N= 10,$ and six runs for stratified MHD cases which are listed in table 2. Due to the close similarities between the results obtained for 
$B_0 = 0.4$ and those obtained for 
$B_0 = 0.2,$ we mainly focus on the former. When useful we also present the results obtained for 
$B_0=0.2$.
Runs for stratified MHD turbulence with of a mean magnetic field. Initial values of the horizontal velocity 
$U_h$ and horizontal and vertical length scales 
$L_h$ and 
$L_v,$ horizontal and vertical Froude numbers 
$F_{rh}$ and 
$F_{rv}$ and buoyancy Reynolds number 
$R_{eb}=F_{rh}2R_e,$ as well as the initial value of the inverse Alfvén–Mach number 
$M^{-1}=B_0/u_{rms}$. The final time of the numerical simulations is 
$t_f=30$.
4.1. Initial isotropic conditions
The approach pursued in this work consisted in performing DNS of the fully nonlinear equations of motion (2.1) here solved numerically by using the pseudospectral SNOOPY code (Lesur & Longaretti Reference Lesur and Longaretti2007). The domain of integration is a triple-periodic cubic box with edge 
$L_0=2{\rm \pi}$ consisting of 
$512^3$ grid points. Aliasing errors are filtered out through the implementation of the so called 2/3 rule and as a result the minimum and maximum wavenumbers are 
$2/3$ dealiasing rule and as a result the minimum and maximum wavenumbers are 
$k_{min} = 1$ and 
$k_{max}\approx 170,$ respectively. The simulations are started with isotropic initial conditions of purely hydrodynamic turbulence generated with a separate run, consisting of velocity fluctuations whose spectral distribution follows the initial kinetic energy spectrum,
\begin{equation} E_{\kappa}(k,0)\propto k^4 \exp\left ({-}k^2/k^2_p\right). \end{equation}
In other words, in order to let turbulence develop a well extended inertial range and strong enough nonlinear energy transfers, a precomputation of the isotropic case is done for approximately one eddy-turnover time (see Salhi et al. Reference Salhi, Jacobitz, Schneider and Cambon2014) to generate the initial condition, then mean magnetic field and mean density gradient are turned on in the production runs. In this preliminary stage only, large-scale forcing is applied until a statistical steady state of classical Kolmogorov-like isotropic turbulence is reached. This forcing is applied in a shell of wavenumbers corresponding to 
$1\leqslant k\leqslant 4.$ The different parameters obtained at the end of this precomputation phase are summarized in table 1. We note that the resulting precomputation kinetic energy radial spectrum (RPS) 
$E_\kappa (k,0)$ (see figure 4a) exhibits indeed a 
$-5/3$ power law scaling in the range 
$3\leqslant k \leqslant 20$. As the dynamics is then let to evolve, the kinetic energy spectrum 
$E_\kappa (t)$ during the freely decaying homogeneous isotropic turbulence (HIT) phase remains close to the Kolmogorov prediction 
$E_\kappa (t)\sim t^{-10/7}$ (see figure 5c).
(a) The radial spectrum 
$E(k,0)\propto k^4\exp (-k^2/k_p^2)$ with 
$k_p=6$ compared with the resulting RPS 
$E(k,0)$. (b) Time evolution of the compensated total energy 
$E(t)$ with respect to 
$t^{-0.94}$ in the case of non-stratified MHD turbulence of 
$B_0=2.$
(a) Time evolution of total (
$\textrm {kinetic}+\textrm {magnetic}+\textrm {potential}$) energy normalized by its initial value for the runs B04N2, B04N5 and B04N10. The DNS results for the non-stratified MHD case of 
$B_0=0.4$ are also reported for comparison. (b) Time evolution of kinetic, magnetic, potential energies normalized by the initial value of the kinetic energy for the run B04N2. (c) Time evolution of total (
$\textrm {kinetic}+\textrm {potential}$) energy normalized by its initial value for PST with 
$N=2,$ 
$N=5$ and 
$N=10.$ The DNS results for HIT are also reported for comparison. (d) Time evolution of kinetic, magnetic, potential and total dissipation rates normalized by the initial value of the kinetic dissipation rate for the run B04N2.
4.2. A weak stratified MHD turbulence regime
In figure 1 we show the time evolution of horizontal Froude number, 
${F_{rh}=U_h/(L_hN)}$ (figure 1a) and Alfvén–Mach number 
$M=U_h/B_0$ (figure 1b) for the runs B04N5 and B04N10, as well as those for (non-magnetized) PST with 
$N=5$ and those for non-stratified MHD case of 
$B_0=0.4.$ As it can be seen, the present DNS results rather correspond to case (II). Here,
$$\begin{gather} U_{h}=\sqrt{2E_{\kappa h}},\quad E_{\kappa h}=\frac{1}{2}\sum_{\boldsymbol{k} \in {\mathcal{T}}} \left(\vert \hat{u}_1\vert^2+\vert \hat{u}_2\vert^2\right), \end{gather}$$
$$\begin{gather}L_h={2{\rm \pi}}E_\kappa^{{-}1}\sum_{\boldsymbol{k} \in {\mathcal{T}}}\frac{E_\kappa(k_\perp)}{k_\perp},\quad L_v={2{\rm \pi}}E_\kappa^{{-}1}\sum_{\boldsymbol{k} \in {\mathcal{T}}}\frac {E_\kappa(k_\parallel)}{k_\parallel}, \end{gather}$$
denotes the r.m.s. horizontal velocity, the horizontal kinetic energy (the same notation is used for the magnetic field), and 
$L_h$ and 
$L_v$ are the horizontal and vertical length scales, respectively.
Indeed, both the horizontal and vertical (
$F_{rv}=U_h/(L_vN))$ Froude numbers decrease rapidly with time, becoming very small at the final time of the numerical simulation because 
$U_h$ decreases while oscillating (see figure 1b) whereas 
$L_h$ (respectively, 
$L_v)$, which initially takes the value 
$L_h = 2.32$ (respectively, 
$L_v=3.16,$ see table 2), increases with time, exhibiting oscillations, and appears by asymptoting to a limit value 
$(< 6).$ For the run B04N5, 
$F_{rh} = 0.01$ and 
$F_{rv}= 0.08$ at 
$t=1,$ and 
$F_{rh} =0.002$ and 
$F_{rv}= 0.003$ at 
$t = 20.$
Regarding the buoyancy Reynolds number 
$R_{eb}=U_h^3/(\nu N^2 L_h)$, which constitutes a primary control parameter in characterizing the regime of purely stratified flows (see, e.g. Billant & Chomaz Reference Billant and Chomaz2001; Brethouwer et al. Reference Brethouwer, Billant, Lindborg and Chomaz2007; Portwood et al. Reference Portwood, de Bruyn Kops, Taylor, Salehipour and Caulfield2016; Chini et al. Reference Chini, Michel, Julien, Rocha and Colm-cille2022), it decreases rapidly becoming less than unity for 
$t>3$ (not shown). In addition, according to our DNS results there is a weak MHD turbulence regime (with or without stratification) because the ratio of the nonlinear eddy turnover time to the Alfvén time, 
$\tau _{nl}/\tau _{a}=(B_0/L_v)/(U_b/L_h),$ which increases with time, exceeds unity for 
$t>4$ as shown by figure 2(a).
Therefore, we may conclude that the regime of the stratified MHD cases studied in the present work rather corresponds to a weak stratified MHD turbulence regime affected by viscosity.
4.3. Temporal evolution of global quantities
4.3.1. Decay of total energy and dissipation rate
In this section, we examine the temporal evolution of total (
$\textrm {kinetic}+\textrm {magnetic}+\textrm {potential}$) energy, 
$E=E_\kappa +E_m+E_p$ and total dissipation rate 
${\varepsilon _t= \varepsilon _\kappa +\varepsilon _m+\varepsilon _p}$ where 
$\varepsilon _\kappa =\nu \sum _{\boldsymbol {k}\in {\mathcal {T}}}k^2\vert \hat {\boldsymbol {u}}\vert ^2$ (the same notation is used for the other fields). Recall that in the present study, the magnetic and thermal Prandtl numbers are taken equal to unity (
$\nu =\eta =\kappa ).$
The role of a uniform magnetic field in MHD turbulence has been widely explored in the literature (see e.g. Verma Reference Verma2004; Bigot, Galtier & Politano Reference Bigot, Galtier and Politano2008a; Briard & Gomez Reference Briard and Gomez2018): it makes MHD turbulent plasmas anisotropic and slows down the energy decay by reducing the transfer along the mean magnetic field. It is worth noticing that two-dimensional turbulence is more representative than isotropic MHD turbulence with no mean magnetic field (or than MHD turbulence with a weak mean magnetic field), when it comes to the investigation of strong MHD turbulence, for which the energy cascade occurs preferentially in the direction perpendicular to the mean magnetic field. In balanced strong MHD turbulence (i.e. with zero cross-helicity) Cho, Lazarian & Vishniac (Reference Cho, Lazarian and Vishniac2002) found that the total energy decay follows a power law 
$E(t)\sim t^{\gamma _m}$ where 
$\gamma _m$ is very close to unity. Bigot, Galtier & Politano (Reference Bigot, Galtier and Politano2008b) predicted the slowing down of the energy decay due to anisotropy in the limit of a strong mean magnetic field, with distinct power laws for energy decay of shear-Alfvén waves 
$({\sim }t^{-2/3})$ and pseudo-Alfvén waves 
$({\sim }t^{-1}).$ According to Banerjee & Jedamzik (Reference Banerjee and Jedamzik2004), the total energy in isotropic MHD turbulence decays as 
${\sim }t^{-1.33}$ for an initial 
$k^4$-spectrum. For the non-stratified MHD cases of 
$B_0=0.2$ and 
$B_0=0.4,$ our numerical simulations indicate that the energy decay exponent is close to 
$1.25$ (see figure 5a), which is lower than for isotropic MHD turbulence. This is due to the presence of the mean magnetic field, having an effect even when its intensity is weak. In contrast, for the non-stratified MHD case with 
$B_0=2$ (not considered here), the total energy decay rate is close to 
$0.94$ (see figure 4b), which is not far from that (
${\sim }1$) found in the case of balanced strong MHD turbulence (see Cho et al. Reference Cho, Lazarian and Vishniac2002).
In figure 5(a) we show the temporal evolution of 
$E(t)$ for the runs B02N2, B02N5 and B02N10 as well as for the non-stratified MHD case with 
$B_0= 0.4.$ As it can be seen, the energy level increases as 
$N$ increases, and the total energy decrease is clearly slower for the stratified MHD cases than for the non-stratified MHD case. Note that the decay of 
$E(t)$ depends on 
$B_0$ and 
$N$ and not only on the ratio 
$B_0/(L_iN).$ Indeed, for runs B02N5 and B04N10, the ratio 
$B_0/(L_iN)= 0.023$ is the same while the energy level for the run B04N10 is higher than that of the run B02N5 (not shown). On the other hand, it appears that the decay of 
$E(t),$ as well as the decay of the kinetic, magnetic and potential energies which show a decaying oscillation (see figure 5b), tend to follow the power law decay 
${\sim }t^{-\gamma _s}$ with 
$\gamma _s=0.95$ thus with the same decaying trend of the total energy (
$\textrm {kinetic}+\textrm {potential}$) for the PST cases shown in figure 5(c). This temporal decay exponent is higher (respectively, lower) in absolute value than that predicted by Davidson (Reference Davidson2010) for Saffman turbulence 
$\gamma _s=4/5=0.8$ (Bachelor turbulence, 
$\gamma _s=8/7\approx 1.14)$ at high-
$R_{eb}.$ In their study of stratified homogeneous turbulence using a two-point closure EDQNM statistical model and DNS, Staquet & Godeferd (Reference Staquet and Godeferd1998) reported that energy decay rates are close to 
$1.$
On the other hand, the fact that for the stratified MHD cases the temporal decay exponents of kinetic, magnetic and potential energies are the same (see figure 5b), in agreement with linear theory (LT) predictions for 
$\nu =\eta =\kappa =1$ (see §§ 3.3 and 4.3.2) would be due to the fact that the magnetic and thermal Prandtl numbers are unity.
In figure 5(d) we plot the time evolution of the kinetic, magnetic and potential dissipation rates 
$\varepsilon _\kappa (t),$ 
$\varepsilon _m(t)$ and 
$\varepsilon _p(t)$ as well as the total dissipation 
$\varepsilon _t (t)$ for the run B04N2. As it can be seen, 
$\varepsilon _m$ (respectively, 
$\varepsilon _p)$ which is initially zero, increases with time, due to the generation of small-scale magnetic (potential) fluctuations, and reaches a maximum and then decreases while oscillating. The kinetic dissipation rate 
$\varepsilon _\kappa$ shows decaying oscillations. Note that 
$\varepsilon _p$ remains less than 
$\varepsilon _m$ and 
$\varepsilon _k$ for all times. This means that the energy dissipation by the viscous effects or by Joules effects are more important than the energy dissipation by thermal diffusive effects. The total dissipation rate 
$\varepsilon _t$ shows a slight growth, which caused by the initial increase in 
$\varepsilon _m,$ before decreasing. It tends to follow the power law decay 
${\sim }t^{-1.95}.$
4.3.2. Linear theory versus DNS for energy components
In this section, we compare the LT predictions for the kinetic energy 
$E_\kappa (t),$ magnetic energy 
$E_m(t)$ and potential energy 
$E_p(t)$ with the DNS results for the stratified MHD case with 
$B_0=0.4$ and 
$N=10,$ which has the smallest value of the initial 
$F_{rh}=U_h/(L_h N)=0.047.$ The time evolution of 
$E_\kappa (t),$ 
$E_m(t)$ and 
$E_m(t)$ yielded by LT
$$\begin{gather} E_\kappa(t)=\frac{1}{4}\int_0^{+\infty}E_\kappa(k,0)\left[\int_0^{\rm \pi} \left (\cos^2(\omega_a t)+\cos^2(\omega_{ag}t)\right) \sin\theta \,{\rm d}\theta \right]\exp({-2\nu k^2t})\,{\rm d}k \end{gather}$$
$$\begin{gather}E_m(t)=\frac{1}{4}\int_0^{+\infty}E_\kappa(k,0)\left[\int_0^{\rm \pi} \left (\sin^2(\omega_a t)+\omega_a^{+2}\sin^2(\omega_{ag}t)\right) \sin\theta \,{\rm d}\theta \right]\exp({-2\nu k^2t})\,{\rm d}k \end{gather}$$
$$\begin{gather}E_p(t)=\frac{1}{4}\int_0^{+\infty}E_\kappa(k,0)\left[\int_0^{\rm \pi} \omega_g^{+2}\sin^2(\omega_{ag} t)\sin\theta \,{\rm d}\theta \right]\exp({-2\nu k^2t})\,{\rm d}k \end{gather}$$
have been obtained numerically by using the RPS (see the beginning of § 4 and figure 4a) for 
$E_\kappa (k,0),$ so that the initial conditions are the same for both DNS and LT.
Figure 6(a) shows the evolution of 
$E_\kappa (t),$ 
$E_m(t)$ and 
$E_p(t)$ yielded by LT versus 
$Nt/{\rm \pi}$. As it can be seen, for 
$Nt/{\rm \pi} > 6,$ the kinetic energy, magnetic energy and potential energy exhibit parallel evolution. The decay rate changes with the elapsed time. For instance, at 
$Nt/{\rm \pi} >100,$ the energy decay tends to follow 
${\sim }t^{-1.25}.$ Note that, according to LT, the decay rate in the long-time limit depends on the shape of 
$E_\kappa (k,0)$ near 
$k = 0$ so that if 
$E_\kappa (k,0)\propto k^n$ then the total energy 
$E(t)=E(0) (1+2\nu k_p^2t)^{-(n+1)/2}$ as well as the components 
$E_\kappa (t),$ 
$E_m(t)$ and 
$E_p(t)$ behave like 
${\sim }t^{-(n+1)/2},$ as shown by the inset of figure 6(a) obtained for 
$n=2$ and 
$n=4$ (see also Hanazaki & Hunt (Reference Hanazaki and Hunt1996) for the pure stratified case). Recall that, for 
$\omega _g\gg \omega _a\gg \eta k^2$ and 
$\nu =\kappa =0,$ LT predicts 
$E_\kappa \sim t^{-1/2}$ and 
$E_m\sim t^{-5/4}$ for 
$t\gg 1$ (see Sreenivasan & Maurya Reference Sreenivasan and Maurya2021).
(a) Linear theory predictions. Time evolution of the kinetic energy 
$E_\kappa (t),$ magnetic energy 
$E_m(t)$ and potential energy 
$E_p(t)$ normalized by 
$E_\kappa (0)$ for the case with 
$B_0=0.4,$ 
$N=10,$ so that the initial Froude number is approximately 
$F_{rh}=U_h/(L_h N)=0.047.$ The inset shows the time evolution of 
$E_\kappa (t)$ yielded by LT where 
$E_\kappa (k,0)\propto k^n\exp ({-(k^2/k_p^2)})$ with 
$n=2,$ 
$n=4$ and 
$k_p=6$. (b–d) Comparison between DNS results (run B04N10) and LT predictions for the time evolution of the kinetic energy, magnetic energy and potential energy, respectively.
In figure 6(b–d) we compare the results of the LT for the energy components with those of the DNS, recalling that for this comparison we have the same initial conditions for LT and for DNS. We observe that, for 
$Nt>2{\rm \pi},$ LT overestimates the energy components. The faster decay of energy components in DNS is due to the nonlinear effect. Note that the differences between the LT predictions and the DNS results for the kinetic energy and magnetic energy are more significant in the non-stratified MHD case than in the stratified MHD case. Therefore, we can conclude that, compared with the non-stratified MHD case, the nonlinear interactions are reduced due to the presence of the stable stratification, at least for the values of 
$N$ considered in the present study (
$N\leqslant 10).$ This explains why the decay of the energy components in the present DNS is not quasilinear.
4.3.3. Vertical components of energies
In figure 7(a) we plot the variation of twice the vertical kinetic energy 
$2E_{kv}=\langle u_3^2\rangle$ versus the dimensionless time 
$Nt/{\rm \pi}$ for the runs B04N2, B04N5 and B04N10. The results obtained for PST with 
$N=2,$ 
$N=5$ and 
$N=10$ are also reported for comparison. We observe that for fixed 
$(Nt/{\rm \pi} )>0),$ 
$E_{kv}(t)$ increases as 
$N$ increases, and the profiles of 
$E_{kv}(t)$ for the stratified MHD cases and PST look similar with decaying oscillations whose periods are approximately the same (
${\sim }{\rm \pi} /N).$ For PST the oscillations are more damped. Note that for the non-stratified MHD case of 
$B_0=0.4$ or 
$B_0=0.2,$ 
$E_{\kappa v}(t)$ exhibits decaying oscillations with an oscillation period approximately 
$L_v(0)/B_0\sim {\rm \pi}/B_0\sim 7.85>{\rm \pi} /N$ for the values of the governing parameters used in the simulation of the present study.
Variation of twice the vertical kinetic energy (a) and potential energy (textitb) versus the dimensionless time 
$Nt/{\rm \pi}$ for the runs B04N2, B04N5 and B04N10. The DNS results for PST with 
$N=2, 5, 10$ are also reported for comparison.
The vertical kinetic and potential energies evolve in similar ways but with strong oscillatory energy exchanges. The potential energy 
$E_p(t),$ which is initially zero, increases with time and reaches a maximum at approximately 
$t={\rm \pi} /(2N)$ then decays while oscillating as shown by figure 7(b). For fixed 
$(Nt/{\rm \pi} )>0,$ 
$E_p$ increases as 
$N$ increases. The decaying oscillations around zero emerging in the buoyancy flux represent net exchanges from the vertical component of the kinetic energy to potential energy (positive values) and back from potential into kinetic energy (negative values).
It clearly appears that, at least, for 
$0< B_0\leqslant 0.4$ and 
$0< B_0/(L_iN) \leqslant 0.12,$ the profiles of the vertical kinetic energy, 
$E_{\kappa v}(t),$ the potential energy, 
$E_p(t)$ and the buoyancy flux, 
$-\langle u_3\vartheta \rangle$ (not shown) for the stratified MHD cases and the purely stratified cases look very similar with decaying oscillations whose periods are approximately the same 
${\sim }{\rm \pi} /N.$
Concerning the vertical magnetic energy 
$E_{mv}(t),$ it is significantly affected by the presence of the buoyancy force. It increases from zero, reaches a maximum and then decays with time while oscillating. The oscillation period is approximately 
${\sim }L_v(0)/B_0\sim {\rm \pi}/B_0$ in the case without stratification and approximately 
${\rm \pi} /N (< L_v(0)/B_0)$ in the presence of the buoyancy force. The effect of stratification on the vertical magnetic energy is conspicuous since 
$E_{mv}$ is drastically reduced as 
$N$ increases (see figure 8a). In the presence of stratification, the contribution of the vertical magnetic energy to the (total) magnetic energy is very low: at large time, it does not exceed 
$1\,\%$ as shown by figure 8(b). The significant reduction of the intensity of vertical magnetic fluctuations in the presence of stratification is also illustrated by figure 9, representing snapshots of physical-space surface plots at 
$t=5.$ In the stratified MHD case, the flow appears smooth, and the decaying magnetic field has a similar structure.
(a) Variation of twice the vertical magnetic energy versus the dimensionless time 
$Nt/{\rm \pi}$ for the runs B04N2, B04N5 and B04N10. (b) Time evolution of the ratio of vertical magnetic energy to (total) magnetic energy 
$E_{mv}(t)/E_m(t)$ for the runs B04N5 and B04N10. The DNS results for the non-stratified MHD case of 
$B_0=0.4$ are also reported for comparison.
Snapshots of physical-space surface plots at 
$t=5$ for vertical magnetic fluctuations: (a) non-stratified MHD case with 
$B_0=0.4$; (b) stratified MHD case with 
$B_0=0.4$ and 
$N=10$ (run B04N10).
4.3.4. Horizontal components of energies
For non-stratified MHD cases, the horizontal kinetic energy 
$E_{\kappa h}(t)$ shows decaying oscillations with an oscillation period of approximately 
$L_v(0)/B_0\sim {\rm \pi}/B_0 .$ The horizontal magnetic energy 
$E_{\kappa h}(t)$, initially zero, increases and reaches a maximum then decreases while oscillating with the same oscillation period 
${\sim }{\rm \pi} /B_0.$ When the buoyancy force and the mean magnetic field are simultaneously present, an increase in the energy levels of 
$E_{\kappa h}(t)$ (and also of 
$E_{mh}(t)$) is observed. However, the profiles of 
$E_{\kappa h}(t)$ (and those of 
$E_{mh}(t)$) with or without stratification look almost the same. For instance, the maxima (respectively, the minima) of decreasing oscillations shown in the development of 
$E_{\kappa h}(t)$ (respectively, 
$E_{mh}(t)$) are reached at times which are the same in the presence or not of the buoyancy force as shown in figure 10(a) (respectively, figure 10b). On the other hand, we note that the time evolution of the kinetic-magnetic flux, which represents the exchanges from the horizontal kinetic energy to horizontal magnetic energy, is not strongly affected by the buoyancy force (not shown). Therefore, we may conclude that the horizontal components of kinetic and magnetic energies are more sensitive to the effect of the mean magnetic field than the vertical components. On the other hand, we note that the presence of the buoyancy force can make the flow strongly anisotropic (at the magnetic fluctuations) because it induces a significant reduction of 
$E_{mv}(t)$ and a small increase of 
$E_{mh}(t).$ This point will be examined further in § 4.4.
Time evolution of the horizontal kinetic (a) and magnetic (b) energies for the runs B04N2, B04N5 and B04N10. The DNS results for the non-stratified MHD case of 
$B_0=0.4$ as well as those for PST with 
$N=10$ are also reported for comparison.
4.3.5. Energy of waves
As shown in §§ 3.1 and 3.2, the modes 
$A_1^\pm$ have the same expression in both the two cases (I) and (II) (see (3.12a) and (3.19a)). In the former case, they are associated with fast (Alfvén) waves while in the latter case they are slow modes. The energy of these two modes is denoted by 
$E^{(aw)}(t)=\frac {1}{2}\sum _{\boldsymbol {k} \in {\mathcal {T}}} ( \vert A_1^-\vert ^2+ \vert A_1^+\vert ^2)$, where the superscript 
$(aw)$ stands for Alfvén waves. In counterpart, the modes 
$A_2^\pm,$ which characterize fast (gravity) waves in case (II) or fast magnetogravity waves in case (I) as well as the mode 
$A_0$ (NP mode) do not have the same expression in the two cases. Because our DNS results rather correspond to case (II), we consider the expression given by (3.19a,b) for the modes 
$A_0$ and 
$A_1^\pm,$ so that, the energy of these modes takes the form
$$\begin{gather} E^{(np)}(t)=\frac{1}{2}\sum_{\boldsymbol{k} \in {\mathcal{T}}}\left \vert A_0\right\vert^2= \frac{1}{2}\sum_{\boldsymbol{k} \in {\mathcal{T}}}\frac{k^2}{k_\perp^2}\left \vert \hat{b}_3\right \vert^2, \end{gather}$$
$$\begin{gather}E^{(gw)}(t)=\frac{1}{2}\sum_{\boldsymbol{k} \in {\mathcal{T}}}\left (\left \vert A_2^-\right\vert^2+\left \vert A_2^+\right\vert^2\right)=\frac{1}{2}\sum_{\boldsymbol{k} \in {\mathcal{T}}}\left (\frac{k^2}{k_\perp^2} \left \vert \hat{u}_3\right\vert^2+\left \vert \hat{\vartheta}\right\vert^2\right). \end{gather}$$
Note that the sum, 
$E^{(np)}(t)+E^{(gw)}(t),$ has the same expression in both cases (I) and (II).
The energies 
$E^{(aw)}$ and 
$E^{(gw)}$ , which are initially equal, decrease with time and tend to follow the power law decay 
${\sim }t^{-0.95}$ as for the total energy 
$E(t)$ (not shown). However, the contribution to the total energy coming from 
$E^{(aw)}$ is more important than that coming from 
$E^{(gw)}$ especially for the cases with 
$N=5$ and 
$N=10$ as shown by figure 11(a).
(a) Time evolution of the ratio between the energy of Alfvén waves and total energy 
$E^{(aw)}(t)/E(t)$ and the ratio between the energy of gravity waves and total energy 
$E^{(gw)}(t)/E(t)$. (b) Variation of the energy of the NP mode, 
$E^{(np)}(t),$ normalized by 
$E^{(aw)}(0)$ for the runs B04N2, B04N5 and B04N10.
In figure 11(b) we show the time evolution of the energy of the NP mode normalized by half of the initial kinetic energy 
$E_\kappa (0).$ Because the magnetic and density fluctuations are initially zero, the energy 
$E^{(np)}(t)$ is then initially zero. The profile of 
$E^{(np)}(t)$ closely resembles that of the vertical magnetic energy (see (4.4)): 
$E^{(np)}(t)$ increases until reaching a maximum at 
$Nt/{\rm \pi} \sim 0.5$ then decreases oscillating with a period of approximately 
${\rm \pi} /N,$ and its maximum of 
${E}^{(np)}(t)$ significantly decreases as 
$N$ increases. Compared with 
$E^{(aw)}$ or with 
$E^{(gw)},$ 
$E^{(np)}(t)$ remains very small. The energy ratios 
$E^{(np)}(t)/E^{(gw)}$ and 
$E^{(np)}(t)/E(t)$ do not exceed 
$0.27$ and 
$0.1,$ respectively. From figure 11(b), it is clear that NP mode is not a stagnant mode. It has probably become dependent on the wave modes, as indicated by one of the referees.
Note that the dependence of 
$E^{(aw)}(t),$ 
$E^{(gw)}(t)$ and 
$E^{(np)}(t)$ on 
$B_0$ and 
$N$ emerging in the DNS results is due to the nonlinear transfer terms (see (3.16)).
We now examine the behaviour of the quadratic part of the 
$L_2$ norm of MIPS,
\begin{equation} \varGamma_q(t)=\frac{N^2}{2}\sum_{\boldsymbol{k} \in {\mathcal{T}}} \left \vert \hat{b}_3+\mathrm{i} F_rM^{{-}1}\left (k_\parallel L\right)\hat{\vartheta}\right\vert^2. \end{equation}
Because the vertical magnetic fluctuations are strongly affected by stratification, so that the vertical magnetic energy remains very small compared with horizontal magnetic energy as well as to potential and kinetic energy components, the two terms on the right-hand side of (4.5) can have the same order of magnitude, and hence, for small 
$F_r$ and finite 
$M$ (case (II)) it is not justified to neglect the second term. This is confirmed by the present DNS results as shown by figure 12(a) obtained for the run B04N5. The figure shows the time development of 
$\varGamma _q(t)$ and its counterpart obtained by neglecting the second term on the right-hand side of (4.5),
\begin{equation} \varGamma^{(a)}_q=\varGamma_q\left (F_r\ll 1, M\sim {O}\left (1\right)\right)= N^2E_{mv}(t). \end{equation}Variation of the quadratic part 
$\varGamma _q$ of the 
$L_2$ norm of MIPS (see (4.5)) and its counterpart 
$\varGamma _q^{(a)}$ obtained in the limit of small 
$F_r$ and 
$M\sim {O}(1)$ (case (II), see (4.6)) versus the dimensionless time 
$Nt/{\rm \pi}$ for the run B04N5. (b) Time-evolution of 
$\varGamma _q(t)$ normalized by 
$\varGamma (t)$ (the 
$L_2$ norm of MIPS) for the run B04N5.
As for the contribution of the cubic and quartic parts to the 
$L_2$ norm of MIPS, these are not very significant compared with the contribution of the quadratic part 
$\varGamma _q,$ at least for the values of the governing parameters used in the simulations of the present study as illustrated by figure 12(b) displaying the time evolution of 
$\varGamma _q(t)/\varGamma (t)$ for the run B04N5.
4.4. Spectra
In this section our analysis focuses on spectra, which provide information on the distribution of energy across the different scales. As shown in previous studies, in stratified turbulence it is more convenient to consider horizontal and vertical spectra separately because of the anisotropy of the flow (see, e.g. Dewan Reference Dewan1997; Lindborg Reference Lindborg2006): stratification is often considered to significantly suppress vertical turbulent motions at scales larger than the Ozmidov scale, 
$L_O=\sqrt {\varepsilon _\kappa /N^3}$ (see, e.g. Sagaut & Cambon Reference Sagaut and Cambon2008).
4.4.1. Total energy
The radial spectrum 
$E(k,t)$ of total energy at 
$t=5$ is plotted in figure 13 for the run B04N10 as well as for the non-stratified MHD case of 
$B_0=0.4$ and for PST of 
$N=10.$ The dashed line indicates the 
$k^{-5/3}$ power-law scaling; our inertial range is, however, a bit too short to perform a fit of this spectrum. Although for 
$t<2$ and at 
$4\leqslant k\leqslant 15,$ 
$E(k)$ looks similar to such power law scaling, it becomes more and more steep as time increases.
(a) Radial spectrum of total energy at 
$t=5$ for the run B04N10. The DNS results for PST with 
$N=10$ and for the non-stratified MHD case of 
$B_0=0.4$ are also reported for comparison. (b) Vertical wavenumber spectra of the horizontal magnetic energy at 
$t=5$ for the runs B02N2 and B04N2. The DNS results for the non-stratified MHD cases of 
$B_0=0.2$ and 
$B_0=0.4$ are also reported for comparison.
Compared with the case without stratification, there is a decrease in 
$E(k)$ at high wavenumbers and an increase at low wavenumbers when the buoyancy force is present. We note that the profiles of 
$E(k)$ for the stratified MHD cases with 
$B_0=0.2$ and 
$B_0=0.4$ remain very close to their counterparts obtained for PST, in agreement with the analysis of the development of the total energy presented in § 4.3.1.
4.4.2. Horizontal kinetic and magnetic energies
In figure 13(b) we plot the vertical wavenumber spectra of the horizontal magnetic energy 
$E_{mh}(k_\perp )=\frac {1}{2}\sum _{k_\perp } ( \vert \hat {b}_1\vert ^2+ \vert \hat {b}_2\vert ^2)$, versus the horizontal wavenumber 
$k_\perp$ at 
$t=5$ for the runs B02N10 and B04N10. The results obtained in the non-stratified MHD cases of 
$B_0=0.2$ and 
$B_0=0.4$ are reported for comparison. A similar behaviour is found for the vertical wavenumber spectra of the horizontal kinetic energy 
$E_{\kappa h}(k_\perp ).$ As we can see, at almost all scales, 
$E_{m h}(k_\perp )(k_\perp )$ (or 
$E_ {\kappa h}(\perp )(k_\perp )$) obtained in the case 
$B_0=0.4$ remains greater than that obtained in the case 
$B_0=0.2.$ Similar results are found for 
$E_{\kappa h}(k_\parallel )$ and 
$E_ {mh }(k_\parallel )$ (not shown). We notice that, at 
$4\leqslant k_\perp \le 10,$ the decay of 
$E_{m h}(k_\perp )$ (or 
$E_{\kappa h}(k_\perp )$) as well as the vertical wavenumber spectra of the vertical kinetic energy 
$E_{\kappa v}(k_\perp )$ (see figure 14a) tend to follow the power law 
${\sim }k^{-3}_\perp$ and this for a wide range of time 
$3 \leqslant t\leqslant 30.$
Vertical wavenumber spectra of the vertical kinetic energy (a) and vertical magnetic energy (b) at 
$t=5$ for the runs B02N2 and B04N2. The DNS results for the non-stratified MHD cases of 
$B_0=0.2$ and 
$B_0=0.4$ are also reported for comparison.
4.4.3. Vertical kinetic and magnetic energies
In agreement with the results presented in figure 8(a) which indicate that the maximum of the vertical kinetic energy increases as 
$N$ increases, at almost all scales, the vertical wavenumber spectra of the vertical kinetic energy 
$E_{\kappa v} (k_\perp )$ increases as 
$N$ increases (not shown). To illustrate the effect of the 
$\boldsymbol {B}_0$ intensity, in figure 14(a), we present 
$E_{\kappa v}(k_\perp )$ versus 
$k_\perp$ for the runs B02N2 and B04N2 and 
$t=5.$ As one can see, at all scales, there is a significant increase of 
$E_{\kappa v}(k_\perp )$ for the case 
$B_0= 0.4$ compared with the case 
$B_0= 0.2.$ Similar results are found for the potential energy spectra.
A significant decrease of the vertical magnetic energy due to the presence of the buoyancy force occurs at all horizontal scales, as shown in figure 14(b) which displays 
$E_{mv}(k_\perp )$ versus 
$k_\perp$ for the runs B02N10 and B04N10. The figure also reveals that the difference between the spectra is more pronounced with stratification than without it. On the other hand, we note that, at all horizontal scales, 
$E_{mv}(k_\perp )$ decreases as 
$N$ increases (not shown). Conversely, the decrease in the vertical wavenumber spectra 
$E_{mv}(k_\parallel )$, as 
$N$ increases, occurs almost at large vertical scales where 
$E_{mv}(k_\parallel )$ flatten out.
4.4.4. Anisotropy scale by scale
To quantify, scale by scale, the anisotropy of the flow, we examine the behaviour of the wavenumber-dependent anisotropy parameters (see, e.g. Sundar et al. Reference Sundar, Verma, Alexakis and Chatterjee2017),
\begin{equation} A_\kappa(k)=\frac{E_{\kappa h}(k)}{2E_{\kappa v}(k)},\quad A_m(k)=\frac{E_{m h}(k)}{2E_{m v}(k)}. \end{equation}
The parameters 
$A_\kappa (k)$ and 
$A_m(k)$ are defined such that 
$A_\kappa = 1$ and 
$A_m=1$ for isotropic MHD flow.
Figure 15(a) displays the plot of 
$A_m(k)$ for the runs B04N2, B04N5 and B04N10. In the case without stratification, 
$A_m> 1$ for 
$k\leqslant 4$ and 
$A_m<1$ for 
$k> 4.$ Similar results are found for 
$A_k(k)$ (not shown). Let us recall that in our numerical simulations the mean magnetic field is not strong enough to produce fully anisotropic dynamics. For the stratified MHD cases, 
$A_m(k)$ is very large at 
$k=1$ and remains greater than one for 
$k > 1.$ The peak at 
$k=1$ for the ratio 
$E_{mh}(k)/(2E_{mv}(k))$ becomes larger and larger as 
$N$ increases. This is caused not by an excess energy of 
$E_{mh}(k)$ but rather due to the fact that 
$E_{mv}(k)$ is very small at almost all the scales. For the velocity field, 
$A_\kappa (k)< 6$ (not shown). However, for 
$k_\perp =0$ (i.e. 
$\boldsymbol {k}\parallel \hat {\boldsymbol {z}})$, both 
$A_\kappa (k_\perp =0)$ and 
$A_\kappa (k_\perp =0)$ are very large. This is due to the fact that, at 
$k_\perp =0$ the vertical components of kinetic and magnetic energies contain much less energy than the horizontal components. On the other hand, the fact that 
$E_{mv}\ll E_{mh}$ almost at all scales does not imply that there exists an inverse cascade of energy, especially since the radial fluxes for
\begin{equation} {\rm \pi}_{mh}(k)=\sum_{p=k}^{k_{max}}{\mathcal{N}}_{nl}^{(mh)}(p),\quad {\rm \pi}_{mv}(k)=\sum_{p=k}^{k_{max}}{\mathcal{N}}_{nl}^{(mv)}(p) \end{equation}
are positive as illustrated by figure 15(b). Here, 
${\mathcal {N}}_{nl}^{(mh}(k)$ and 
${\mathcal {N}}_{nl}^{(mv}(k)$ denote the radial spectra of the nonlinear transfer terms for the horizontal and vertical magnetic energies, respectively. From figure 15(b), we observe that 
${\rm \pi} _{mv}(k)$ is very small compared with 
${\rm \pi} _{mh}(k).$ We note that for Alfvenic turbulence, where there is a forward cascade only, it is observed that 
$E_{mh}\gg E_{mv}$ (see Alexandrova et al. Reference Alexandrova, Carbone, Veltri and Sorriso-Valvo2008; TenBarge et al. Reference TenBarge, Podesta, Klein and Howes2012).
(a) Scale by scale, the ratio between the horizontal magnetic energy and twice the vertical magnetic energy at 
$t=5$ for the runs B04N2, B04N5 and B04N10. The DNS results for the non-stratified MHD case of 
$B_0=0.4$ are also reported for comparison. (b) Radial fluxes of vertical kinetic and magnetic energies at 
$t=5$ for the run B04N5.
4.4.5. Waves and NP mode
In figure 16(a) we show the energy spectra of Alfvén waves, gravity waves and NP mode versus the horizontal wavenumber 
$k_\perp$ for the run B04N5 and 
$t=5.$ The energy spectra 
$E^{(aw)}(k_\perp )$ and 
$E^{(gw)}(k_\perp )$ are very similar, except that the energy of Alfvén waves 
$E^{(aw)}(k_\perp )$ has larger magnitude than the energy of gravity waves 
$E^{(gw)}(k_\perp )$, at all horizontal scales, in agreement with the energy partition shown in figure 16(a). Similar results are found for vertical wavenumber spectra (not shown). The high wavenumber energy of 
$E^{(aw)}(k_\perp )$ (or 
$E^{(gw)}(k_\perp ))$ decays much more rapidly in time than the low wavenumber energy, as shown in figure 16(b). Furthermore, for 
$4 < k_\perp < 15,$ 
$E^{(aw)}(k_\perp )$ tends to follow the power law 
${\sim }k_\perp ^{-2}$ for 
$0< t\leqslant 3,$ while in a large time range, 
$3< t \leqslant 20,$ it behaves like 
$k_\perp ^{-3}$ as shown by figure 16(b).
(a) Horizontal wavenumber spectra of energies of Alfvén waves 
$(E^{(aw)}(k_\perp ))$ gravity waves 
$(E^{(gw)}(k_\perp ))$ and the NP mode 
$(E^{(np)}(k_\perp ))$ at 
$t=5$ for the run B04N5. (b) Horizontal wavenumber spectra 
$E^{(aw)}(k_\perp,t)$ wave energy at several values of time 
$t$ for the run B04N5.
Note that in the (non-magnetized) rotating and stratified case, Kafiabad, Savva & Vanneste (Reference Kafiabad, Savva and Vanneste2019) give an explanation for the transient slope 
$k_\perp ^{-2}$ based on the scattering of inertia-gravity waves by the QG mode. As emphasized here, the analogy between our case and the rotating stratified one is essentially formal, whereas the physical meaning of the NP mode is very different. Especially the role of system rotation is probably crucial in the dynamics of the QG mode, whereas rotation does not affect the MIPS.
On the other hand, the energy of the NP mode remains low, particularly at large horizontal scales. The significant reduction of the energy of NP mode as 
$N$ increases rather occurs at large vertical (or horizontal) scales, that become less an less energetic as 
$N$ increases. This is produced by the reduction of the nonlinear transfer as 
$N$ increases (not shown). At large vertical (or horizontal) scales, 
$E^{(np)}(k_\parallel )$ (or 
$E^{(np)}(k_\perp )$) exhibits a flat shape even at early time.
5. Concluding remarks
Linear theory and DNS were used to study the effects of stable ambient density stratification and a weak imposed mean magnetic field on initially isotropic turbulence, for an electrically conducting Boussinesq fluid with unitary thermal and magnetic Prandtl numbers. In all the simulations the mean magnetic field (
$\boldsymbol {B}_0=B_0\hat {\boldsymbol {z}})$ aligns with the buoyancy gradient of strength 
$N.$
In the inviscid linear limit, a normal mode decomposition has been performed, indicating the existence of three types of motion depending on whether both the Froude (
$F_r)$ and the Alfvén–Mach (
$M)$ numbers are small (case (I)) or only 
$Fr$ is small but 
$M$ is finite (case (II)). In the former case, there is an NP mode, Alfvén waves with frequency 
$\omega _a=B_0k_\parallel$ and magnetogravity waves with frequency 
$\omega _{ag}$ such that 
$\omega _{ag}^2={\omega _a^2+ \omega _g^2}.$ In case (II), there are slow modes (i.e. associated to zero eigenvalues of the fast linear operator 
${\mathcal {L}}_F$ given by (2.19a,b)) and fast gravity waves of pulsation 
$\omega _g.$ In this case, the Alfvén waves are slow modes. Recall that 
$\omega _g=Nk_\perp /k$ is the buoyancy frequency, and 
$k_\parallel$ and 
$k_\perp$ are the vertical and horizontal wavenumbers, respectively.
All the simulations were performed with an improved version of the SNOOPY code and solve numerically the MHD Boussinesq equations for incompressible decaying stratified MHD turbulence under a weak mean magnetic field (see (2.1)). The DNS were started from initial isotropic conditions, with zero initial magnetic and density fluctuations, and unitary magnetic and thermal Prandtl numbers. These initial conditions constitute severe constraining conditions for the generation of the NP mode. It is worthwhile to recall that its definition (see (3.12), (3.19)) excludes the real vortex part, in contrast with the non-magnetized stratified case. To our knowledge, the present numerical simulations are among the first simulations that study the effects of the buoyancy force and the mean magnetic field, when they are simultaneously present, on MHD turbulence. For this, we have opted for a decaying turbulence mainly to avoid any artefact due to the external forcing (see also Meyrand et al. Reference Meyrand, Galtier and Kiyani2016). Eleven runs, at a moderate initial Reynolds number (see table 1), including two runs for the non-stratified MHD cases of 
$B_0=0.2$ and 
$B_0=0.4$ and three runs for PST with 
$N=2,$ 
$N=5$ and 
$N=10$ have been performed.
The development of 
$F_{rh}$ (horizontal Froude number), 
$\tau _{nl}/\tau _a$ (ratio of nonlinear time to Alfvén time) and 
$R_{eb}$ (buoyancy Reynolds number) indicates that the regime of the stratified MHD cases studied in the present work rather corresponds to a weak stratified MHD turbulence regime affected by viscosity.
We studied the temporal behaviour of several global quantities to characterize the dynamics of stratified MHD flows and the influence of the intensity of the buoyancy force and the mean magnetic field on it. While for the non-stratified MHD cases, the total energy 
$E(t)$ tends to follow the power law decay 
${\sim }t^{-1.25},$ for the stratified MHD cases as well as for the purely stratified cases 
$E( t)$ tends to behave like 
${\sim }t^{-0.95}.$ For the same initial Reynolds number, the decay of 
$E(t)$ depends on 
$N$ and 
$B_0$ and not only on the ratio 
$B_0/(L_iN)$ and the energy level increases as 
$N$ increases even though the increase observed between the cases of 
$N=5$ and 
$N=10$ is very slight. We found that vertical motions are more affected by the effect of stratification than by the effect of the mean magnetic field, while it is the opposite for horizontal motions. Indeed, the profiles of the vertical kinetic energy, potential energy and buoyancy flux for the stratified MHD cases and the purely stratified cases look very similar with decaying oscillations whose periods are approximately the same 
${\sim }{\rm \pi} /N.$ The effect of stratification on the development of the vertical magnetic energy is very significant in the sense that, at 
$0< Nt/{\rm \pi} <2,$ there is a drastic decrease of the energy level as 
$N$ increases (see figure 8). In counterpart, the horizontal kinetic and magnetic energies are more impregnated by the mean magnetic field although there is a slight increase of energy level as 
$N$ increases (see figure 10). Compared with the non-stratified MHD cases, stable stratification generates significant anisotropy by restricting vertical motions in favour of horizontal fluid motions. This is particularly reflected in the behaviour of the ratio between horizontal and vertical magnetic energies, which takes a higher value almost at all scales: for 
$N=10$ the peak at 
$k=1$ of the ratio 
$E_{mh}(k)/(2E_{mv}(k))$ exceeds 
$10^3$. This is caused not due to excess energy of 
$E_{mh}$ but rather to the fact that 
$E_{mv}$ is very small in almost all scales. However, the fact that 
$E_{mv}\ll E_{mh}$ at almost all scales does not imply that there is an inverse cascade of energy, especially since radial fluxes for horizontal magnetic energies and verticals have a positive sign.
Concerning the time-development of the energy of Alfvèn waves and magnetogravity waves, it tends to follow the power law decay 
${\sim }t^{-0.95}$ as for the total energy. The contribution coming from Alfvén waves to the total energy is more important than that of magnetogravity waves: at high horizontal (or vertical) wavenumbers the main contribution comes from Alfvén waves. This is mainly due to the fact that, at these scales, 
$\omega _a=B_0\vert k_\parallel \vert$ is more important than 
$\omega _g=(k_\perp /k)N,$ so that, Alfvén waves carry higher energy in that range. On the other hand, it is shown that for 
$4\leqslant k_\perp \leqslant 15,$ the energy spectra of magnetogravity waves and Alfvén waves tend to follow the power law 
${\sim }k_\perp ^{-3}$ for 
$3< t<20.$ At large and intermediate horizontal (or vertical) scales, the energy of NP mode, which is more sensitive on both 
$B_0$ and 
$N$ than the energy of waves, exhibits a flat shape.
Funding
R.M., F.F. and R.R. acknowledge support from the project ‘EVENTFUL’ (ANR-20-CE30-0011), funded by the French ‘Agence Nationale de la Recherche’ – ANR through the program AAPG-2020. The computing resources utilized in this work were provided by PMCS2I at the École Centrale de Lyon.
Declaration of interests
The authors report no conflict of interest.
