# Waves and non-propagating mode in stratified and rotating magnetohydrodynamic turbulence

W. Mouhali, A. Salhi, T. Lehner, C. Cambon

► To cite this version:

W. Mouhali, A. Salhi, T. Lehner, C. Cambon. Waves and non-propagating mode in stratified and rotating magnetohydrodynamic turbulence. Physics of Fluids, 2024, 36, pp.12519. 10.1063/5.0243689. insu-04858036

HAL Id: insu-04858036

https://insu.hal.science/insu-04858036v1

Submitted on 29 Jan 2025

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

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

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

Distributed under a Creative Commons Attribution 4.0 International License
# Waves and non-propagating mode in stratified and rotating magnetohydrodynamic turbulence

W. Mouhali,¹ A. Salhi,² T. Lehner,³ and C. Cambon⁴

¹) *LyRIDS, ECE-Paris Engineering School, 10 Rue Sextius-Michel, 75015 Paris, France.*

²) *Faculté des Sciences de Tunis, Université Tunis El Manar, 2092 El Manar, Tunisia.*ª)

³) *LUTH, UMR 8102 CNRS, Observatoire de Paris-Meudon, 5 place de Janssen, F-92195 Meudon, France.*

⁴) *CNRS, École Centrale de Lyon, INSA Lyon, Université Claude Bernard Lyon 1, Laboratoire de Mécanique des Fluides et d'Acoustique, UMR 5509, F-69134 Ecully, France.*

(Revised 29 January 2025)

In this study we consider a freely decaying, stably stratified and rotating homogeneous magneto-hydrodynamic (MHD) turbulent plasma with a vertical background magnetic field ($B_0 = B_0\hat{z}$), aligned with the density gradient (with a constant Brunt-Väisälä frequency $N$) viewed in a frame rotating uniformly around the vertical axis ($\Omega_0 = \Omega_0\hat{z}$). Quasi-linear theory is used to analyse the flow dynamics for an inviscid and non diffusive Boussinesq fluid. We perform a normal mode decomposition emphasizing three types of motions: a non-propagating (NP) mode, which is no longer a vortex mode, and slow and fast magneto-inertia-gravity waves. The total energy as well as the $L_2$ norm, say $\Gamma$, of the magnetic induction potential scalar (MIPS), which remains similar to the potential enstrophy for non-magnetized rotating and stratified flows, are inviscid invariants. In contrast with the potential vorticity for non-magnetized rotating and stratified flows, the MIPS is not affected by system rotation in the quasi-linear limit, and this is the effect of rotation which presumes an inverse cascade of energy in the equilibrium statistical mechanics. We characterized the system setting up our investigation from the point of view of equilibrium statistical mechanics in the limit of small Froude number and small Alfvén-Mach number. In this limit, the non-propagating quantity $\Gamma$ can be approximated by its quadratic part which explicitly depends only on the vertical component of the fluctuating magnetic field and the density fluctuations. We demonstrate that the partition function restricted to the non-propagating manifold does not indicate an inverse cascade of energy.

## I. INTRODUCTION

Magneto-inertia-gravity (MIG) waves can occur in an electrically conducting fluid submitted to the Lorentz, buoyancy, and Coriolis forces. They are ubiquitously present in several geophysical and astrophysical systems (like the Earth's core, Sun's interior, solar corona, astrophysical accretion disks¹⁻⁶). The MIG waves are used to explain the secular variation of the Earth's magnetic field over the course of hundreds of years¹ and the redistribution of angular momentum in the Sun². They also play an important role in the dynamo theory³⁻⁵.

In addition to the presence of Alfvén waves, acoustic waves and MHD Rossby waves in the solar atmosphere⁷⁻¹⁰, several observations provided strong evidence for the presence of internal gravity waves (IGWs)⁸,¹¹⁻¹³. The IGWs 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¹⁴.

In this paper, we study the interaction of fast and slow MIG waves as well as non-propagating modes occurring

in stably stratified magneto-hydrodynamic (MHD) turbulence subject to a uniform mean magnetic field for an electrically conductive, Boussinesq fluid. The mean (or base) flow is generated by constants $N$ (Brunt-Väisälä frequency), $B_0$ (magnetic field intensity) and $\Omega_0$ (rotation rate), according to the figure 1 where the mean magnetic field $B_0 = B_0\hat{z}$, the rotation vector $\Omega_0 = \Omega_0\hat{z}$ as well as the mean buoyancy scalar gradient $N\hat{z}$ are parallel to the unit vertical vector $\hat{z}$. This simple model allows to study the dynamics of a rotating and stratified (statistically) homogeneous MHD turbulence subjected to a mean magnetic field. It is relevant to characterize turbulent motions occurring in several astrophysical and geophysical systems (e.g., turbulent motions close to the visible solar surface, the Earth's core).

The flow is expressed as a linear combination of modes, each of which is a solution of the equations of motion in the absence of nonlinearity and dissipation. In the non-magnetized case, modes have the form $\exp(i(k \cdot x + s_k \omega(k))$, where $k$ is the wave vector, $s_k = 0, \pm 1$, and $\omega(k) = \sqrt{\omega_r^2 + \omega_g^2}$ is the dispersion relation of inertial-gravity waves, where $\omega_r = 2\Omega|k_\parallel|/k$ and $\omega_g = Nk_\perp/k$ denote the inertial and gravity frequencies, respectively. Here, $k_\parallel = k \cdot \hat{z}$, $k_\perp = \|k \times z\|$ and $k^2 = k_\perp^2 + k_\parallel^2$. The $s_k = \pm 1$ component of the flow represents such waves, whereas the $s_k = 0$ one, here referred to

<footnote>ª) Electronic mail: abdelaziz.salhi@fst.utm.tn</footnote>
<figure><img src="image_2.png" /><figcaption>FIG. 1. A schematic drawing of the mean flow in a frame rotating uniformly about the vertical axis ($\Omega_0 = \Omega_0\hat{z}$) in the presence of a vertical uniform magnetic field ($B = B_0\hat{z}$) and vertical stable stratification ($N$ being the Brunt-Väisälä frequency). The gravity vector is given by $g = -g\hat{z}$ with $g > 0$.</figcaption></figure>

as the NP (non-propagating) component, does not. The NP mode (without magnetic effect) is often described in the geophysical literature as the vortical mode or the potential vorticity (PV) mode¹⁵,

$$ \tilde{\Pi}_{\kappa} = fN + N\omega_z + f\partial_z\vartheta + \boldsymbol{\omega}\cdot\nabla\vartheta, \quad (1) $$

where $f = 2\Omega_0$ being the Coriolis parameter, $\omega_z = \boldsymbol{\omega}\cdot\hat{z}$ is the vertical component of the vorticity vector, $\boldsymbol{\omega} = \nabla \times \boldsymbol{u}$ and $\vartheta$ denotes the fluctuating buoyancy scalar. Note that the constant part $fN$ does not participate in the dynamics and therefore can be ignored. The above modes form a complete set for incompressible flow and so, even in the presence of non-linearity and dissipation, the flow can be expressed as a combination of modes. In purely stratified flows, the wave-vortex decomposition of Riley et al.¹⁶, offered a simple definition of the NP 'vortex' mode in physical space, as horizontal flow (perpendicular to gravity) generated by a streamfunction (as in two-dimensional) turbulence. In contrast with 2D turbulence, however, the stream-function was dependent on the vertical coordinate in order to generate a really three-dimensional vortex mode. In the rotating stratified case, the modal decomposition is easily recovered, with a NP mode combining the abovementioned 'vortex' part and a contribution from buoyancy, whereas gravity waves are replaced by inertia-gravity waves.

For rotating and stratified MHD turbulence there are NP mode as well as fast and slow MIG waves with frequencies $\omega_{f,s}$ ($f$ and $s$ for fast and slow, respectively) such that⁶

$$ \omega_{f,s}^2(\mathbf{k}) = \frac{1}{2} \left( 2\omega_a^2 + \omega_r^2 + \omega_g^2 \pm \sqrt{(\omega_r^2 + \omega_g^2)^2 + 4\omega_r^2\omega_a^2} \right) \quad (2) $$

where $\omega_a = B_0|k_\parallel|$ denotes the frequency of Alfvén waves. Therefore the nomenclature used for non-magnetized stratified and rotating turbulence ought to be radically revised in MHD turbulence, since the Alfvén waves affect all the components of the flow, rendering wavy the components, that were coined vortical in the absence of mean magnetic field. This does not question the technique of modal decomposition (see §3 and Salhi et al.¹⁷) but the NP mode in MHD has no longer a possible

connection with a vortex component. For instance, the NP mode in stratified MHD turbulence is connected with the magnetic induction potential scalar (MIPS), which is a Lagrangian invariant for a non-diffusive, electrically conducting Boussinesq fluid, introduced by Salhi et al.¹⁸ (see also Salhi & Cambon¹⁹),

$$ \tilde{\Pi} = B_0 N + N b_z + B_0 \partial_z \vartheta + \boldsymbol{b} \cdot \nabla \vartheta, \quad (3) $$

where $b_z$ is the vertical component of the fluctuating magnetic field. Note that MIPS does not explicitly depend on rotation. The only analogy with the vortex mode in (non-magnetized) stratified turbulence is that the vertical component of the fluctuating magnetic field $b_z$ replaces the vertical vorticity component $\omega_z$ for the MIPS.

It appears useful to briefly recall some previous studies characterizing the role played by $\Pi_\kappa$ in the dynamics of geophysical flows and stratified turbulence. For strong rotation and strong stratification, so that $Ro = U/(fL) \to 0$ and $Fr = U/(NL) \to 0$ where $Ro$ and $Fr$ are the Rossby and Froude numbers, the inertia-gravity are eliminated in the leading-order solution of the rotating-stratified Boussinesq equations giving leading-order quasi-geostrophic (QG) flow (see, e.g.,²⁰,²¹). Here, $U$ and $L$ are characteristic velocity and length scales, respectively. It is well known that turbulent QG flow produces an inverse energy cascade as a result (see, e.g.,²²,²³) with total (kinetic+potential) energy spectrum scaling as $E(k) \propto k^{-5/3}$ for the low wavenumbers, and a forward cascade of potential entrophy, analogous to the enstrophy cascade in two-dimensional turbulence, with $E(k) \propto k^{-3}$ at high wavenumbers²⁴⁻²⁶. 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 of $N/f$ develop inverse cascade²⁷⁻³⁰. Note that in the recent study by Scott & Cambon³¹, weak turbulence theory was applied to rotating stratified turbulence, with the NP mode, that is not possible in principle in conventional wave-turbulence theory. The turbulence was supposed weak, i.e. one or other (or both) of the Rossby and Froude numbers is small, so nonlinearity is negligible over time spans comparable with the periods of inertial gravity waves, but can be effective over longer spans, provided visco-diffusive dissipation does not kill the turbulence beforehand. It was shown that, for weak turbulence and at leading order, the NP component evolves independently of the wave one and can be described numerically using a specifically adapted DNS scheme. On the other hand, the wave spectrum evolves according to wave-turbulence theory.

The inviscid statistical tendency at low $Ro$ and $Fr$ has been explored theoretically and numerically in some previous studies. In the limit $Fr \to 0$ and/or $Ro \to 0$, Bartello³² employed statistical mechanics to consider the effects of combined conservation of total (kinetic+potential) energy and potential entrophy (i.e., the $L_2$ norm of PV). Indeed, in the limit $Fr \to 0$ and/or
$Ro \to 0$, the potential entrophy can be approximated by its quadratic part, and hence, the two inviscid invariants can be expressed in terms of normal modes that are the non-propagating (or slow) mode and the fast modes. In the study by Herbert et al.³³, who also considered the statistical mechanical equilibrium, it was shown that restricting the equilibrium probability distribution to the slow manifold produces an inverse cascade. In contrast, when taking into account the whole phase space including the waves results in a direct cascade.

We propose here to apply in the magnetized case the quasi-linear approach of Herbert et al.³³, using statistical physics and absolute equilibrium. The analysis involves only the NP mode, that is the MIPS mode in the magnetized case instead of the QG mode. The effect of system rotation on the QG mode yielded the presumption of an inverse cascade, in connection with a negative temperature. Since the MIPS is not affected by rotation, at least in the quasi-linear limit, we can exclude the occurrence of an inverse cascade in that limit. Even if system rotation is taken into account, the normal mode analysis shows that the distinction between vortex mode and wave modes is drastically modified in the magnetized case: Alfvén waves affect all the components of the flow and render wavy some modes, which were considered as vortex without external mean magnetic field. Note that in the recent DNS study by Salhi et al.¹⁷ of (non-rotating) stratified MHD turbulence at low $Fr$ and finite Alfvén mach number $M = U/B_0$, it is shown that the $L_2$ norm of MIPS can be approximated by its quadratic part which dominates the cubic and quartic parts allowing for the use of equilibrium statistical mechanics.

The paper is organised as follows. Section 2 describes the equations for the fluctuating MHD flow (velocity, magnetic field and buoyancy scalar) subjected to Coriolis force, buoyancy force, and Lorentz force, in the presence of the canonic flow in figure 1. The normal mode decomposition is performed in section 3. For this purpose the method of Salhi et al.¹⁷ is streamlined, with the use of the orthonormal Craya-Herring frame from the very beginning, as in several publications from our team (see, e., Salhi & Cambon³⁴). A special care is given to the complex way of deriving the case at vanishing magnetic field, with some singular definition of the NP mode. Equilibrium statistical mechanics is investigated in section 4. Conclusions and perspectives are given in section 5.

## II. MATHEMATICAL FORMULATION

### A. The Boussinesq ideal MHD equations

We consider a freely decaying, rotating and stratified homogeneous magneto-hydrodynamic (MHD) turbulence with a background magnetic field in a frame rotating uniformly about a vertical fixed axis. Both the background magnetic field and the gradient of the background density align with the vertical axis. In that frame, the Boussi-

nesq ideal MHD equations can be written as (see e.g., Wilczyński et al.⁹),

$$ (\partial_t + \mathbf{u} \cdot \nabla) \mathbf{u} = -\nabla p - 2\Omega_0 \hat{z} \times \mathbf{u} + ((B_0 \hat{z} + \mathbf{b}) \cdot \nabla) \mathbf{b} + N\vartheta \hat{z}, \quad (4a) $$

$$ (\partial_t + \mathbf{u} \cdot \nabla) \mathbf{b} = ((B_0 \hat{z} + \mathbf{b}) \cdot \nabla) \mathbf{u}, \quad (4b) $$

$$ (\partial_t + \mathbf{u} \cdot \nabla) \vartheta = -N \hat{z} \cdot \mathbf{u}, \quad (4c) $$

$$ \nabla \cdot \mathbf{u} = 0, \quad \nabla \cdot \mathbf{b} = 0, \quad (4d) $$

where $\mathbf{u}$ and $\mathbf{b}$ respectively denote the fluctuating velocity and magnetic field. The total magnetic field $B_0 \hat{z} + \mathbf{b}$, where $B_0$ is a positive constant, is scaled using Alfvén velocity units: it is divided by $\sqrt{\rho_0 \mu_0}$ where $\rho_0$ and $\mu_0$ are the constant density and the magnetic permeability of the fluid and $\hat{z}$ denotes a unit vertical vector such that $(\hat{x}, \hat{y}, \hat{z})$ is the canonical basis of $\mathbb{R}^3$. Here, $\Omega_0$ is a positive constant that denotes the rotation rate and $p$ is the total fluctuating pressure (including the centrifugal potential and the magnetic part) divided by the constant density $\rho_0$.

Indicate that the Boussinesq approximation ignores density differences except where they appear in terms multiplied by $g$ (the acceleration due to gravity, see, e.g., Spiegel & Veronis³⁵). Accordingly, the continuity equation reduces to $\nabla \cdot \mathbf{u} = 0$ as for an incompressible fluid. The instantaneous density $\tilde{\rho}$ is split in the Boussinesq approximation as

$$ \tilde{\rho}(\mathbf{x}, t) = \rho_0 + \varrho(z) + \rho(\mathbf{x}, t) \quad (5) $$

where $\varrho(z)$ is the variation in the absence of motion; and $\rho(\mathbf{x}, t)$ is the fluctuation resulting from motion³⁵. For a (statistically) stratified homogeneous turbulence, a linear profile for $\varrho(z)$ is considered, so that, the (squared) Brunt-Väisälä frequency

$$ N^2 = -\frac{g}{\rho_0} \frac{d\varrho}{dz} \ge 0, \quad \vartheta = -\frac{g}{N\rho_0}\rho, \quad (6) $$

is a positive constant (stable stratification). Here, the buoyancy scalar $\vartheta(\mathbf{x}, t)$ is rescaled as a velocity.

The first equation in the above system is the momentum equation under the Boussinesq approximation, the second is the induction equation for the magnetic field and the third is the buoyancy scalar equation. The RHS of equation (4) represents the material (or Lagrangian) derivative where the term $\partial_t \mathbf{u}$ (respectively, $\partial_t \mathbf{b}$ and $\partial_t \vartheta$) reflects the change in the fluctuating velocity (respectively, magnetic field and buoyancy scalar), whereas the term $(\mathbf{u} \cdot \nabla) \mathbf{u}$ (respectively, $(\mathbf{u} \cdot \nabla) \mathbf{b}$ and $(\mathbf{u} \cdot \nabla) \vartheta$) reflects the change in velocity (respectively, magnetic field and buoyancy scalar) in space, as the fluid particle moves one point to another in the flow field in infinitesimal time. In the momentum equation, the terms $-2\Omega_0 \hat{z} \times \mathbf{u}$ and $N\vartheta \hat{z}$ represent the Coriolis and buoyancy forces per unit mass. The term $((B_0 \hat{z} + \mathbf{b}) \cdot \nabla) \mathbf{b}$ represents the magnetic tension per unit mass. Indeed, the Lorentz force can be split into two parts: one part which corresponds to the
magnetic pressure force and the other part which corresponds to the magnetic tension. For a perfect electrically conducting fluid moving at the velocity $\mathbf{u}$, one deduces from the Maxwell-Faraday equation the induction equation (4b). Note that the magnetic field is always divergence free, i.e., $\nabla \cdot \mathbf{b} = 0$. The buoyancy scalar equation (4c) results from the fact that, in the Boussinesq approximation, the material derivative of the internal energy is zero for a non diffusive fluid.

Note that the pressure $p$ can be eliminated by computing the divergence of the momentum equation in (4) together with $\nabla \cdot \mathbf{u} = 0$ and $\nabla \cdot \mathbf{b} = 0$, and then solving the resulting elliptic equation for the pressure $p$ to get

$$p(\mathbf{x}, t) = -\nabla^{-2} (\nabla \cdot (\mathbf{u} \cdot \nabla \mathbf{u} - \mathbf{b} \cdot \nabla \mathbf{b})) + \nabla \cdot (2\Omega_0 \hat{\mathbf{z}} \times \mathbf{u}) - N \nabla \vartheta \cdot \hat{\mathbf{z}}. \quad (7)$$

The equation for the absolute vorticity vector, $w_a = \nabla \times \mathbf{u} + 2\Omega_0 \hat{\mathbf{z}}$, is obtained by taking the curl of the momentum balance and has the form

$$\partial_t w_a = \nabla \times (\mathbf{u} \times w_a + \mathbf{j} \times \mathbf{b}) + \nabla \vartheta \times \hat{\mathbf{z}} \quad (8)$$

In addition, we report here the equation for the gradient of the buoyancy scalar,

$$(\partial_t + \mathbf{u} \cdot \nabla) (\nabla \vartheta) = -(\nabla \mathbf{u})^T \cdot (\nabla \vartheta) \quad (9)$$

which allows us to deduce the equations for the potential vorticity (PV, see Pedlowsky [15]) and the equation for the so-called “magnetic induction potential scalar” (MIPS, see Salhi et al.¹⁸). Here, $T$ denotes transpose.

## B. Conserved quantities for inviscid system

The inviscid quadratic invariant of the system (4) is the total (kinetic+magnetic+potential) energy,

$$E(t) = \frac{1}{2} \int_{\mathcal{D}} (\|\mathbf{u}(\mathbf{x}, t)\|^2 + \|\mathbf{b}(\mathbf{x}, t)\|^2 + \vartheta(\mathbf{x}, t)^2) \, d^3\mathbf{x}, \quad (10a)$$

$$\partial_t E(t) = 0 \quad (10b)$$

the domain $\mathcal{D}$ corresponds to a triple-periodic domain (of side length $L_0$, taking $L_0 = 2\pi$) and $\mathbf{j}(\mathbf{x}, t)$ denotes the normalized current density. In contrast, both the cross-helicity, $H^{(c)}$, and the magnetic helicity, $H^{(m)}$,

$$H^{(c)} = \frac{1}{2} \int_{\mathcal{D}} (\mathbf{u} \cdot \mathbf{b}) \, d^3\mathbf{x}, \quad H^{(m)} = \frac{1}{2} \int_{\mathcal{D}} (\mathbf{a} \cdot \mathbf{b}) \, d^3\mathbf{x}, \quad (11a)$$

$$\partial_t H^{(c)} = \Omega_0 \hat{\mathbf{z}} \cdot \int_{\mathcal{D}} (\mathbf{b} \times \mathbf{u}) \, d^3\mathbf{x} + N \hat{\mathbf{z}} \cdot \int_{\mathcal{D}} (\vartheta \mathbf{b}) \, d^3\mathbf{x}, \quad (11b)$$

$$\partial_t H^{(m)} = B_0 \hat{\mathbf{z}} \cdot \int_{\mathcal{D}} (\mathbf{b} \times \mathbf{u}) \, d^3\mathbf{x} \quad (11c)$$

are not conserved. Here, $\mathbf{a}$ is the potential vector, $\mathbf{b} = \nabla \times \mathbf{a}$, which satisfies the Coulomb gauge, $\nabla \cdot \mathbf{a} = 0$. We note that in the non-stratified MHD case, ($N=0$), the hybrid helicity, $H = H^{(c)} - (\Omega_0/B_0)H^{(m)}$ is an inviscid

quadratic invariant (provided that the rotation vector $\Omega_0$ and the background magnetic field $\mathbf{B}_0$) are collinear, see Shebalin³⁶; Galtier³⁷).

On the other hand, from Eq. (8) and (9) we deduce the equation for the potential vorticity (PV) which is very useful as an invariant in stratified geophysical flows (see, e.g., Pedlowsky¹⁵),

$$\tilde{\Pi}_{\kappa} = w_a \cdot (N \hat{\mathbf{z}} + \nabla \vartheta) = (2\Omega_0 \hat{\mathbf{z}} + \nabla \times \mathbf{u}) \cdot (N \hat{\mathbf{z}} + \nabla \vartheta), \quad (12a)$$

$$(\partial_t + \mathbf{u} \cdot \nabla) \tilde{\Pi}_{\kappa} = (N \hat{\mathbf{z}} + \nabla \vartheta) \cdot (\nabla \times (\mathbf{j} \times \mathbf{b})) \quad (12b)$$

It follows that, in the case of a magnetized Boussinesq ideal fluid, $\tilde{\Pi}_{\kappa}$ is not a Lagrangian invariant 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, Salhi et al.¹⁸),

$$\tilde{\Pi} = (B_0 \hat{\mathbf{z}} + \mathbf{b}) \cdot (N \hat{\mathbf{z}} + \nabla \vartheta), \quad (13a)$$

$$(\partial_t + \mathbf{u} \cdot \nabla) \tilde{\Pi} = 0 \quad (13b)$$

is conserved following the flow in the absence of magnetic and thermal diffusion. We note that $\tilde{\Pi}$ does not depend explicitly on the velocity field or rotation. We also note that, in the inviscid limit, the untruncated Boussinesq MHD equations are further constrained by the Lagrangian conservation of MIPS, which generates an infinite number of invariants, in particular the $L_2$ norm of MIPS,

$$(\partial_t + \mathbf{u} \cdot \nabla) \left( \int_{\mathcal{D}} \left( \tilde{\Pi} \right)^2 d^3\mathbf{x} \right) = 0, \quad (14)$$

which is analogous to the potential entrophy (i.e., the $L_2$ norm of PV) in the non-magnetised stratified case⁶. Both PV and MIPS are non-propagating quantities.

## C. Non-dimensionalised 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 hor-
izontal and the vertical directions and U is the velocity
scale. As in the study of the case of non-magnetised ro-
tating and stratified system by Embid & Majda²⁰, a non-
dimensionalised version of the Boussinesq MHD equa-

tions can be written in abstract form as

$$
\partial_t' \mathbf{Y} + \mathcal{L}(\mathbf{Y}) + \mathcal{B}(\mathbf{Y}, \mathbf{Y}) = \mathbf{0} \quad (15)
$$

where $t' = tU/L$ denotes the dimensionless time, $\mathbf{Y}$ is given by

$$
\mathbf{Y} = \begin{pmatrix} \mathbf{u}' \\ \mathbf{b}' \\ \vartheta' \end{pmatrix} = U^{-1} \begin{pmatrix} \mathbf{u} \\ \mathbf{b} \\ \vartheta \end{pmatrix}, \quad (16)
$$

the linear operator is given by

$$
\mathcal{L}(\mathbf{Y}) = \begin{pmatrix} -Fr^{-1}\vartheta'\hat{z} - M^{-1}\partial_{x_3'}\mathbf{b}' - Fr^{-1}\nabla' (\nabla'^{-2}\partial_{x_3'}\vartheta') + Ro^{-1}\nabla' (\nabla'^{-2}(\nabla' \cdot (\hat{z} \times \mathbf{u}')) \\ -M^{-1}\partial_{x_3'}\mathbf{u}' \\ Fr^{-1}\mathbf{u}'\cdot\hat{z} \end{pmatrix}, \quad (17)
$$

the quadratic operator is given by

$$
\mathcal{B}(\mathbf{Y}, \mathbf{Y}) = \begin{pmatrix} \mathbf{u}' \cdot \nabla' \mathbf{u}' - \mathbf{b}' \cdot \nabla' \mathbf{b}' - \nabla' (\nabla'^{-2} (\nabla' \cdot (\mathbf{u}' \cdot \nabla' \mathbf{u}' - \mathbf{b}' \cdot \nabla' \mathbf{b}')) \\ \mathbf{u}' \cdot \nabla' \mathbf{b}' - \mathbf{b}' \cdot \nabla' \mathbf{u}' \\ \mathbf{u}' \cdot \nabla' \vartheta' \end{pmatrix}. \quad (18)
$$

Here, $\nabla' = L^{-1}\nabla$ and the non-dimensional numbers

$$
Fr = \frac{U}{NL}, \quad M^{-1} = \frac{B_0}{U}, \quad Ro = \frac{U}{2\Omega_0 L} \quad (19)
$$

are the Froude, inverse Alfvén-Mach and Rossby num-
bers, respectively.

The constant part $NB_0$ of MIPS does not participate
in the dynamics and, therefore, can be ignored. Accord-
ingly, in a non-dimensional form, the expression of MIPS
reduces to

$$
\Pi' = \frac{\Pi}{NU} = b_3' + FrM^{-1}\partial_{x_3'}\vartheta' + Fr\mathbf{b}'\cdot\nabla'\vartheta', \quad (20)
$$

which does not explicitly depend on the Rossby number.

We consider the case of small Froude number, small Alfvén-Mach number and small Rossby number and use the scaling

$$
Fr = \epsilon, \quad M = C_m^{-1}\epsilon, \quad Ro = C_r^{-1}\epsilon, \quad \text{as } \epsilon \to 0, \quad (21)
$$

where $C_m$ and $C_r$ are the proportionality coefficients.
Accordingly, the equation (15) can be rewritten as

$$
\partial_{t'} \mathbf{Y} + \epsilon^{-1} \mathcal{L}_F(\mathbf{Y}) + \mathcal{B}(\mathbf{Y}, \mathbf{Y}) = \mathbf{0}, \quad (22a)
$$

$$
\mathbf{Y}(t=0) = \mathbf{Y}(\mathbf{x}) \tag{22b}
$$

where the linear operator $\mathcal{L} = \epsilon^{-1}\mathcal{L}_F$ has only a fast
component $\mathcal{L}_F$, while the expression of $\Pi'$ is of the form

$$
\Pi' = b_3' + C_m \partial_{x_3'} \vartheta' + \epsilon \mathbf{b}' \cdot \nabla' \vartheta'. \quad (23)
$$

Equation (20) clearly indicates that, in the asymptotic
limit of small Froude number and small Alfvén-Mach
number, the $L_2$ norm of MIPS can be approximated by
its quadratic part $\Gamma_q$,

$$
\Gamma(t) = \frac{1}{2} \int_D (\tilde{\Pi})^2 d^3 x = \Gamma_q + (\Gamma - \Gamma_q) \approx \Gamma_q, \quad (24a)
$$

$$
\Gamma_q(t) = \frac{1}{2} \int_D (B_0 \partial_z \vartheta + Nb_3)^2 d^3 x. \quad (24b)
$$

Therefore, we can consider that, the two inviscid ($\nu =$
$\eta = \kappa = 0$) quadratic invariants of the above system are
the total energy $E$, which is valid uniformly in $\epsilon$, and the
quadratic part $\Gamma_q$ of the $L_2$ norm of MIPS, allowing for
the use of statistical mechanics arguments in this limit
(see §4).

We note that in some previous studies of non-
magnetised rotating and stratified systems with Fr ≪ 1
and Ro ≪ 1, the potential enstrophy Γ(κ) is approxi-
mated by its quadratic part (see, e.g., Bartello³², Waite
& Bartello³⁸, Herbert et al.³³),

$$
\Gamma^{(\kappa)}(t) = \frac{1}{2} \int_{\mathcal{D}} \left( \tilde{\Pi}_{\kappa} \right)^2 d^3 \mathbf{x} \approx \Gamma_{q}^{(\kappa)}, \quad (25a)
$$

$$
\Gamma_q^{(\kappa)}(t) = \frac{1}{2} \int_D (2\Omega_0 \partial_z \vartheta + N(\partial_x u_2 - \partial_y u_1))^2 d^3 \mathbf{x} \quad (25\text{b})
$$

where $u_1 = \hat{\mathbf{x}} \cdot \mathbf{u}$ and $u_2 = \hat{\mathbf{y}} \cdot \mathbf{u}$, so that, the inviscid invariants of the above systems are the total (ki-
netic+potential) energy and the quadratic part of the potential entrophy.

### III. A NORMAL MODE DECOMPOSITION

The normal modes form a complete basis for the whole fluctuating field. They are built from the latter set of variables, resulting in a non-propagating mode (one component) and slow and rapid waves (two-component each). Accordingly, in the linear limit, wavy motion and non-wavy motion are disentangled, with a minimal number of variables. Finally, the amplitudes of the normal modes are decoupled and time-independent in the linear limit, but they can describe a fully nonlinear problem, if the time-dependency is restored, as used in wave turbulence theory and in direct numerical simulations.

We consider the linear equation associated with (22),

$$ \partial_t \mathbf{Y} + \varepsilon^{-1} \mathcal{L}_F(\mathbf{Y}) = \mathbf{0}, \quad (26) $$

and introduce the Fourier decomposition

$$ \mathbf{Y}(\mathbf{x}, t) = \sum_{\mathbf{k} \in \mathcal{T}} U \hat{\mathbf{Y}}(\mathbf{k}, t) e^{i \mathbf{x} \cdot \mathbf{k}} \quad (27) $$

where $\hat{\mathcal{T}} = \mathbb{Z}^3$ is the set of wave vectors and $i^2 = -1$. The condition for the velocity, buoyancy and magnetic fields to be real is that $\tilde{\hat{u}}(\mathbf{k}) = \hat{u}(-\mathbf{k})$, $\tilde{\hat{\vartheta}}(\mathbf{k}) = \hat{\vartheta}(-\mathbf{k})$ and $\tilde{\hat{b}}(\mathbf{k}) = \hat{b}(-\mathbf{k})$, where $\hat{\cdot}$ denotes complex conjugate. In Fourier space, the incompressibility constraint and the divergence-free condition for the magnetic field read $\mathbf{k} \cdot \hat{\mathbf{u}} = 0$ and $\mathbf{k} \cdot \hat{\mathbf{b}} = 0$, respectively, which signifies that both $\hat{\mathbf{u}}$ and $\hat{\mathbf{b}}$ belong to the plane perpendicular to the wave vector $\mathbf{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 3D Fourier space, in the plane normal to $\mathbf{k}$. The pressure is removed since the Poisson equation (see equation (7)) is equivalent to an algebraic relation along $\mathbf{k}$. Accordingly, in the three-dimensional Fourier space, only five components remain, two for the velocity, two for the magnetic field, and one for the buoyancy scalar.

It is worth mentioning that phase information is lost in some spectra of second-order correlations, like the energy, because of the product of a spectral component by its complex conjugate, but Fourier decomposition is relevant in homogeneous turbulence, far from physical boundaries. On the other hand, phase information is essential for wavy components, via the time-dependent, rapidly oscillating terms, in the phase $\mathbf{k} \cdot \mathbf{x} - \omega(\mathbf{k})t$. Fourier decomposition and related normal mode decomposition are very relevant in our study for disentangling wavy motion, with rapid phases, and non-propagating motion. Of

course, the modulation is seen in homogeneous turbulence, so that we agree to quote and discuss the reference suggested by the Referee for an alternative approach$^{39}$

#### A. Differential system for the toroidal and poloidal components

For a non zero horizontal wavenumber, $k_{\perp} \equiv \|\hat{\mathbf{z}} \times \mathbf{k}\| \neq 0$, one can use a local orthonormal frame (i.e., the Cray-Herring frame$^{34}$)

$$ e^{(3)} = k^{-1}\mathbf{k}, \quad e^{(1)} = k_{\perp}^{-1}\mathbf{k} \times \hat{\mathbf{z}}, \quad e^{(2)} = e^{(3)} \times e^{(1)} \quad (28) $$

in which the two conditions $\mathbf{k} \cdot \hat{\mathbf{u}} = 0$ and $\mathbf{k} \cdot \hat{\mathbf{b}} = 0$, are satisfied by construction and $\hat{\mathbf{u}}$ (or $\hat{\mathbf{b}}$) has only two components: the toroidal component (along $e^{(1)}$) and the poloidal component (along $e^{(2)}$),

$$ \hat{\mathbf{u}} = u^{(1)}e^{(1)} + u^{(2)}e^{(2)}, \quad \hat{\mathbf{b}} = u^{(3)}e^{(1)} + u^{(4)}e^{(2)}, \quad (29a) $$

$$ u^{(1)} = k_{\perp}^{-1} (k_2 \hat{u}_1 - k_1 \hat{u}_2), \quad u^{(2)} = -k_{\perp}^{-1} k \hat{u}_3, \quad (29b) $$

$$ u^{(3)} = k_{\perp}^{-1} (k_2 \hat{b}_1 - k_1 \hat{b}_2), \quad u^{(4)} = -k_{\perp}^{-1} k \hat{b}_3 \quad (29c) $$

where $k = \|\mathbf{k}\|$ is the radial wavenumber.

For $k_{\perp} = 0$, so that $k_1 \equiv \hat{\mathbf{x}} \cdot \mathbf{k} = 0$ and $k_2 \equiv \hat{\mathbf{y}} \cdot \mathbf{k} = 0$, the vertical components of $\hat{\mathbf{u}}$ and $\hat{\mathbf{b}}$ vanish because $\mathbf{k} \cdot \hat{\mathbf{u}} = \mathbf{k} \cdot \hat{\mathbf{b}} = 0$. To facilitate calculations in that case, we set

$$ e^{(1)} = \frac{1}{\sqrt{2}} (\hat{\mathbf{x}} - \hat{\mathbf{y}}), \quad e^{(2)} = \frac{1}{\sqrt{2}} (\hat{\mathbf{x}} + \hat{\mathbf{y}}), \quad e^{(3)} = \hat{\mathbf{z}}, \quad (30) $$

and hence,

$$ \hat{\mathbf{u}} = u^{(1)}e^{(1)} + u^{(2)}e^{(2)}, \quad (31a) $$

$$ u^{(1)} = \frac{1}{\sqrt{2}} (\hat{u}_1 - \hat{u}_2), \quad u^{(2)} = \frac{1}{\sqrt{2}} (\hat{u}_1 + \hat{u}_2), \quad (31b) $$

where the same notation is used for the magnetic field $\mathbf{b}$.

Given the equation (17) the substitution of (27) into (26) leads to the linear differential system

$$ d_t \hat{\mathbf{Y}} = \hat{\mathcal{L}} \cdot \hat{\mathbf{Y}} \quad (32) $$

or equivalently,

$$ d_t \begin{pmatrix} u^{(1)} \\ u^{(2)} \\ b^{(1)} \\ b^{(2)} \\ \hat{\vartheta} \end{pmatrix} = \begin{pmatrix} 0 & s\omega_r & is\omega_a & 0 & 0 \\ -s\omega_r & 0 & 0 & is\omega_a & -\omega_g \\ is\omega_a & 0 & 0 & 0 & 0 \\ 0 & is\omega_a & 0 & 0 & 0 \\ 0 & \omega_g & 0 & 0 & 0 \end{pmatrix} \cdot \begin{pmatrix} u^{(1)} \\ u^{(2)} \\ b^{(1)} \\ b^{(2)} \\ \hat{\vartheta} \end{pmatrix} \quad (33) $$

where $\omega_r$, $\omega_a$ and $\omega_g$ are the pulsations of inertial, Alfvén and gravity waves, respectively,

$$ \omega_r = 2\Omega_0 k^{-1} |k_{\parallel}|, \quad \omega_a = B_0 |k_{\parallel}|, \quad \omega_g = N k^{-1} k_{\perp}, \quad (34) $$
and

$$s = k_{\parallel}^{-1} |k_{\parallel}| = \pm 1, \quad k_{\parallel} = k \cdot \hat{z}. \qquad (35)$$

For these waves, the Coriolis force, magnetic tension and buoyancy respectively play the role of resorting force.

The matrix $\hat{\mathcal{L}}(\mathbf{k})$ is a normal matrix: $\hat{\mathcal{L}} \cdot \hat{\mathcal{L}}^{\dagger} = \hat{\mathcal{L}}^{\dagger} \cdot \hat{\mathcal{L}}$, and is anti-hermitian: $\hat{\mathcal{L}}^{\dagger} = -\hat{\mathcal{L}}$ (while $i\hat{\mathcal{L}}$ is hermitian),

where $\dagger$ denotes the transpose and complex conjugation, $(\hat{\mathcal{L}}^{\dagger})_{lm} = (\tilde{\mathcal{L}})_{ml}$. Then, there exists a unitary matrix $\mathbf{P} \in \mathcal{U}(5)$ (i.e., $\mathbf{P} \cdot \mathbf{P}^{\dagger} = \mathbf{P}^{\dagger} \cdot \mathbf{P} = \mathbf{I}_5$) such that

$$\hat{\mathcal{L}} = \mathbf{P} \cdot \mathrm{diag}(\lambda_j) \cdot \mathbf{P}^{\dagger} \qquad (36)$$

where $\lambda_j = -i\omega_j$ ($j = 1, 2, \dots, 5$) are the eigenvalues of the matrix $\hat{\mathcal{L}}$,

$$\omega_{np} \equiv i\lambda_1 = \omega_1 = 0, \qquad (37a)$$

$$\omega_s \equiv i\lambda_2 = \omega_2 = -i\lambda_3 = -\omega_3 = \sqrt{\frac{1}{2}\left(2\omega_a^2 + \omega_r^2 + \omega_g^2 - \sqrt{(\omega_r^2 + \omega_g^2)^2 + 4\omega_a^2\omega_r^2}\right)}, \qquad (37b)$$

$$\omega_f \equiv i\lambda_4 = \omega_4 = -i\lambda_5 = -\omega_5 = \sqrt{\frac{1}{2}\left(2\omega_a^2 + \omega_r^2 + \omega_g^2 + \sqrt{(\omega_r^2 + \omega_g^2)^2 + 4\omega_a^2\omega_r^2}\right)}, \qquad (37c)$$

It clearly appears that, in addition to the non-propagating mode corresponding to $\omega_{np} = 0$, there are slow and fast magneto-inertia-gravity waves with pulsations $\omega_s$ and $\omega_f$, respectively. Likewise, for the non-rotating stratified MHD case, so that $\omega_r(\mathbf{k}) = 0$, there are a non-propagating mode corresponding to $\omega_{np} = 0$, Alfvén waves of pulsation $\omega_s = \omega_a$ and magneto-gravity waves of pulsation $\omega_f = \omega_{ag} = \sqrt{\omega_a^2 + \omega_g^2}$. (see Salhi et al.¹⁷).

values $\lambda_{2,3} = \pm i\omega_s$ and $\mathbf{Z}_f^{\mp}$ corresponding to the eigenvalues $\lambda_{4,5} = \pm i\omega_f$, so that,

$$\hat{\mathcal{L}} \cdot \mathbf{Z}_0 = 0, \quad \hat{\mathcal{L}} \cdot \mathbf{Z}_s^{\mp} = \pm i\omega_s \mathbf{Z}_s^{\mp}, \quad \hat{\mathcal{L}}_0 \cdot \mathbf{Z}_f^{\mp} = \pm i\omega_f \mathbf{Z}_f^{\mp}.$$

The set of these vectors form an orthonormal frame. In the canonical basis of $\mathbb{C}^5$, they are the columns of the matrix $\mathbf{P}$ (see Appendix), so that,

## B. Normal modes

We introduce the eigenvectors $\mathbf{Z}_0(\mathbf{k})$ corresponding to the eigenvalue $\lambda_1 = 0$, $\mathbf{Z}_s^{\mp}(\mathbf{k})$ corresponding to the eigen-

where the normal modes are of the form

$$\hat{\mathbf{Y}} = A_0 \mathbf{Z}_0 + A_s^{\mp} \mathbf{Z}_s^{\mp} + A_f^{\mp} \mathbf{Z}_f^{\mp} \qquad (38)$$

$$A_0 = \hat{\mathbf{Y}} \cdot \mathbf{Z}_0^{\dagger} = \frac{1}{\sqrt{\omega_a^2 + \omega_g^2}} (\omega_g b^{(2)} - i \sigma_a \hat{\vartheta}), \qquad (39a)$$

$$A_s^{\mp} = \hat{\mathbf{Y}} \cdot \mathbf{Z}_s^{\mp\dagger} = \frac{1}{\sqrt{2(\sigma_s^2 + \omega_a^2)}} (\omega_s u^{(1)} \mp i\sigma_s u^{(2)} \mp s\omega_a b^{(1)} + \omega_s^{-1}\sigma_s (i\omega_a b^{(2)} - \omega_g \hat{\vartheta})), \qquad (39b)$$

$$A_f^{\mp} = \hat{\mathbf{Y}} \cdot \mathbf{Z}_f^{\mp\dagger} = \frac{1}{\sqrt{2(\sigma_f^2 + \omega_a^2 + \omega_g^2)}} (-i\sigma_f u^{(1)} \mp \omega_f u^{(2)} \mp \omega_f^{-1}\sigma_f s\omega_a b^{(1)} - s\omega_a b^{(2)} - i\omega_g \hat{\vartheta}), \qquad (39c)$$

with

$$\sigma_s = \frac{\omega_s^2 - \omega_a^2}{s\omega_r}, \quad \sigma_f = \frac{s\omega_r\omega_f^2}{\omega_f^2 - \omega_a^2}, \qquad (40)$$

which reduce to $\sigma_s = \omega_s$ and $\sigma_f = \omega_f$ for the case where

$k_{\perp} = 0$ for which $\omega_g = 0$. Note that the above normal modes provide us with insights on the structure of the phase space (see §4).

For the non-rotating stratified MHD case the normal modes are easily deduced from (39) because $\omega_r(\mathbf{k}) = 0$,
$\omega_s = \omega_a$, $\omega_f = \sqrt{\omega_a^2 + \omega_g^2}$, $\sigma_f = 0$ and $\lim_{\omega_r \to 0} \sigma_s = 0$ (see Appendix). For the non-stratified case ($N=0$), the differential system (33) reduces to a differential system of dimension 4 because $\omega_g = 0$ and $\hat{\vartheta} = 0$. In this flow case, there are only slow and fast magneto-inertia waves of pulsation,

$$ \omega_s = \omega_1 = -\omega_2 = \frac{\omega_r}{2} \left( \sqrt{1 + (B_0 k / \Omega_0)^2} - 1 \right), \quad (41a) $$

$$ \omega_f = \omega_3 = -\omega_4 = \frac{\omega_r}{2} \left( \sqrt{1 + (B_0 k / \Omega_0)^2} + 1 \right) \quad (41b) $$

and of a circular polarisation (see, e.g., Galtier [37]), and not a non-propagating mode, $A_0 = 0$. The non-zero normal modes are obtained by replacing in (39) $\sigma_s$ by $\omega_s$ and $\sigma_f$ by $\omega_f$, with $A_0 = 0$. The deduction of the normal modes for the case of non-magnetized rotating and stratified turbulence is considered in §3.4. In the limit of small Froude number, $Fr \to 0$ and small Alfvén-Mach number, $M \to 0$, $\Gamma(t)$ is approximated by its quadratic part $\Gamma_q(t)$. Accordingly, $E(t)$ and $\Gamma_q(t)$ can be expressed in terms of the above normal modes as follows:

$$ E(t) = E^{(np)}(t) + E^{(sw)}(t) + E^{(fw)}(t), \quad (42a) $$

$$ E^{(np)} = \frac{1}{2} \sum_{\mathbf{k} \in \hat{\tau}} |A_0|^2, \quad (42b) $$

$$ E^{(sw)} = \frac{1}{2} \sum_{\mathbf{k} \in \hat{\tau}} \left( |A_s^-|^2 + |A_s^+|^2 \right), \quad (42c) $$

$$ E^{(fw)} = \frac{1}{2} \sum_{\mathbf{k} \in \hat{\tau}} \left( |A_f^-|^2 + |A_f^+|^2 \right), \quad (42d) $$

$$ \Gamma_q(t) = \frac{1}{2} \sum_{\mathbf{k} \in \hat{\tau}} (\omega_a^2 + \omega_g^2) |A_0|^2, \quad (42e) $$

where $E^{(fw)}$, $E^{(sw)}$ and $E^{(np)}$ are the energy of fast waves, slow waves and non-propagating (or slow) modes. Recall that $\Gamma_q(t)$ does not explicitly depend on the rotation or the fluctuating velocity. Its spectral density is equal to that of the non-propagating mode multiplied by $\omega_{ag}^2 = (\omega_a^2 + \omega_g^2)$ (i.e., the square of the pulsation of magneto-gravity waves).

### C. Case of non-magnetized rotating and stratified turbulence

For the case of non-magnetized rotating and stratified turbulence, the system of linear differential equation (32) reduces to a differential system of dimension 3, because $\omega_a(\mathbf{k}) = 0$ and $\hat{b} = 0$. Accordingly, the eigenvalues given by (37) reduce to

$$ \lambda_1 = 0, \quad \lambda_2 = -\lambda_3 = -i\omega_{ig} = -i\sqrt{\omega_r^2 + \omega_g^2}, \quad (43) $$

and hence, in this flow case, there are a non-propagating (NP) mode (associated to the zero eigenvalues) and

inertia-gravity waves of pulsation $\omega_{ig}$. In that case, the expressions of the normal modes can be found as (see, e.g., Scott & Cambon$^{31}$)

$$ A_{np} = \frac{1}{\sqrt{(\omega_r^2 + \omega_g^2)}} (\omega_g u^{(1)} - s\omega_r \hat{\vartheta}), \quad (44a) $$

$$ A^{\pm} = -\frac{i}{\sqrt{2(\omega_r^2 + \omega_g^2)}} \left( s\omega_r u^{(1)} \pm i\sqrt{\omega_r^2 + \omega_g^2} u^{(2)} + \omega_g \hat{\vartheta} \right). \quad (44b) $$

These expressions can be recovered from equation (39) by setting $\mathbf{b} = 0$, $\omega_a(\mathbf{k}) = 0$ and using the following limits (see Appendix)

$$ \lim_{\omega_a \to 0} \omega_s = 0, \quad (45a) $$

$$ \lim_{\omega_a \to 0} \sigma_s = 0, \quad (45b) $$

$$ \lim_{\omega_a \to 0} \frac{\omega_a^2}{\omega_s^2} = 1 + \frac{\omega_r^2}{\omega_g^2}, \quad (45c) $$

$$ \lim_{\omega_a \to 0} \frac{\omega_s}{\sqrt{\sigma_s^2 + \omega_a^2}} = \frac{\omega_g}{\sqrt{\omega_r^2 + \omega_g^2}}, \quad (45d) $$

$$ \lim_{\omega_a \to 0} \frac{\sigma_s \omega_g}{\omega_s \sqrt{\sigma_s^2 + \omega_a^2}} = \frac{\omega_r}{\sqrt{\omega_r^2 + \omega_g^2}}. \quad (45e) $$

In fact, we find

$$ A_0 = 0, \quad A_s^- = A_s^+ = \frac{1}{\sqrt{2}} A_{np}, \quad A_f^\pm = A^\pm, \quad (46) $$

so that, $|A_{np}|^2 = |A_s^-|^2 + |A_s^+|^2$.

In this flow case, the expression of total (kinetic + potential) energy is of the form

$$ E(t) = E^{(np)} + E^{(iw)}, \quad (47a) $$

$$ E^{(np)} = \frac{1}{2} \sum_{\mathbf{k} \in \hat{\tau}} |A_{np}|^2, \quad (47b) $$

$$ E^{(iw)} = \frac{1}{2} \sum_{\mathbf{k} \in \hat{\tau}} \left( |A^-|^2 + |A^+|^2 \right), \quad (47c) $$

where $E^{(np)}$ and $E^{(iw)}$ are the energy of the non-propagating (or slow) mode and inertia-gravity waves, respectively. On the other hand, the expression of the quadratic part of the potential entrophy is given by

$$ \Gamma_q^{(\kappa)}(\mathbf{k}, t) = \frac{1}{2} \sum_{\mathbf{k} \in \hat{\tau}} k^2 (\omega_r^2 + \omega_g^2) |A_{np}|^2. \quad (48) $$

Recall that, in the limit $Fr \ll 1$ and $Ro \ll 1$, the potential entrophy can be approximated by its quadratic part.

## IV. EQUILIBRIUM STATISTICAL MECHANICS

A direct cascade of energy is dominant in incompressible 3-D isotropic turbulence, and yields a nonlinear energy transfer from largest scales to smallest ones. This
energy is dissipated via viscous effects, so that the decay of energy is significantly larger than in the linear case. On the other hand, the inverse cascade is present in, e.g., 2-D turbulence and amounts to a transfer of energy from smallest to largest scales. Accordingly, the energy is less and less dissipated by viscous effects. An inverse energy cascade is presumed in rotating stratified turbulence without MHD, but not in the present MHD case as shown in this section. Because the theoretical analysis only concerns the non propagating mode, we think that the very different nature of this mode (with a vortex part only in the case without MHD) and the impact of rotation on it (only in the case without MHD) explains the differences.

### A. Restricted partition function

Inviscid unforced truncated stratified MHD turbulence can be studied as a problem of equilibrium statistical mechanics. The approach consists in constructing a proba-

bility density function on the phase space of a dynamical system based only on the invariants of the dynamics and then using it to obtain ensemble-average statistics. This approach has proven useful in the past, for example to account for the forward energy cascade in three-dimensional turbulence and the inverse energy cascade in two-dimensional turbulence (see Holloway$^{40}$).

As reported in §1, Herbert et al. 33 studied the restricted equilibrium and energy cascade in rotating and stratified flows considering that the potential entrophy can be approximated by its quadratic part. The other quadratic invariant being the total energy. They showed that, when the statistics is restricted to the contributions of the slow modes, there exists a regime of negative temperature featuring an infrared divergence in both the radial and the axisymmetric average energy spectrum, characteristic of an inverse cascade regime.

For the MHD Boussinesq flows under study here, we consider that $\Gamma$ (i.e., the $L_2$ norm of MIPS) is approximated by its quadratic part. The other quadratic invariant being the total energy $E$

We define a phase space $\Lambda$ in Fourier space as follows (see, e.g., Herbert et al. 33),

$$ \Lambda = \{ (\hat{u}, \hat{b}, \hat{\vartheta}) \in L_2(\hat{T}) \times L_2(\hat{T}) \times L_2(\hat{T}), \forall \mathbf{k} \in \hat{T}, \hat{Y}(\mathbf{k}) = \bar{\hat{Y}}(-\mathbf{k}), \mathbf{k} \cdot \hat{u} = 0, \mathbf{k} \cdot \hat{b} = 0 \}. \quad (49) $$

Given a truncation to a finite number of wavevectors, $\mathcal{B} = \{\mathbf{k} \in \hat{T}, k_{min} \le k \le k_{max}\}$, the macroconical ensemble probability density function is

$$ f_{\alpha, \beta} (A_0, A_s^{\pm}, A_f^{\pm})_{\mathbf{k} \in \mathcal{B}} = \frac{1}{\mathcal{Z}} \exp(-\beta E - \alpha \Gamma_q) \quad (50) $$

where $\beta$ and $\alpha$ are the Lagrange multipliers associated to conservation of $E$ and $\Gamma_q$, respectively, and $\mathcal{Z}$ is the partition function,

$$ \mathcal{Z} = \prod_{\mathbf{k} \in \mathcal{B}} \int_{\mathbb{R}^{10}} \exp \left\{ -\frac{\Psi}{2} \right\} d\Lambda, \quad (51a) $$

$$ \Psi = (\beta + \alpha \omega_{ag}^2) |C_1|^2 + \beta \sum_{\ell=2}^{5} |C_{\ell}|^2, \quad (51b) $$

$$ d\Lambda = \prod_{\ell=1}^{5} d\Re\{C_{\ell}\} \prod_{\ell=1,5} d\Im\{C_{\ell}\}, \quad (51c) $$

in which $C_1 = A_0$, $C_2 = A_s^-$, $C_3 = A_s^+$, $C_4 = A_f^-$, and $C_5 = A_f^+$. Because the integral (51) can be rewritten as

the product of two factors $\mathcal{Z} = \mathcal{Z}^{(np)}\mathcal{Z}^{(w)}$ where

$$ \mathcal{Z}^{(np)} = \prod_{\mathbf{k} \in \mathcal{B}} \frac{2\pi}{(\beta + \alpha \omega_{ag}^2)}, \quad (52a) $$

$$ \mathcal{Z}^{(w)} = \prod_{\mathbf{k} \in \mathcal{B}} \frac{16\pi^4}{\beta^4}, \quad (52b) $$

(52c)

the factors $\mathcal{Z}^{(np)}$ and $\mathcal{Z}^{(w)}$ can be considered as being the partition functions restricted to the submanifolds $\Lambda^{(np)}(\subset \Lambda)$ and $\Lambda^{(w)}(\subset \Lambda)$, respectively, where $\Lambda^{(np)}$ is defined by the condition

$$ A_s^-(\mathbf{k}) = A_s^+(\mathbf{k}) = A_f^-(\mathbf{k}) = A_f^+(\mathbf{k}) = 0 \quad (53) $$

while the manifold $\Lambda^{(w)}$ is defined the condition $A_0(\mathbf{k}) = 0$, so that $\Lambda = \Lambda^{(np)} \oplus \Lambda^{(w)}$. Thus,

$$ \mathcal{Z} = \mathcal{Z}^{(np)} \mathcal{Z}^{(w)} = \prod_{\mathbf{k} \in \mathcal{B}} \frac{(2\pi)^5}{\beta^4 \left( \beta + \alpha \left( B_0^2 k_{\parallel}^2 + N^2 k^{-2} k_{\perp}^2 \right) \right)} \quad (54) $$

The two restricted partition functions correspond to the restricted macro-conical ensemble probability density functions

$$ f_{\beta, \alpha}^{(np)} = \frac{1}{\mathcal{Z}^{(np)}} \exp(-\beta E - \alpha \Gamma_q) \chi_{\Lambda^{(np)}}, \quad (55a) $$

$$ f_{\beta, \alpha}^{(w)} = \frac{1}{\mathcal{Z}^{(w)}} \exp(-\beta E - \alpha \Gamma_q) \chi_{\Lambda^{(w)}}, \quad (55b) $$
respectively, where $\chi_Y$ is the indicator function defined as follows

$$ \chi_Y(y) = \begin{cases} 1 & \text{if } y \in \mathbf{Y} \\ 0 & \text{otherwise.} \end{cases} \qquad (56) $$

The necessary and sufficient condition for the convergence of the Gaussian integrals is that the restriction of the quadratic form $\beta E + \alpha \Gamma_q$ to $\Lambda^{(np)}$ or $\Lambda^{(w)}$ be positive definite: $\beta \ge 0$ is the condition for $\mathcal{Z}^{(w)}$ to be defined, while the condition for $\mathcal{Z}^{(np)}$ to be defined is that

$$ \forall \mathbf{k} \in \mathcal{B}, \quad \beta + \alpha \omega_{ag}^2 = \beta + \alpha \left( B_0^2 k_{\parallel}^2 + N^2 k^{-2} k_{\perp}^2 \right) \ge 0. \qquad (57) $$

Recall that for the non-magnetized stratified and rotating case, the total (kinetic+potential) energy as well as the quadratic part of the potential entrophy are viscoid invariants in the limit $Fr \to 0$ and $Ro \to 0$ (see (48), (47)). In that case, there are a non-propagating (NP) mode and inertia-gravity waves (igw), as already indicated, and the expression of the partition functions $\mathcal{Z}^{(np)}$ and $\mathcal{Z}^{(igw)}$ given by (51a,b) reduces to

$$ \mathcal{Z}^{(np)} = \prod_{\mathbf{k} \in \mathcal{B}} \frac{2\pi}{(\beta_r + \alpha_r k^2 \omega_{ig}^2)}, \quad \mathcal{Z}^{(w)} = \prod_{\mathbf{k} \in \mathcal{B}} \frac{4\pi^2}{\beta^2}, \qquad (58) $$

where $\alpha_r$ and $\beta_r$ are the Lagrange multipliers associated to the total energy $E(t)$ and the quadratic part of the potential entrophy, $\Gamma_q^{(\kappa)}$ and $\omega_{ig}^2 = \omega_r^2 + \omega_g^2$. The condition for $\mathcal{Z}^{(np)}$ to be defined is that

$$ \forall \mathbf{k} \in \mathcal{B}, \quad \beta_r + k^2 \alpha_r \omega_{ig}^2 = \beta_r + \alpha_r \left( 4\Omega_0 k_{\parallel}^2 + N^2 k_{\perp}^2 \right) \ge 0. \qquad (59) $$

and $\beta_r \ge 0$ is the condition for $\mathcal{Z}^{(w)}$ to be defined. We notice that the realisability condition (59) can be recovered from (57) by replacing the square of Alfvén waves $B_0^2 k_{\parallel}^2$ by the square of inertial waves $4\Omega_0^2 k_{\parallel}^2 k^{-2}$ and the Lagrange multiplier $\alpha$ by $k^2 \alpha_r$. As indicated by Herbert et al. [33], the system that has phase space restricted to non-propagating modes can reach negative temperature states. In that case, they have identified three regimes by analogy with 2D flows (see, e.g., Kraichnan & Montgomery$^{41}$) and 3-D helically constrained flows (see, e.g., Herbert$^{42}$):

$$ \text{regime (I):} \quad \alpha_r > 0 \text{ and } \beta_r < 0, \qquad (60a) $$

$$ \text{regime (II):} \quad \alpha_r > 0 \text{ and } \beta_r > 0, \qquad (60b) $$

$$ \text{regime (III):} \quad \alpha_r < 0 \text{ and } \beta_r > 0, \qquad (60c) $$

see their figure 1. For regimes (I) and (III), the condition (59) implies

$$ \beta_r + \alpha_r k_{max}^2 \max(N^2, 4\Omega_0^2) > 0, \qquad (61a) $$

$$ \beta_r + \alpha_r k_{min}^2 \min(N^2, 4\Omega_0^2) > 0, \qquad (61b) $$

respectively. In the next section, we show that for the stratified MHD case with $B_0 \ne 0$, the Lagrange multiplier $\alpha$ must be positive ($\alpha \ge 0$).

## B. Average value of energy at statistical equilibrium

The average value of energy at statistical equilibrium with respect to the full and restricted canonical probability distributions are given by

$$ \langle E \rangle = -\partial_{\beta} \ln \mathcal{Z} = \sum_{\mathbf{k} \in \mathcal{B}} \left[ \frac{4}{\beta} + \frac{1}{\beta + \alpha \left( B_0^2 k_{\parallel}^2 + N^2 k^{-2} k_{\perp}^2 \right)} \right], \qquad (62) $$

and

$$ \langle E \rangle_{np} = -\partial_{\beta} \ln \mathcal{Z}^{(np)} = \sum_{\mathbf{k} \in \mathcal{B}} \frac{1}{\beta + \alpha \left( B_0^2 k_{\parallel}^2 + N^2 k^{-2} k_{\perp}^2 \right)}, \qquad (63) $$

respectively. We can similarly obtain the average value of $\Gamma_q$ with respect to the restricted canonical probability distribution:

$$ \langle \Gamma_q \rangle_{np} = -\frac{\partial}{\partial \alpha} \ln \mathcal{Z}^{(np)} = \sum_{\mathbf{k} \in \mathcal{B}} \frac{\omega_{ag}^2}{\beta + \alpha \omega_{ag}^2} \qquad (64) $$

Note that $\langle \Gamma_q \rangle = -\partial_{\alpha} \mathcal{Z} = \langle \Gamma_q \rangle_{np}$. For convenience, in the above expressions, we replace discrete sums by integrals, so that,

$$ \langle E \rangle_{np} = \int_{\mathcal{F}} \frac{1}{\beta + \alpha \left( B_0^2 k_{\parallel}^2 + N^2 k^{-2} k_{\perp}^2 \right)} d^3 \mathbf{k}, \qquad (65) $$

where $\mathcal{F} = \{\mathbf{k} \in \mathbb{R}^3 \text{ such that } k_{min} \le k \le k_{max}\}$. For instance, the resulting radial spectrum $\langle E(k) \rangle_{np}$ such that

$$ \langle E \rangle_{np} = \int_{k_{min}}^{k_{max}} \langle E(k) \rangle_{np} dk, \qquad (66a) $$

$$ \langle E(k) \rangle_{np} = 2\pi k^2 \int_0^\pi \frac{\sin \theta}{\beta + \alpha (B_0^2 k^2 \cos^2 \theta + N^2 \sin^2 \theta)} d\theta \qquad (66b) $$

is found as
$$ \langle E(k) \rangle_{np} = \begin{cases} \frac{2\pi k^2}{\sqrt{\alpha(\beta+\alpha N^2)(N^2-B_0^2k^2)}} \ln \left( \frac{1+\sqrt{\alpha(N^2-B_0^2k^2)/(\beta+\alpha N^2)}}{1-\sqrt{\alpha(N^2-B_0^2k^2)/(\beta+\alpha N^2)}} \right) & \text{if } B_0k < N \\ \frac{4\pi k^2}{(\beta+\alpha N^2)} & \text{if } B_0k = N \\ \frac{4\pi k^2}{\sqrt{\alpha(\beta+\alpha N^2)(B_0^2k^2-N^2)}} \tan^{-1} \sqrt{\frac{\alpha(B_0^2k^2-N^2)}{(\beta+\alpha N^2)}} & \text{if } B_0k > N. \end{cases} \quad (67) $$

In couterpaert, the variation of the average value of the radial spectrum of the $L_2$ norm of MIPS at statistical equilibrium restricted to the submanifold of non-propagating modes $\Lambda^{(np)}$ is given by

$$ \langle \Gamma_q(k) \rangle_{np} = \frac{16\pi k^2}{\alpha} - \frac{\beta}{\alpha} \langle E(k) \rangle_{np}. \quad (68) $$

For the non-magnetized stratified and rotating case, the integral form of $\langle E(k) \rangle_{np}$ is obtained by replacing in

(66b) the term $(\beta + \alpha\omega_{ag}^2)$ by $(\beta_r + \alpha_r k^2 \omega_{ig}^2)$, so that,

$$ \langle E(k) \rangle_{np} = 2\pi k^2 \int_0^\pi \frac{\sin \theta}{\beta_r + \alpha_r k^2 (4\Omega_0^2 \cos^2 \theta + N^2 \sin^2 \theta)} d\theta. \quad (69) $$

The integration yields (see also, Herbert et al.³³),

$$ \langle E(k) \rangle_{np} = \begin{cases} \frac{2\pi k}{\sqrt{\alpha_r(\beta_r+\alpha_r k^2 N^2)(N^2-4\Omega_0^2)}} \ln \left( \frac{1+k\sqrt{\alpha_r(N^2-4\Omega_0^2)/(\beta_r+\alpha_r k^2 N^2)}}{1-k\sqrt{\alpha_r(N^2-4\Omega_0^2)/(\beta_r+\alpha_r k^2 N^2)}} \right) & \text{if } 2\Omega_0 < N \\ \frac{4\pi k^2}{(\beta_r+\alpha_r k^2 N^2)} & \text{if } 2\Omega_0 = N \\ \frac{4\pi k}{\sqrt{\alpha_r(\beta_r+\alpha_r k^2 N^2)(4\Omega_0^2-N^2)}} \tan^{-1} \sqrt{\frac{\alpha_r k^2 (4\Omega_0^2-N^2)}{(\beta_r+\alpha_r k^2 N^2)}} & \text{if } 2\Omega_0 > N \end{cases} \quad (70) $$

and hence,

$$ \langle \Gamma_q^{(\kappa)}(k) \rangle_{np} = \frac{4\pi k^2}{\alpha_r} - \frac{\beta_r}{\alpha_r} \langle E(k) \rangle_{np}. \quad (71) $$

Consider the non-magnetized stratified and rotating case with $N = 2\Omega_0 = 1$, for example. The solution (70) for $\langle E(k) \rangle_{np}$ becomes divergent if $k^2 = -\beta_r/\alpha_r \ge k_{min}^2$. Consequently, for the choice $\beta_r = -0.9999 < 0$ and $\alpha_r = 100 > 0$ (i.e., regime I) there is an infrared (IR) divergence for $k_{min}$ close to 0.1 (see figure 2a). This indicates that a condensation of energy in the gravest modes should be expected, and suggests the presence of an inverse cascade for the vortex modes (Herbert et al.³³) similar to the inverse cascade of 2D turbulence (Kraichnan & Montgomery⁴¹). On the other hand, for the choice $\beta_r = 1000 > 0$ and $\alpha_r = -0.1 < 0$ (regime III) there is an ultraviolet (UV) divergence for $k_{max}$ close to 100, as illustrated by figure 2b. Note that the $2\Omega_0 < N$ and $2\Omega_0 > N$ cases exhibit similar behaviour.

For the stratified MHD case, the analysis of the equation (67) indicates that $\langle E(k) \rangle_{np}$ is not defined for $\alpha < 0$, independently of the wavenumber $k$. Indeed, suppose that $\alpha < 0$, then the quantity $\beta + \alpha N^2$, which is independent of $k$, must be of negative sign if $k \ne N/B_0$ in order to satisfy the first and third equality in (67), while at $k = N/B_0$, it must be of positive sign in order to satisfy the second equality in (67). Thus, for the magnetized system and unlike the non-magnetized system, the regime (III) (i.e., $\alpha < 0$ and $\beta > 0$) does not exist.

On the other hand, for the admissible values of the Lagrange multipliers $(\alpha, \beta)$ such that $\alpha > 0$ and

$$ -\min(B_0^2 k_{min}^2, N^2) \alpha < \beta < 0 \quad (72) $$

(i.e., regime I)), there is no an IR divergence. Indeed, for a given value of $(N, B_0)$ there exists a wavenumber $k$ such that $k = k_c = N/B_0$ for which $\langle E(k_c) \rangle_{np} \to +\infty$ if $\beta \to (-\alpha/N^2)^+ < 0$, but, at the same time, equation (67) implies that $\langle E(k) \rangle_{np} \to +\infty$ for all $k > k_c$, and $\langle E(k) \rangle_{np}$ is not defined for $k < k_c$. For $\alpha > 0$ and $\beta > 0$ (regime II) the condition (57) is automatically satisfied.

As an illustration, figure 3 shows $\langle E(k) \rangle_{np}$ versus $k$ for $B_0 = 10$, $N = 5$, $\alpha = 100$ and $\beta = -0.9$ (regime I, figure 3(a)), and for $B_0 = 1$, $N = 1$, $\alpha = 1$ and $\beta = 2$ (regime II, figure 3(b)). For regime I (see figure 2(a)), $\langle E(k) \rangle_{np}$ behaves like $k$ at large $k(> N/B_0)$, whereas at small $kB_0/N$ it scales like $k^{1.7}$. For regime (II), the energy scales like $k^2$ (equipartition in three dimensions) for $k \le N/B_0$, and increases linearly with $k$ at large $k(> N/B_0)$. We also show in figure 2 the variation of $\langle \Gamma_q(k) \rangle_{np}$, which behaves like $k^2$ at small and large scales.

Accordingly, we may conclude that, in contrast to the case of stratified and rotating flows, in stratified MHD turbulence with $B_0 \ne 0$, the equilibrium energy spectra rather points to a forward cascade of energy even in the restricted to the submanifold of non-propagating modes $\Lambda^{(np)}$.
<figure><img src="image_3.png" /><figcaption>FIG. 2. Radial spectra of energy $\langle E(k) \rangle_{np}$ (solid black lines) and the quadratic part of potential enstrophy $\langle \Gamma_q^{(\kappa)}(k) \rangle_{np}$ (solid blue lines) for non-magnetized stratified and rotating turbulence ($N = 2\Omega_0 = 1$) at statistical equilibrium restricted to the submanifold of non-propagating modes $\Lambda^{(np)}$, in different ($\alpha_r, \beta_r$) regimes: (a) $\alpha_r = 100 > 0$ and $\beta_r = -0.9999 < 0$) (regime I) for which there is an infrared divergence); (b) for $\alpha_r = -0.1 < 0$ and $\beta_r = 1000 > 0$ (regime III) for which there is an ultraviolet divergence).</figcaption></figure>

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

V. CONCLUDING REMARKS

Quasi-linear theory has been used to study the effects of stable ambient density stratification and imposed mean magnetic field on homogeneous turbulence viewed in frame rotating uniformly about a vertical axis for an electrically conducting Boussinesq fluid. Here, the mean magnetic field ($B_0 = B_0\hat{z}$) aligns with the buoyancy gradient of strength $N$, as well with the rotation axis ($\Omega_0 = \Omega_0\hat{z}$).

The main results of the present study can be quoted as follows.

i) For the stratified and rotating MHD case, the $L_2$ norm of MIPS, which is an inviscid invariant, can be approximated by its quadratic part $\Gamma_q$ if the Froude and Alfvén-Mach numbers are sufficiently small, $Fr \to 0$, $M \to 0$. The $L_2$ norm of MIPS is similar to the potential enstrophy (i.e. the $L_2$ norm of PV) which is an inviscid invariant for the non-magnetized rotating and stratified case, and it can be approximated by its quadratic part $\Gamma_q^{(\kappa)}$ if the Froude and Rossby numbers are sufficiently small, $Fr \to 0$, $R_0 \to 0$. For each one of these systems

<figure><img src="image_5.png" /><figcaption>FIG. 3. Radial spectra of energy $\langle E(k) \rangle_{np}$ (solid black lines) and $L_2$ norm of MIPS $\langle \Gamma_q(k) \rangle_{np}$ (solid blue lines) for stratified and rotating MHD turbulence with a mean magnetic field at statistical equilibrium restricted to the submanifold of non-propagating modes $\Lambda^{(np)}$, in different ($\alpha, \beta$) regimes: (a) $\alpha = 100 > 0$ and $\beta = -0.9$ (regime I) with $N = 5$ and $B_0 = 10$; (b) $\alpha = 1 > 0$ and $\beta = 2 > 0$ (regime II) with $N = 1$ and $B_0 = 1$.</figcaption></figure>
(magnetized or not) the total energy is an inviscid invariant.

ii) A normal mode decomposition has been performed in the case of stratified and rotating MHD turbulence. It indicates the existence of three types of motion: slow and fast magneto-inertia-gravity waves with frequency $\omega_s$ and $\omega_f$ (see Eq. (34)) and a non-propagating mode which does not explicitly depends on rotation. From the expression of these normal modes we recover that of those given in previous study of the non-magnetized rotating and stratified case for which the non-propagating mode explicitly depends on rotation.

iii) We have characterized the two systems (magnetised or not) setting up our investigation from the point of view of equilibrium statistical mechanics where the integral defining the partition function ($\mathcal{Z}^{(np)}(\mathbf{k})$) is restricted to the slow manifold. From the expression of $\mathcal{Z}^{(np)}(\mathbf{k})$, we have derived the radial spectra $\langle E(k) \rangle_{np}$, $\langle \Gamma_q(k) \rangle_{np}$ and $\langle \Gamma_q^{(\kappa)}(k) \rangle_{np}$ (see Eqs. (67) and (70)).

iv) In the case of the non-magnetised system there is an infrared divergence which indicates that a condensation of energy in the gravest modes should be expected. This suggests the presence of an inverse cascade for the vortex modes (see, Herbert et al. 33). In contrast, in the case of the stratified and rotating MHD system with $B_0 \neq 0$, we showed that the partition function restricted to the slow manifold does not indicate an inverse cascade of energy. Note that recent DNS study by Salhi et al. 17 for the non-rotating stratified MHD turbulence with a weak mean magnetic field also indicates a forward cascade even for the magnetic field for which the vertical magnetic energy is drastically reduced by stratification.

## AUTHOR DECLARATIONS

### Conflict of Interest

The authors have no conflicts to disclose.

## APPENDIX

### A.1. Eigenvectors of the matrix $\mathcal{L}$

The eigenvectors $\mathbf{Z}_0$, $\mathbf{Z}_s^{\pm}$ and $\mathbf{Z}_f^{\pm}$ associated to the eigenvalues $\lambda_1 = 0$, $\lambda_{2,3} = \pm \omega_s$, $\lambda_{4,5} = \pm \omega_f$ (see equation

(37)), respectively, are of the form

$$ \mathbf{Z}_0(\mathbf{k}) = \frac{1}{\Psi_0} \begin{pmatrix} 0 \\ 0 \\ 0 \\ \omega_g \\ \text{is} \omega_a \end{pmatrix}, \qquad (73a) $$

$$ \mathbf{Z}_s^{\pm}(\mathbf{k}) = \frac{1}{\Psi_s} \begin{pmatrix} \pm \omega_s^2 \\ -i\omega_s^2 \\ -s\omega_a\omega_s \\ \pm \text{is} \omega_a \sigma_s \\ \pm \omega_g \sigma_s \end{pmatrix} \qquad (73b) $$

$$ \mathbf{Z}_f^{\pm}(\mathbf{k}) = \frac{1}{\Psi_f} \begin{pmatrix} i\omega_f \sigma_f \\ \pm \omega_f^2 \\ \mp \text{is} \omega_a \sigma_f \\ -s\omega_a \omega_f \\ i\omega_g \sigma_f \end{pmatrix} \qquad (73c) $$

where

$$ \Psi_0 = \sqrt{\omega_a^2 + \omega_g^2}, \qquad (74a) $$

$$ \Psi_s = \sqrt{2\omega_s^2 (\sigma_s^2 + \omega_a^2)}, \qquad (74b) $$

$$ \Psi_f = \sqrt{2\omega_f^2 (\sigma_f^2 + \omega_a^2 + \omega_g^2)}. \qquad (74c) $$

These vectors are normalised and are orthogonal two by two, and are the columns of the matrix $\mathbf{P}$ (see equation (36)). They form an orthonormal basis of $\mathbb{C}^5$.

### A.2. The limit case: $\omega_a(\mathbf{k}) \to 0, \hat{b}(\mathbf{k}) \to 0, \forall \mathbf{k} \in \mathbb{R}^3$.

In the limit $\omega_a(\mathbf{k}) \to 0, \hat{b}(\mathbf{k}) \to 0, \forall \mathbf{k} \in \mathbb{R}^3$, the pulsations $\omega_s$ (see Eq. (37b)) and $\sigma_s$ (see Eq. (40)) as well as $\Psi_s(\mathbf{k})$ tend towards zero

$$ \lim_{\omega_a \to 0} \omega_s(\mathbf{k}) = 0, \quad \lim_{\omega_a \to 0} \sigma_s(\mathbf{k}) = 0, \quad \lim_{\omega_a \to 0} \Psi_s(\mathbf{k}) = 0. $$

In this case, there are indeterminate forms for the unit vectors $\mathbf{Z}_s^{\pm}$. To overcome this indeterminate form, we consider the following equivalent form:

$$ \mathbf{Z}_s^{\pm} = \frac{1}{\|\mathbf{Z}_a^{\pm}\|} \mathbf{Z}_a^{\pm}, \qquad (75a) $$

$$ \mathbf{Z}_a^{\pm} = \frac{1}{(\omega_a^2 + \omega_g^2 - \omega_s^2)} \begin{pmatrix} (\omega_a^2 + \omega_g^2 - \omega_s^2) \\ \pm \text{is} \omega_r \omega_s \\ \mp \frac{s\omega_a}{\omega_s} (\omega_a^2 + \omega_g^2 - \omega_s^2) \\ -i\omega_a \omega_r \\ -s\omega_r \omega_g \end{pmatrix}, \qquad (75b) $$

where

$$ \|\mathbf{Z}_a^{\pm}\|^2 = 1 + \frac{\omega_a^2}{\omega_s^2} + \frac{(\omega_a^2 + \omega_g^2 + \omega_s^2)}{(\omega_a^2 + \omega_g^2 - \omega_s^2)^2} \omega_r^2. \qquad (76) $$
On the other hand, we take into account the orthogonality between $Z_a^+$ and $Z_a^-$, so that

$$Z_a^+ \cdot Z_a^- = 1 - \frac{\omega_a^2}{\omega_s^2} + \frac{\omega_r^2}{(\omega_a^2 + \omega_g^2 - \omega_s^2)} = 0. \quad (77)$$

It follows that

$$\|Z_a^\pm\|^2 = 2 \left( 1 + \frac{(\omega_a^2 + \omega_g^2)\omega_r^2}{(\omega_a^2 + \omega_g^2 - \omega_s^2)^2} \right) \quad (78a)$$

$$= 2 \left( \frac{\omega_a^2}{\omega_s^2} + \frac{\omega_r^2 \omega_s^2}{(\omega_a^2 + \omega_g^2 - \omega_s^2)^2} \right) \quad (78b)$$

and

$$\lim_{\omega_a \to 0} \|Z_a^\pm\|^2 = 2 \left( 1 + \frac{\omega_r^2}{\omega_g^2} \right), \quad (79a)$$

$$\lim_{\omega_a \to 0} \frac{\omega_a^2}{\omega_s^2} = 1 + \frac{\omega_r^2}{\omega_g^2}. \quad (79b)$$

Accordingly, we deduce the limits given by Eq. (45b).

$$\lim_{\omega_a \to 0} \frac{\sigma_s^2 \omega_g^2}{\omega_s^2 (\sigma_s^2 + \omega_a^2)} = \lim_{\omega_a \to 0} \frac{(1 - \omega_a^2/\omega_s^2)^2 \omega_g^2}{(\omega_r^2 \omega_a^2/\omega_s^2 + \omega_s^2 (1 - \omega_a^2/\omega_s^2))^2} \quad (80a)$$

$$= \frac{\omega_r^2}{\omega_r^2 + \omega_g^2}, \quad (80b)$$

$$\lim_{\omega_a \to 0} \frac{\omega_s^2}{(\sigma_s^2 + \omega_a^2)} = \lim_{\omega_a \to 0} \frac{\omega_r^2}{(\omega_r^2 \omega_a^2/\omega_s^2 + \omega_s^2 (1 - \omega_a^2/\omega_s^2))^2} \quad (80c)$$

$$= \frac{\omega_g^2}{\omega_r^2 + \omega_g^2}. \quad (80d)$$

REFERENCES

¹R. Hide, Free hydromagnetic oscillations of the Earth's core and the theory of the geomagnetic secular variation, Philos. Trans. R. Soc. London A: Math. Phys. Eng. Sci. 259, 615 (1966).

²P. Kumar, S. Talon, and J.-P. Zahn, Angular momentum transport by gravity waves and its effect on the rotation of the solar interior, Astrophys. J. 520, 859 (1999).

³P. H. Roberts and A. M. Soward, Magnetohydrodynamics of the Earth's core, Annu. Rev. Fluid Mech. 4, 117 (1972).

⁴H. K. Moffatt, Field Generation in Electrically Conducting Fluids (Cambridge University Press, Cambridge/London/New York/Melbourne, 1978).

⁵C. C. Finlay, Waves in the presence of magnetic fields, rotation, and convection, in Dynamos, Vol. LXXXVIII of Lecture Notes of the Les Houches Summer School, edited by P. Cardin and L. Cugliandolo (Elsevier, New York, 2008), Chap. 8.

⁶A. Salhi, F. S. Baklouti, F. Godeferd, T. Lehner & C. Cambon (2017). Energy partition, scale by scale, in magnetic Archimedes Coriolis weak wave turbulence. Phys. Rev. E, 95(2), 023112.

⁷S. D. Grant, D. B. Jess, T. V. Zaqarashvili, C. Beck, H. Socas-Navarro, M. J. Aschwanden, ... & R. L. Hewitt (2018). Alfvén wave dissipation in the solar chromosphere. Nature Physics, 14(5), 480-483.

⁸T. Straus, B. Fleck, S. M. Jefferies, G. Cauzzi, S. W. McIntosh, K. Reardon,... & M. Steffen, M. (2008). The energy flux of internal gravity waves in the lower solar atmosphere. The Astrophysical Journal, 681(2), L125.

⁹F. Wilczyński, D. W. Hughes & E. Kersalé. Magnetic buoyancy instability and the anelastic approximation: regime of validity and relationship with compressible and Boussinesq descriptions. J. Fluid Mech., 942 (2022) A46.

¹⁰T. V. Zaqarashvili, S. Lomineishvili, P. Leitner, A. Hanslmeier, P. Gömör & M. Roth (2021). Kink instability of triangular jets in the solar atmosphere. Astronomy & Astrophysics, 649, A179.

¹¹F. Kneer & N. B. González (2011). On acoustic and gravity waves in the solar photosphere and their energy transport. Astronomy & Astrophysics, 532, A111.

¹²K. Nagashima, B. Löptien, L. Gizon, A. C. Birch, A. C., R. Cameron, S. Couvidat,... & Stein, R. (2014). Interpreting the helioseismic and magnetic imager (HMI) multi-height velocity measurements. Solar Physics, 289, 3457-3481.

¹³G. Vigeesh & M. Roth (2020). Synthetic observations of internal gravity waves in the solar atmosphere. Astronomy & Astrophysics, 633, A140.

¹⁴A. Hague & R. Erdélyi (2016). Buoyancy-driven magnetohydrodynamic waves. The Astrophysical Journal, 828(2), 88.

¹⁵J. Pedlosky (2013). Geophysical fluid dynamics. Springer Science & Business Media (2013).

¹⁶J. J. Riley, R. W. Metcalfe & M. A. Weissman (1981). Direct numerical simulations of homogeneous turbulence in density stratified fluids. In Nonlinear Properties of Internal Waves: La Jolla Institute (Vol. 76, No. 1, pp. 79-112). AIP Publishing.

¹⁷A. Salhi, A. Khlifi, R. Marino, F. Feraco, R. Foldes & C. Cambon (2024). Waves and non-propagating modes in stratified MHD turbulence subject to a weak magnetic field. Accepted (J. Fluid Mech).

¹⁸A. Salhi, T. Lehner, F. Godeferd & C. Cambon. Magnetized stratified rotating shear waves. Phys. Rev. E, 85(2), (2012) 026301.

¹⁹A. Salhi & C. Cambon. Magneto-gravity-elliptic instability. J. Fluid Mech., 963, (2023): A9.

²⁰P. F. Embid & A. J. Majda (1998). Low Froude number limiting dynamics for stably stratified flow with small or finite Rossby numbers. Geophysical & Astrophysical Fluid Dynamics, 87(1-2), 1-50.

²¹S. Kurien, B. Wingate & M. A. Taylor (2008). Anisotropic constraints on energy distribution in rotating and stratified turbulence. Europhysics letters, 84(2), 24003.

²²J. G. Charney (1971). Geostrophic turbulence. Journal of the Atmospheric Sciences, 28(6), 1087-1095.

²³A. Vallgren & E. Lindborg (2010). Charney isotropy and equipartition in quasi-geostrophic turbulence. J. Fluid Mech., 656, 448-457.

²⁴G. D. Nastrom & K. S. Gage (1983). A first look at wavenumber spectra from GASP data. Tellus A, 35(5), 383-388.

²⁵D. K. Lilly (1989). Two-dimensional turbulence generated by energy sources at two scales. Journal of Atmospheric Sciences, 46(13), 2026-2030.

²⁶M. E. Maltrud & G. K. Vallis (1991). Energy spectra and coherent structures in forced two-dimensional and beta-plane turbulence. J. Fluid Mech., 228, 321-342.

²⁷Smith, L. M., & Waleffe, F. (2002). Generation of slow large scales in forced rotating stratified turbulence. J. Fluid Mech., 451, 145-168.

²⁸R. Marino, D. Rosenberg, C. Herbert & A. Pouquet (2015). Interplay of waves and eddies in rotating stratified turbulence and the link with kinetic-potential energy partition. Europhysics Letters, 112(4), 49001.

²⁹C. Herbert, R. Marino, D. Rosenberg & A. Pouquet (2016). Waves and vortices in the inverse cascade regime of stratified turbulence with or without rotation. J. Fluid Mech., 806, 165-204.
<footnote n="30"><sup>30</sup>A. van Kan & A. Alexakis (2022). Energy cascades in rapidly rotating and stratified turbulence within elongated domains. J. Fluid Mech., **933**, A11.</footnote>

<footnote n="31"><sup>31</sup>J. F. Scott & C. Cambon. Evolution of weak, homogeneous turbulence with rotation and stratification. J. of Fluid Mech. 979 (2024): A17.</footnote>

<footnote n="32"><sup>32</sup>P. Bartello (1995). Geostrophic adjustment and inverse cascades in rotating stratified turbulence. Journal of Atmospheric Sciences, 52(24), 4410-4428.</footnote>

<footnote n="33"><sup>33</sup>C. Herbert, A. Pouquet & R. Marino. Restricted equilibrium and the energy cascade in rotating and stratified flows. J. Fluid Mech., **758**, (2014) 374-406.</footnote>

<footnote n="34"><sup>34</sup>A. Salhi & C. Cambon (1997). An analysis of rotating shear flow using linear theory and DNS and LES results. J. Fluid Mech., **347**, 171-195.</footnote>

<footnote n="35"><sup>35</sup>E. A. Spiegel & G. Veronis (1960). On the Boussinesq approximation for a compressible fluid. Astrophys. J., **131**, 442-447.</footnote>

<footnote n="36"><sup>36</sup>J. V. Shebalin (2006). Ideal homogeneous magnetohydrodynamic turbulence in the presence of rotation and a mean magnetic field. J. Plasma Physics, **72**(4), 507-524.</footnote>

<footnote n="37"><sup>37</sup>S. Galtier. Weak turbulence theory for rotating magnetohydrodynamics and planetary flows. J. Fluid Mech., **757** (2014) 114-154.</footnote>

<footnote n="38"><sup>38</sup>M. L. Waite & P. Bartello (2006). The transition from geostrophic to stratified turbulence. J. Fluid Mech., **568**, 89-108.</footnote>

<footnote n="39"><sup>39</sup>A. Altuntaş, L. Davidson & S. H. Peng (2019). A new approximation to modulation-effect analysis based on empirical mode decomposition. Phys. Fluids **31**, 025117.</footnote>

<footnote n="40"><sup>40</sup>G. Holloway (1986). Eddies, waves, circulation, and mixing: Statistical geofluid mechanics. Annual Rev. Fluid Mech., **18**(1), 91-147.</footnote>

<footnote n="41"><sup>41</sup>R. H. Kraichnan & D. Montgomery. Two-dimensional turbulence. Reports on Progress in Physics, **43**(5), (1980) 547.</footnote>

<footnote n="42"><sup>42</sup>C. Herbert. Restricted partition functions and inverse energy cascades in parity symmetry breaking flows. Phys. Rev. E, **89**(1), (2014) 013010.</footnote>