2.1. The Boussinesq MHD equations

The Boussinesq approximation ignores density differences except where they appear in terms multiplied by the acceleration due to the buoyancy force (see, e.g. Spiegel   Veronis Reference Spiegel and Veronis1960; Salhi et al. Reference Salhi, Lehner, Godeferd and Cambon2013). Accordingly, the continuity equation reduces to $\boldsymbol{\nabla }{\boldsymbol{\cdot }} \tilde {\boldsymbol{u}}=0$ as for an incompressible fluid.

In a frame rotating uniformly about a fixed axis at the rate ${{\boldsymbol{\varOmega }}}_{\!p}$ , the Boussinesq MHD equations for the instantaneous velocity field $\tilde {\boldsymbol{u}}({\boldsymbol{x}},t),$ magnetic field $\tilde {\boldsymbol{b}}({\boldsymbol{x}},t)$ and buoyancy scalar $\tilde \vartheta ({\boldsymbol{x}},t)$ can be written as (see, e.g. Davidson Reference Davidson2013)

(2.1a) \begin{align} \left (\partial _t +\tilde {{\boldsymbol{u}}}{\boldsymbol{\cdot }}{\boldsymbol{\nabla }}\right )\!{\tilde {{\boldsymbol{u}}}}&= -\boldsymbol{\nabla }{\tilde p}-2{{\boldsymbol{\varOmega }}}_{\!p}\times \tilde {{\boldsymbol{u}}} +\big (\tilde {{\boldsymbol{b}}}\,{\boldsymbol{\cdot }} {\boldsymbol{\nabla }}\big )\tilde {{\boldsymbol{b}}} +\tilde {\vartheta }{\boldsymbol{n}}+\nu {\nabla} ^2 \tilde {{\boldsymbol{u}}}, \end{align}

(2.1b) \begin{align} \left (\partial _t +\tilde {{\boldsymbol{u}}}{\boldsymbol{\cdot }}{\boldsymbol{\nabla }}\right ){\tilde {{\boldsymbol{b}}}}&= \big (\tilde {{\boldsymbol{b}}}\,{\boldsymbol{\cdot }}{\boldsymbol{\nabla }}\big )\tilde {{\boldsymbol{u}}}+\eta {\nabla} ^2\tilde {{\boldsymbol{b}}}, \end{align}

(2.1c) \begin{align} \left (\partial _t +\tilde {{\boldsymbol{u}}}{\boldsymbol{\cdot }}{\boldsymbol{\nabla }}\right )\!{\tilde {\vartheta }}&=\kappa {\nabla} ^2\tilde {\vartheta }, \end{align}

(2.1d) \begin{align} \boldsymbol{\nabla }{\boldsymbol{\cdot }}\, {\tilde {{\boldsymbol{u}}}}&=0, \end{align}

(2.1e) \begin{align} \boldsymbol{\nabla }{\boldsymbol{\cdot }}\, {\tilde {{\boldsymbol{b}}}}&=0. \end{align}

Here, the magnetic field 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, $\tilde p$ being the total pressure (including the centrifugal potential and the magnetic part) divided by the constant density $\rho _0,$ and ${\boldsymbol{\varOmega }}_{\!p}$ is the rotation vector of the reference frame which is constant and $\nu ,$ $\eta$ and $\kappa$ are positive constants that denote the kinematic viscosity, magnetic diffusivity and thermal diffusivity, respectively, and $\boldsymbol{n}$ is a unit vector that aligns with the basic buoyancy scalar gradient. As indicated by an anonymous referee, the natural systems possess a non-negligible adiabatic gradient, sometimes very strong (unlike in experimental situations), and hence, $\tilde \vartheta$ should be really defined as an excess over the adiabatic profile (influencing also the basic pressure field), so that the equations indeed look like (2.1).

The equation for the absolute vorticity vector, $\tilde {{\boldsymbol{w}}}= \boldsymbol{\nabla }\times \tilde {{\boldsymbol{u}}}+2{\boldsymbol{\varOmega }}_{\!p},$ is obtained by taking the curl of the momentum equation

(2.2) \begin{equation} \partial _t {\tilde {{\boldsymbol{w}}}}= \boldsymbol{\nabla }\times \big (\tilde {{\boldsymbol{u}}}\times \tilde {{\boldsymbol{w}}}+\tilde {{\boldsymbol{j}}}\times \tilde {{\boldsymbol{b}}}\big )+\boldsymbol{\nabla }\tilde \vartheta \times {\boldsymbol{n}}+\nu {\nabla} ^2\tilde {{\boldsymbol{w}}}, \end{equation}

where $\tilde {{\boldsymbol{j}}}=\boldsymbol{\nabla }\times \tilde {{\boldsymbol{b}}}$ is the normalized current density. In addition, we report here the equation for the gradient of the buoyancy scalar,

(2.3) \begin{equation} \left (\partial _t +\tilde {{\boldsymbol{u}}}{\boldsymbol{\cdot }}{\boldsymbol{\nabla }}\right )\big (\boldsymbol{\nabla }{\tilde {\vartheta }}\big )=-\big (\boldsymbol{\nabla }\tilde {{\boldsymbol{u}}}\big )^T{\boldsymbol{\cdot }} \big (\boldsymbol{\nabla }\tilde {\vartheta }\big )+\kappa {\nabla} ^2\big (\boldsymbol{\nabla }\tilde \vartheta \big ), \end{equation}

where $T$ denotes the transpose. Accordingly, from (2.2) and (2.3) we deduce the equation for the potential vorticity, $\tilde \varPi _\kappa = \tilde {{\boldsymbol{w}}}{\boldsymbol{\cdot }}{\boldsymbol{\nabla }} \tilde {\vartheta },$ which is very useful as an invariant in stratified geophysical flows (see, e.g. Pedlosky Reference Pedlosky2013),

(2.4) \begin{equation} \left (\partial _t +\tilde {{\boldsymbol{u}}}{\boldsymbol{\cdot }}{\boldsymbol{\nabla }}\right )\!{\tilde \varPi _\kappa }= \big (\boldsymbol{\nabla }\tilde \vartheta \big ){\boldsymbol{\cdot }} \big ( \boldsymbol{\nabla }\times \big (\tilde {{\boldsymbol{j}}}\times \tilde {{\boldsymbol{b}}}\big )\big )+ \nu \big (\boldsymbol{\nabla }\tilde \vartheta \big ){\boldsymbol{\cdot }} {\nabla} ^2\tilde {\boldsymbol{w}}+ \kappa \tilde {\boldsymbol{w}}{\boldsymbol{\cdot }}{\nabla} ^2\big (\boldsymbol{\nabla }\tilde \vartheta \big ). \end{equation}

It follows that, in the case of a magnetized Boussinesq ideal fluid, $\tilde \varPi _\kappa$ is not a Lagrangian invariant: while it removes the baroclinic torque, it does not remove the Lorentz torque. In contrast, the so-called MIPS, $\tilde \varPi _m=\tilde {\boldsymbol{b}}{\boldsymbol{\cdot }} \boldsymbol{\nabla }\tilde \vartheta$ (see, Salhi et al. Reference Salhi, Lehner, Godeferd and Cambon2012), which is a solution of the following transport equation:

(2.5) \begin{equation} \left (\partial _t +\tilde {{\boldsymbol{u}}}{\boldsymbol{\cdot }}{\boldsymbol{\nabla }}\right )\!{\tilde \varPi _m}= \eta \big (\boldsymbol{\nabla }\tilde \vartheta \big ){\boldsymbol{\cdot }} {\nabla} ^2\tilde {\boldsymbol{b}}+ \kappa \tilde {\boldsymbol{b}}{\boldsymbol{\cdot }}{\nabla} ^2\big (\boldsymbol{\nabla }\tilde \vartheta \big ), \end{equation}

deduced from (2.1b ) and (2.3), is a Lagrangian invariant for a magnetized Boussinesq ideal fluid. As will be shown later, the use of the magnetic invariant (i.e. MIPS) makes it possible to reduce the linear system of ordinary differential equations, which represents a Floquet problem of the inviscid stability problem, from a five-dimensional system to a four-dimensional one only.

2.2. The unperturbed system

The solutions of system (2.1) and (2.2) are conveniently decomposed into a ‘basic flow’ ( ${\boldsymbol{U}},P,{\boldsymbol{W}}, {\boldsymbol{B}},\varTheta )$ and a `disturbance’ $({\boldsymbol{u}},p,{\boldsymbol{w}}, {\boldsymbol{b}},\vartheta )$ as

(2.6) \begin{equation} \tilde {\boldsymbol{u}}={\boldsymbol{U}}+{\boldsymbol{u}},\quad \tilde p=P+ p,\quad \tilde {\boldsymbol{w}}= {{\boldsymbol{W}}}+{\boldsymbol{w}},\quad \tilde {\boldsymbol{b}}={\boldsymbol{B}}_0+{\boldsymbol{b}},\quad \tilde \vartheta =\varTheta +\vartheta , \end{equation}

where the basic flow is a solution of system (2.1). As indicated in § 1, we consider the unbounded basic flows KBF and MBF given by (1.1) and (1.3), respectively, that can be thought to represent a small patch of a stratified precessing magnetized flow subject to a background magnetic field.

To discuss the admissibility conditions, i.e. the basic flow must be a solution of the system (2.1) (see, Craik Reference Craik1989), it is convenient to consider the equation for the basic absolute vorticity (2.2), the induction (2.1b ) and the equation for the basic buoyancy scalar gradient (2.3) supplemented by (2.1d , e ),

(2.7a) \begin{align} \left (\partial _t+{\boldsymbol{U}}{\boldsymbol{\cdot }} \boldsymbol{\nabla }\right )\!{\boldsymbol{W}}&={\boldsymbol{W}}{\boldsymbol{\cdot }}\boldsymbol{\nabla }{\boldsymbol{U}}+ \boldsymbol{\nabla }\times \left (\left ({\boldsymbol{B}}_0{\boldsymbol{\cdot }}\boldsymbol{\nabla }\right )\!{\boldsymbol{B}}_0\right )+\nu {\nabla} ^2{\boldsymbol{W}}, \end{align}

(2.7b) \begin{align} \left (\partial _t +{{\boldsymbol{U}}}{\boldsymbol{\cdot }}{\boldsymbol{\nabla }}\right )\!{{{\boldsymbol{B}}_0}}&={\boldsymbol{B}}_0{\boldsymbol{\cdot }}\boldsymbol{\nabla }{\boldsymbol{U}}+\eta {\nabla} ^2\!{{\boldsymbol{B}}_0}, \end{align}

(2.7c) \begin{align} \left (\partial _t +{{\boldsymbol{U}}}{\boldsymbol{\cdot }}{\boldsymbol{\nabla }}\right )\!\left (\boldsymbol{\nabla }{{\varTheta }}\right )&=- \!\left (\boldsymbol{\nabla }{\varTheta }\right )\!{\boldsymbol{\cdot }} \boldsymbol{\nabla }\!\left ({\boldsymbol{U}}\right )^T+\kappa {\nabla} ^2\!\left (\boldsymbol{\nabla }\varTheta \right )\! , \end{align}

with $\boldsymbol{\nabla }{\boldsymbol{\cdot }} {\boldsymbol{U}}=0$ and $\boldsymbol{\nabla }{\boldsymbol{\cdot }} {\boldsymbol{B}}_0=0.$ Recall that the unit vector $\boldsymbol{n}$ aligns with the basic buoyancy scalar gradient, so that, $\boldsymbol{\nabla }\varTheta \times {\boldsymbol{n}}={{\boldsymbol{0}}}.$

2.2.1. Basic velocity

Recall that the KBF, i.e. the velocity field of the precessing flow described by (1.1), approximates the laminar Poincaré flow far from the boundaries of a spheroid if $\varepsilon \ll 1,$ and the tidal deformation is neglected (see, Poincaré Reference Poincaré1910; Kerswell Reference Kerswell1993; Salhi   Cambon Reference Salhi and Cambon2009). In this case, the basic absolute vorticity vector, ${\boldsymbol{W}}=\boldsymbol{\nabla }\times {\boldsymbol{U}}+2\varOmega _{\!p}=2\varOmega _0 (2\varepsilon \hat {\boldsymbol{x}}+\hat {\boldsymbol{z}} ),$ does not align with the solid body rotation. In counterpart, the MBF, i.e. the velocity field of the precessing flow described by (1.3), characterizes an infinite column in which a tilted (sheared) streamline solution can exist under precession (see, Mahalov Reference Mahalov1993; Salhi   Cambon Reference Salhi and Cambon2009). In this case, ${\boldsymbol{W}}=2\varOmega _0\hat {\boldsymbol{z}}$ aligns with the solid body rotation.

Accordingly, in the pure hydrodynamic case, so that ${\boldsymbol{B}}_0$ and $\boldsymbol{\nabla }\varTheta$ are zero, both KBF and MBF are exact solutions of (2.7a ). Note that, the main difference between these two base flows is the cross-gradient direction of the plane shear: for the KBF it is along $\hat {\boldsymbol{y}}$ and varies linearly with respect to the vertical space coordinate $z,$ while for the MBF it is along the vertical axis $\hat {\boldsymbol{z}}$ and varies linearly with respect to the horizontal space coordinate  $y.$ Note that the trajectory of a fluid particle can be determined by resolving the equation $d_t{\boldsymbol{x}}=\boldsymbol{\nabla }{\boldsymbol{U}}{\boldsymbol{\cdot }}{\boldsymbol{x}}.$ Accordingly, for the KBF, for example, we find

(2.8) \begin{equation} x(\tau )=x_{0}\cos \tau -y_{0}\sin \tau ,\ y(\tau )=x_{0}\sin \tau +y_{0}\cos \tau ,\ z(\tau )=z_0+2\varepsilon (x-x_0), \end{equation}

where $\tau =\varOmega _0t$ is a dimensionless time and ${\boldsymbol{x}}_0=(x_0,y_0,z_0)^T$ is the initial (at $t=0)$ vector position. It follows that the fluid particle describes an ellipse in the $(x^*,y^*)\hbox{-}$ plane,

(2.9a) \begin{align}&\qquad\qquad {\boldsymbol{x}}=x\hat {\boldsymbol{x}}+y\hat {\boldsymbol{y}}+z\hat {\boldsymbol{z}}=x^*\hat {\boldsymbol{x}}^*+y^*\hat {\boldsymbol{y}}^*+z^*\hat {\boldsymbol{z}}^*, \end{align}

(2.9b) \begin{align} x^*&=\frac {1}{\sqrt {1+4\varepsilon ^2}}\left (x+2\varepsilon z\right )\! ,\quad y^*=y,\quad z^*=\frac {1}{\sqrt {1+4\varepsilon ^2}}\left (-2\varepsilon x+z\right )\! . \end{align}

Indeed, one easily verifies that $z^*$ (given by (2.9b )) is time-independent and

(2.10) \begin{equation} \frac {1}{\sqrt {1+4\varepsilon ^2}}\left (x^*-2\varepsilon z^*\right )^2+y^{*2}=\text{const}. \end{equation}

On the other hand, we also note that in some previous studies (see, Mason   Kerswell Reference Mason and Kerswell2002; Barker Reference Barker2016b ; Pizzi et al. Reference Pizzi, Mamatsashvili, Barker, Giesecke and Stefani2022) the `mantle frame’ has been used: this consists in applying the following transformation to the solution (1.1):

(2.11) \begin{equation} x_1=x \cos t+y\sin t,\quad x_2=-x\sin t+y\cos t,\quad x_3=z. \end{equation}

2.2.2. Basic magnetic field and buoyancy scalar

Given the solutions (1.1) and (1.3), a constant basic magnetic field ${\boldsymbol{B}}_0,$ which is a trivial solution of (2.7a ), must also satisfy (2.7b ), so that, ${\boldsymbol{B}}_0{\boldsymbol{\cdot }} \boldsymbol{\nabla }{\boldsymbol{U}}=0,$ and hence, ${\boldsymbol{B}}_0$ aligns with the basic absolute vorticity:

(2.12a) \begin{align} {\boldsymbol{B}}_0&=B_0\!\left (2\varepsilon \hat {\boldsymbol{x}}+\hat {\boldsymbol{z}}\right )\quad \text{for the KBF}, \end{align}

(2.12b) \begin{align} {\boldsymbol{B}}_0&=B_0\hat {\boldsymbol{z}}\quad \text{for the MBF}. \end{align}

Here, $B_0$ is a positive constant.

Given the above solutions for ( ${\boldsymbol{U}}, {\boldsymbol{B}}_0),$ in a stratified MHD case under the Boussinesq approximation the basic buoyancy scalar must be a solution of the (2.7c ).

For the KBF (given by (1.1)), a vertical constant basic density gradient, $\boldsymbol{\nabla }\varrho =({\rm d}\varrho /{\rm d}z)\hat {\boldsymbol{z}},$ so that

(2.13) \begin{equation} \varTheta =-\frac {g}{\rho _0}\left (\frac {{\rm d}\varrho }{{\rm d}z}\right ) z=N^2z,\quad N^2=-\frac {g}{\rho _0}\frac {{\rm d}\varrho }{{\rm d}z},\quad {\boldsymbol{n}}=\frac {\boldsymbol{\nabla }\varTheta }{\Vert {\boldsymbol{\nabla }\varTheta }\Vert }=\hat {\boldsymbol{z}} \end{equation}

is an admissible solution. Here $N$ being the Brunt–Väisälä frequency ( $N^2\gt 0$ for a stable stratification).

For the MBF (given by (1.3)), we consider that, a priori, the basic buoyancy scalar varies linearly with the space coordinates (see, e.g. Craik Reference Craik1989), $\varTheta =(\partial _x\varTheta )x+ (\partial _y\varTheta )y+N^2 {\boldsymbol{z}},$ so that, without precession, it reduces to $\varTheta =N^2 {\boldsymbol{z}}.$ By substituting the above form into (2.7c ), we obtain, $ (\boldsymbol{\nabla }{\boldsymbol{U}} )^T{\boldsymbol{\cdot }} (\boldsymbol{\nabla }\varTheta )={{\boldsymbol{0}}},$

(2.14) \begin{equation} \left (\!\!\begin{array}{ccc} 0&1&0\\ -1&0&-2\varepsilon \\ 0&0&0\end{array}\!\!\right )\!{\boldsymbol{\cdot }}\left ( \!\!\begin{array}{c} \partial _x\varTheta \\ \partial _y\varTheta \\ N^2\end{array}\!\!\right )=\left (\!\!\begin{array}{c} 0\\0\\0 \end{array}\!\!\right ) \end{equation}

or equivalently,

(2.15) \begin{equation} \partial _y\varTheta =0\quad \text{and}\quad \partial _x\varTheta =-2\varepsilon N^2. \end{equation}

Accordingly, the basic buoyancy scalar takes the form

(2.16) \begin{equation} \varTheta ({\boldsymbol{x}})={N^2}\! \left (-2\varepsilon x+z\right )\! ,\quad {\boldsymbol{n}}=\frac {\boldsymbol{\nabla }\varTheta }{\Vert {\boldsymbol{\nabla }\varTheta }\Vert }=\frac {1}{\sqrt {1+4\varepsilon ^2}}\left (-2\varepsilon \hat {\boldsymbol{x}}+\hat {\boldsymbol{z}}\right )\! . \end{equation}

The presence of the additional horizontal component of the basic buoyancy scalar, which vanishes in the case without precession $(\varepsilon =0),$ is thereby due to the gyroscopic torque as for the additional mean shear. In a cylindrical coordinate system $(r,\varphi ,z)$ the solution (2.16) can be rewritten as

(2.17) \begin{equation} \varTheta ({\boldsymbol{x}})={N^2}\! \left (-2\varepsilon \left (\cos \varphi \right )\! r+z\right )\! ,\quad {\boldsymbol{n}}=\frac {1}{\sqrt {1+4\varepsilon ^2}}\left [-2\varepsilon \left ((\cos \varphi )\hat {\boldsymbol{r}}-(\sin \varphi )\hat {\boldsymbol{\varphi }}\right )+\hat {\boldsymbol{z}}\right ]\!. \end{equation}

Thus, we recover the solution proposed by Benkacem et al. (Reference Benkacem, Salhi, Khlifi, Nasraoui and Cambon2022) who extended Mahalov’s theoretical study (Mahalov Reference Mahalov1993) by taking into account the effect of stable stratification on precessing columns: they proposed a first-order correction to the profile of the basic buoyancy scalar because the Coriolis force alters the base flow (see their Appendix A).

For the sake of clarity, we repeat here the base flows considered in the present study,

(2.18) \begin{equation} \boldsymbol{\nabla }{\boldsymbol{U}}=\varOmega _0\!\left (\!\!\begin{array}{ccc} 0&-1&0\\ 1&0&-2\varepsilon \\ 0&0&0 \end{array}\!\!\right )\! , \ {\boldsymbol{\varOmega }}_{\!p}= \varOmega _0\!\left (\!\!\begin{array}{c} \varepsilon \\ 0\\ 0 \end{array}\!\!\right )\! , \ {\boldsymbol{B}}_0= {B_0}\left (\!\!\begin{array}{c} 2\varepsilon \\ 0\\ 1 \end{array}\!\!\right )\! ,\ \boldsymbol{\nabla }\varTheta = N^2\!\left (\!\!\begin{array}{c} 0\\ 0\\ 1 \end{array}\!\!\right )\! ,\end{equation}

for the KBF and

(2.19) \begin{equation} \boldsymbol{\nabla }{\boldsymbol{U}}=\varOmega _0\!\left (\!\!\begin{array}{ccc} 0&-1&0\\ 1&0&0\\ 0&-2\varepsilon &0 \end{array}\!\!\right )\! , \ {\boldsymbol{\varOmega }}_{\!p}= \varOmega _0\!\left (\!\!\begin{array}{c} \varepsilon \\ 0\\ 0 \end{array}\!\!\right )\! , \ {\boldsymbol{B}}_0= B_0\!\left (\!\!\begin{array}{c} 0\\ 0\\ 1 \end{array}\!\!\right )\! ,\ \boldsymbol{\nabla }\varTheta = {N^2}\left (\!\!\begin{array}{c} -2\varepsilon \\ 0\\ 1 \end{array}\!\!\right )\! ,\end{equation}

for the MBF.

3.1. Linearized system in physical space

We substitute the solutions (2.6) into the system (2.1) and linearize. The resulting perturbed equations are

(3.1a) \begin{align} \left (\partial _t +{{\boldsymbol{U}}}{\boldsymbol{\cdot }}{\boldsymbol{\nabla }}\right )\!{{\boldsymbol{u}}}&= -\boldsymbol{\nabla }{p}-{\boldsymbol{u}}{\boldsymbol{\cdot }}\boldsymbol{\nabla }{\boldsymbol{U}}-2\varepsilon \varOmega _0\hat {\boldsymbol{x}}\times {{\boldsymbol{u}}} +\left ({\boldsymbol{B}}_0{\boldsymbol{\cdot }} {\boldsymbol{\nabla }}\right )\!{{\boldsymbol{b}}} +{\vartheta }{\boldsymbol{n}}, \end{align}

(3.1b) \begin{align} \left (\partial _t +{{\boldsymbol{U}}}{\boldsymbol{\cdot }}{\boldsymbol{\nabla }}\right )\!{{\boldsymbol{b}}}&={\boldsymbol{b}}{\boldsymbol{\cdot }}\boldsymbol{\nabla }{\boldsymbol{U}}+ \left ({\boldsymbol{B}}_0{\boldsymbol{\cdot }}{\boldsymbol{\nabla }}\right )\!{{\boldsymbol{u}}}, \end{align}

(3.1c) \begin{align} \left (\partial _t +{{\boldsymbol{U}}}{\boldsymbol{\cdot }}{\boldsymbol{\nabla }}\right )\!{\vartheta }&=-{\boldsymbol{u}}{\boldsymbol{\cdot }}\boldsymbol{\nabla }\varTheta, \end{align}

(3.1d) \begin{align} \boldsymbol{\nabla }{\boldsymbol{\cdot }}\, {{{\boldsymbol{u}}}}&=0, \end{align}

(3.1e) \begin{align} \boldsymbol{\nabla }{\boldsymbol{\cdot }}\, {{{\boldsymbol{b}}}}&=0. \end{align}

For an inviscid and non-diffusive Boussinesq fluid, the potential induction magnetic scalar (MIPS) $\tilde \varPi _m=\varPi _m+\varpi _m+\varpi _m^{(nl)}$ is a Lagrangian invariant, as already indicated (see (2.5)), and its linear part takes the form

(3.2) \begin{equation} \varpi _m={\boldsymbol{b}}{\boldsymbol{\cdot }}\boldsymbol{\nabla }\varTheta +{\boldsymbol{B}}_0{\boldsymbol{\cdot }} \boldsymbol{\nabla }\vartheta \, \end{equation}

while $\varpi _m^{(nl)}={\boldsymbol{b}}{\boldsymbol{\cdot }} \boldsymbol{\nabla }\vartheta$ denotes its nonlinear part which is omitted in the present study.

3.2. Floquet system in wave space

The disturbances are expressed in terms of plane waves, for which the direction and the speed of propagation depend on time (see, Bayly Reference Bayly1986; Craik   Criminale Reference Craik and Criminale1986),

(3.3) \begin{equation} \big ({\boldsymbol{u}},p,{\boldsymbol{b}},N^{-1}\vartheta \big )({\boldsymbol{x}},t)=\big (\hat {\boldsymbol{u}},\hat p,\hat {\boldsymbol{b}},\hat \vartheta \big )({\boldsymbol{k}},t)\exp \!\left (\textrm{i}{\boldsymbol{k}}(t){\boldsymbol{\cdot }}{\boldsymbol{x}}\right ) \end{equation}

where $\textrm{i}^2=-1,$ ${\boldsymbol{k}}(t)$ is the wavevector and $\hat \vartheta ({\boldsymbol{k}})$ has a dimension of a velocity. Accordingly, the material derivative of the fluctuating velocity $\boldsymbol{u}$ can be rewritten as

(3.4) \begin{equation} \left (\partial _t +{{\boldsymbol{U}}}{\boldsymbol{\cdot }}{\boldsymbol{\nabla }}\right )\!{{\boldsymbol{u}}}({\boldsymbol{x}},t)=\left [\partial _t\hat {\boldsymbol{u}}+\textrm{i}\!\left (\left (d_t{\boldsymbol{k}}+\left (\boldsymbol{\nabla }{\boldsymbol{U}}\right )^T{\boldsymbol{\cdot }}\,{\boldsymbol{k}}\right )\!{\boldsymbol{\cdot }}{\boldsymbol{x}}\right )\!\hat {\boldsymbol{u}}\right ]\exp \!\left (\textrm{i}{\boldsymbol{k}}(t){\boldsymbol{\cdot }}{\boldsymbol{x}}\right )\! . \end{equation}

To remove the explicit dependence on $\boldsymbol{x}$ in the resulting equations for the Fourier amplitudes $\hat {\boldsymbol{u}},$ $\hat {\boldsymbol{b}}$ and $\hat \vartheta$ one has to ensure that ${\boldsymbol{k}}(t )$ varies in time according to the eikonal equation (see, Bayly Reference Bayly1986; Craik   Criminale Reference Craik and Criminale1986),

(3.5) \begin{equation} d_t{\boldsymbol{k}}=-\!\left (\boldsymbol{\nabla }{\boldsymbol{U}}\right )^T{\boldsymbol{\cdot }} {\boldsymbol{k}} \end{equation}

where $d_t ({\boldsymbol{\cdot }})\equiv d ({\boldsymbol{\cdot }})/{\rm d}t.$

Substituting the plane waves solution (3.3) into the system (3.1) and taking into account the eikonal (3.5), we obtain

(3.6a) \begin{align} d_t{\hat {\boldsymbol{u}}}&=-\textrm{i} \hat p {\boldsymbol{k}} -\hat {\boldsymbol{u}}{\boldsymbol{\cdot }} \boldsymbol{\nabla }{\boldsymbol{U}}-2\varepsilon \varOmega _0\hat {\boldsymbol{x}}\times {\hat {\boldsymbol{u}}} +\textrm{i} \!\left ({\boldsymbol{k}}{\boldsymbol{\cdot }}{\boldsymbol{B}}_0\right )\!{\hat {\boldsymbol{b}}} +N\hat \vartheta {\boldsymbol{n}}, \end{align}

(3.6b) \begin{align} d_t{\hat {\boldsymbol{b}}}&=\hat {\boldsymbol{b}}{\boldsymbol{\cdot }} \boldsymbol{\nabla }{\boldsymbol{U}}+ \textrm{i} \!\left ({\boldsymbol{k}}{\boldsymbol{\cdot }}{\boldsymbol{B}}_0\right )\!{\hat {\boldsymbol{u}}}, \end{align}

(3.6c) \begin{align} d_t{\hat \vartheta }&=-N^{-1}\hat {\boldsymbol{u}}{\boldsymbol{\cdot }}\boldsymbol{\nabla }\varTheta, \end{align}

together with $ {\boldsymbol{k}}{\boldsymbol{\cdot }} {{\hat {\boldsymbol{u}}}}=0$ and ${\boldsymbol{k}} {\boldsymbol{\cdot }} {{\hat {\boldsymbol{b}}}}=0.$ The use of the latter conditions allows one to eliminate the Fourier amplitude of fluctuating pressure,

(3.7) \begin{equation} \hat p=\textrm{i} k^{-2}\left (\left (\hat {\boldsymbol{u}}{\boldsymbol{\cdot }} \boldsymbol{\nabla }{\boldsymbol{U}}\right )\!{\boldsymbol{\cdot }}\,{\boldsymbol{k}}+\left (\left (\boldsymbol{\nabla }{\boldsymbol{U}}\right )^T{\boldsymbol{\cdot }}\,{\boldsymbol{k}}\right )\!{\boldsymbol{\cdot }}\,\hat {\boldsymbol{u}}+2\varepsilon \varOmega _0 \left (\hat {\boldsymbol{x}}\times \hat {\boldsymbol{u}}\right )\!{\boldsymbol{\cdot }}\,{\boldsymbol{k}}-N\left ({\boldsymbol{n}}{\boldsymbol{\cdot }}\,{\boldsymbol{k}}\right )\!\hat \vartheta \right ) \end{equation}

and to reduce the above seven-component Floquet system to a five-component version, where $k=\Vert {\boldsymbol{k}}\Vert .$ The spectral counterpart of the linear part of MIPS given by (3.2) reads

(3.8) \begin{equation} \hat \varpi _m=\hat {\boldsymbol{b}}{\boldsymbol{\cdot }}\boldsymbol{\nabla }\varTheta +\textrm{i} \!\left ({\boldsymbol{k}}{\boldsymbol{\cdot }}{\boldsymbol{B}}_0\right )\!\hat \vartheta =\text{const.} \end{equation}

It will be used to reduce the five-dimensional Floquet system to a inhomogeneous four-dimensional Floquet system.

3.3. The KBF case

The developments presented in this section as well as those presented in § 3.4 concern the KBF. Those which concern the MBF are reported in Appendix A for greater clarity for the reader.

3.3.1. Time-periodic wavevector

For the KBF (described by (2.18)), (3.5) reduces to

(3.9) \begin{equation} d_tk_1=-\varOmega _0k_2,\quad d_tk_2=\varOmega _0k_1,\quad d_tk_3=2\varepsilon \varOmega _0 k_2=-2\varepsilon d_tk_1, \end{equation}

so that

(3.10) \begin{equation} d_tk^2=d_t\big(k_1^2+k_2^2+k_3^2\big)= 4\varOmega _0\varepsilon k_2k_3, \end{equation}

where $k=\Vert {\boldsymbol{k}}\Vert .$ Equation (3.9) can be solved to give

(3.11) \begin{equation} k_1=k_{\!p}\cos (\tau -\phi ),\quad k_2=k_{\!p}\sin (\tau -\phi ),\quad k_3=k_0-2\varepsilon k_1, \end{equation}

where $\tau =\varOmega _0t$ is a dimensionless time, $\tan \phi =-k_{20}/k_{10}$ and

(3.12) \begin{equation} k_{\!p}=\sqrt {k_1^2+k_2^2}=\sqrt {k_{10}^2+k_{20}^2},\quad k_0=k_{30}+2\varepsilon k_{10}=k_3+2\varepsilon k_1 ,\end{equation}

in which ${\boldsymbol{K}}_0=(k_{10},k_{20},k_{30})^T$ is the initial wavevector. Both $k_{\!p}$ and $k_0$ are time-independent.

For purposes of studying stability, we may set $\phi =0.$ This is easily seen by making the substitution $\varOmega _0t'= \varOmega _0 t - \phi ,$ which eliminates $\phi$ from the equation.

3.3.2. One-dimensional two component and 2-D three component flow cases

For $k_{\!p}=0,$ so that $k_1=k_2=0$ : the wavevector aligns with the vertical unit vector $\hat {\boldsymbol{z}},$ ${\boldsymbol{k}}=k_{3}\hat {\boldsymbol{z}}=k_{30}\hat {\boldsymbol{z}}.$ In this case, the constraints ${\boldsymbol{k}}{\boldsymbol{\cdot }}\,\hat {\boldsymbol{u}}=0$ and ${\boldsymbol{k}}{\boldsymbol{\cdot }}\,\hat {\boldsymbol{b}}=0$ give $\hat u_3=0$ and $\hat b_3=0,$ respectively. Thus, this flow configuration corresponds to one-dimensional two component (1D2C) flows. System (3.6) reduces to a four-dimensional system for the velocity and magnetic field horizontal fluctuations which are not affected by the effect of stratification (at least for the linear regime, i.e. the nonlinear interactions are omitted),

(3.13) \begin{equation} d_t \left (\!\!\begin{array}{c} \hat u_1\\ \hat u_2\\ \hat b_1\\ \hat b_2 \end{array}\!\!\right )= \left (\!\! \begin{array}{cccc} 0&\varOmega _0&\textrm{i} k_3B_0&0\\ -\varOmega _0&0&0&\textrm{i} k_3B_0\\ \textrm{i} k_3B_0&0&0&-\varOmega _0\\ 0&\textrm{i} k_3B_0&\varOmega _0&0\\ \end{array}\!\!\right ) {\boldsymbol{\cdot }} \left (\!\!\begin{array}{c} \hat u_1\\ \hat u_2\\ \hat b_1\\ \hat b_2 \end{array}\!\!\right ) \end{equation}

with $\hat \vartheta =\text{const.},$ in agreement with (3.8). It is not difficult to show that there can be no instability associated with the solution of the above differential system.

For the case where $k_0=0,$ so that $k_3=-2\varepsilon k_1$ and ${\boldsymbol{k}}=k_1 (\hat {\boldsymbol{x}}-2\varepsilon \hat {\boldsymbol{z}} )+k_2\hat {\boldsymbol{y}},$ and the constraints ${\boldsymbol{k}}{\boldsymbol{\cdot }}\,\hat {\boldsymbol{u}}=0$ and ${\boldsymbol{k}}{\boldsymbol{\cdot }}\, \hat {\boldsymbol{b}}=0$ reduce to

(3.14) \begin{equation} k_1\big(\hat u_1+2\varepsilon \hat u_3\big)+k_2\hat u_2=0,\quad k_1\big (\hat b_1+2\varepsilon \hat b_3\big )+k_2\hat b_2=0. \end{equation}

Moreover, in the linear regime, there is no effect of the basic magnetic field on the velocity and density fluctuations. The magnetic field fluctuations behave like a passive field in a precessing stratified flow, since the system (3.1) reduces to

(3.15a) \begin{align} d_t{\hat {\boldsymbol{u}}}&=-\textrm{i} \hat p {\boldsymbol{k}} -\hat {\boldsymbol{u}}{\boldsymbol{\cdot }} \boldsymbol{\nabla }{\boldsymbol{U}}-2\varepsilon \varOmega _0\hat {\boldsymbol{x}}\times {\hat {\boldsymbol{u}}} +N\hat \vartheta {\boldsymbol{n}}, \end{align}

(3.15b) \begin{align} d_t{\hat {\boldsymbol{b}}}&=\hat {\boldsymbol{b}}{\boldsymbol{\cdot }}\boldsymbol{\nabla }{\boldsymbol{U}}, \end{align}

(3.15c) \begin{align} d_t{\hat \vartheta }&=-N^{-1}\hat {\boldsymbol{u}}{\boldsymbol{\cdot }}\boldsymbol{\nabla }\varTheta , \end{align}

or equivalently,

(3.16a) \begin{align} \hat \omega _{3}&=\text{const.},\quad \hat b_3=\text{const.},\quad k_1\hat b_2-k_2\hat b_1 +2\varepsilon k_2\hat b_3= \text{const.}, \end{align}

(3.16b) \begin{align} d_\tau \hat \psi &=-2\textrm{i}\varepsilon k_1\hat \omega _3+\left (N/\varOmega _0\right )k_{\!p}^2\hat \vartheta , \end{align}

(3.16c) \begin{align} d_\tau \vartheta &=-\!\left (N/\varOmega _0\right )\!\hat u_3, \end{align}

where $\hat \psi =k^2\hat u_3,$ $k^2=k_{\!p}^2+k_3^2=k_{\!p}^2+4\varepsilon ^2k_1^2= k_{\!p}^2 (1+4\varepsilon ^2\cos ^2\tau )$ and $\hat \omega _3\equiv \textrm{i}(k_1\hat u_2-k_2\hat u_1)$ denotes the Fourier amplitude of the vertical vorticity component ( ${\boldsymbol{\omega }} =\boldsymbol{\nabla }\times {\boldsymbol{u}})$ . Accordingly, we deduce the following inhomogeneous second-order differential equation (i.e. a Hill’s equation):

(3.17) \begin{equation} d_{\tau \tau }\hat \psi +{{\left (N/\varOmega _0\right )}^2}\big (1+2\varepsilon ^2\!\left (1+2\cos 2\tau \right )\big )\hat \psi =-2\varepsilon k_{\!p}\hat \omega _3\sin \tau . \end{equation}

The stability analysis of (3.17) indicates that there is no instability at order ${ O}(\varepsilon ).$ In counterpart, on times cales of order $\varepsilon ^{-2},$ the solutions are unstable at and near ${N/\varOmega _0}=1$ with a maximal growth rate, $\sigma _{m}={\varepsilon ^2}/{2}$ (more details are given in the Supplementary material).

3.4. Floquet system for the 3-D three component flow case

For perturbations with wavenumber $k_0\ne 0,$ we transform the resulting Floquet system in terms of the following variables to facilitate subsequent calculations (see, Salhi et al. Reference Salhi, Lehner and Cambon2010):

(3.18a) \begin{align} \varOmega _0 \hat c_1&=\left (k_1\hat u_2-k_2\hat u_1\right )+2\varepsilon k_2\hat u_3, \quad \varOmega _0 \hat c_2=-k_0\hat u_3, \end{align}

(3.18b) \begin{align} \varOmega _0\hat c_3&=\big (k_1\hat b_2-k_2\hat b_1\big )+2\varepsilon k_2\hat b_3, \quad \varOmega _0 \hat c_4=-k_0\hat b_3,\quad \varOmega _0 \hat c_5=-k_0\hat \vartheta . \end{align}

Given (3.7), which reduces to

(3.19) \begin{equation} \hat p=-\textrm{i} k^{-2}\big (2\varOmega _0\!\left (k_1\hat u_2-k_2\hat u_1\right )+6\varepsilon \varOmega _0k_2\hat u_3-2\varepsilon \varOmega _0 k_3\hat u_2+k_3N\hat \vartheta \big ), \end{equation}

we transform the system (3.6) according to these new variables,

(3.20a) \begin{align} d_\tau \hat c_1&=-4\varepsilon \left (\frac {k_2k_3}{k^2}+\varepsilon \frac {k_1k_2}{k^2}\right )\!\hat c_1-2\!\left (1+\varepsilon \frac {k_1}{k_0}+4\varepsilon ^3\frac {k_1k_2^2}{k_0k^2}\right )\!\hat c_2+\textrm{i} {\mathcal B}\mu \hat c_3-2\varepsilon {\mathcal N} \frac {k_2k_{\!p}^2}{k_0k^2}\hat c_5, \end{align}

(3.20b) \begin{align} d_\tau \hat c_2&=2\frac {k_0\!\left (k_0-\varepsilon k_1\right )}{k^2}\hat c_1+4\varepsilon ^2\frac {k_1k_2}{k^2}\hat c_2+ \textrm{i} {\mathcal B}\mu \hat c_4+{\mathcal N} \frac {k_{\!p}^2}{k^2}\hat c_5, \end{align}

(3.20c) \begin{align} d_\tau \hat c_3&=\textrm{i} {\mathcal B}\mu \hat c_1, \end{align}

(3.20d) \begin{align} d_\tau \hat c_4&=\textrm{i} {\mathcal B}\mu \hat c_2, \end{align}

(3.20e) \begin{align} d_\tau \hat c_5&=-{\mathcal N}\hat c_2 \end{align}

where $d_\tau ({\boldsymbol{\cdot }})=\varOmega _0^{-1}d_t({\boldsymbol{\cdot }})$ and

(3.21) \begin{equation} -\!1\leqslant \mu =\frac {k_0}{\sqrt {k_{\!p}^2+k_0^2}}\leqslant 1,\quad {\mathcal B}=\frac {B_0\sqrt {k_{\!p}^2+k_0^2}}{\varOmega _0},\quad {\mathcal N}=\frac {N}{\varOmega _0}, \end{equation}

which are the parameters of the Floquet system in addition to the pression parameter $\varepsilon .$ We note that, at small scales, i.e. $\sqrt {k_{\!p}^2+k_0^2}\gg 1,$ $\mathcal B$ is large even if $B_0/\varOmega _0\sim 1.$

On the other hand, combining (3.20d ) and (3.20e ), we deduce the following relation:

(3.22) \begin{equation} {\mathcal N}\hat c_4+\textrm{i} {\mathcal B}\mu \hat c_5= -\varOmega _0^{-1}k_0\!\left (N\hat b_3+\textrm{i} B_0k_0\hat \vartheta \right )=-\varOmega _0^{-1}k_0\hat \varpi = \text{const.,} \end{equation}

which represents the spectral counterpart of the linear part of MIPS (see (3.8)). Accordingly, system (3.20) can be further reduced to a fourth-order inhomogeneous Floquet system for $\hat {\boldsymbol{c}}= (\hat c_1,\hat c_2,\hat c_3,\hat c_4 )^T$ ,

(3.23) \begin{equation} d_\tau \hat {\boldsymbol{c}}=\boldsymbol{D}{\boldsymbol{\cdot }}\, \hat {\boldsymbol{c}} +\hat {\boldsymbol{\varphi }}, \end{equation}

where the only non-zero elements of $\boldsymbol{D}(\tau )$ are

(3.24a) \begin{align} D_{11}&=-4\varepsilon \left (\frac {k_2k_3}{k^2}+\varepsilon \frac {k_1k_2}{k^2}\right )\! , \end{align}

(3.24b) \begin{align} D_{12}&=-2\!\left (1+\varepsilon \frac {k_1}{k_0}+4\varepsilon ^3\frac {k_1k_2^2}{k_0k^2}\right )\! , \end{align}

(3.24c) \begin{align} D_{13}&=D_{31}=D_{42}=\textrm{i} {\mathcal B}\mu , \end{align}

(3.24d) \begin{align} D_{14}&=-2\textrm{i}\varepsilon \frac {{\mathcal N}}{{\mathcal B}\mu ^2}\frac {k_2k_{\!p}^2}{k^3}, \end{align}

(3.24e) \begin{align} D_{21}&=2\frac {k_0\!\left (k_0-\varepsilon k_1\right )}{k^2}, \end{align}

(3.24f) \begin{align} D_{22}&= 4\varepsilon ^2\frac {k_1k_2}{k^2}, \end{align}

(3.24g) \begin{align} D_{24}&=\textrm{i} {\mathcal B}\mu +\textrm{i} \frac {{\mathcal N}}{{\mathcal B}\mu }\frac {k_{\!p}^2}{k^2}. \end{align}

Therefore, the homogeneous system is

(3.25) \begin{equation} d_\tau \hat {\boldsymbol{c}}=\boldsymbol{D}{\boldsymbol{\cdot }}\,\hat {\boldsymbol{c}}, \end{equation}

while the non-zero components of the inhomogeneous term in (3.23) take the form

(3.26) \begin{equation} \hat \varphi _1=-2\textrm{i}\varepsilon \frac {N\varOmega _0}{B_0} \frac {k_2k_{\!p}^2}{k_0k^2}\hat \varpi _m,\quad \hat \varphi _2= \textrm{i}\frac {N\varOmega _0}{B_0} \frac {k_{\!p}^2}{k^2}\hat \varpi _m, \end{equation}

and it can be seen as a time-varying forcing excitation. The linear system (3.23) has the properties $\boldsymbol{D}(\tau +T) =\boldsymbol{D}(\tau )$ and $\hat {\boldsymbol{\varphi }}(\tau +T)=\hat {\boldsymbol{\varphi }}(\tau )$ where $T = 2\pi$ is the period common to both the matrix $\boldsymbol{D}(\tau )$ and the vector $\hat {\boldsymbol{\varphi }}(\tau ).$ Floquet theory does not address stability of the inhomogeneous system described by (3.20) where the ‘forcing excitation’ $\hat {\boldsymbol{\varphi }}(\tau )$ is present. However, the $T\hbox{-}$ periodic nature of $\hat {\boldsymbol{\varphi }}(\tau )$ allows an extension to the theory (see, Slane   Tragesser Reference Slane and Tragesser2011). By adding a forcing excitation, the general behaviour predicted by Floquet theory for the homogeneous system changed only for the case of semisimple Floquet multipliers that are identically equal to one, i.e. one is the only multiplier on the unit circle and is semisimple. In this case, the presence of the forcing excitation caused the solution to diverge (see, Slane   Tragesser Reference Slane and Tragesser2011). Incidentally, for several values of the parameters $({\mathcal B},{\mathcal N},\varepsilon , \mu ),$ we have numerically determined the stability and instability regions by integrating the $5\times 5$ Floquet system (without involving MIPS, i.e. system (3.20)) as well as the resulting $4\times 4$ Floquet system with $\hat \varpi _m=0$ (i.e. system (3.25)). We obtained the same results for both systems (up to numerical errors). We note that the determinant of the fundamental matrix solution is unity at $\tau =2\pi$ for both Floquet systems (see (4.3)). This indicates that the neutral case of system (3.25) could be characterized by the fact that unity is not the only multiplier on the unit circle.

3.5. Resonant cases of MIG waves

In this section, we start from the case without precession, so that both the background shear and the Coriolis force vanish. In that case, there are MIG waves. The resonant cases of MIG waves have been characterized by SC23 in their study of the MGE instability . Because some resonant cases of MIG waves become destabilizing in stratified magnetized precessing flows $(\varepsilon \ne 0),$ we briefly recall here the analysis made by SC23.

We denote by $\boldsymbol{D}_0$ the matrix $\boldsymbol{D}$ for $\varepsilon =0$ (i.e. circular streamlines)

(3.27) \begin{equation} \boldsymbol{D}_0=\varOmega _0\!\left (\!\!\begin{array}{cccc} 0&-2\varOmega _0&\textrm{i} \omega _a&0\\ (2\varOmega _0)^{-1}\omega _r^2&0&0&\textrm{i} \big (\omega _a^2+\omega _g^2\big )\omega _a^{-1}\\ \textrm{i} \omega _a&0&0&0\\ 0&\textrm{i} \omega _a&0&0 \end{array}\!\!\right ) \end{equation}

where $\omega _r,$ $\omega_a$ and $\omega _g$ are the frequencies of inertial, Alfvén and gravity waves, respectively,

(3.28a) \begin{align} \omega _r&=\frac {2\vert {\boldsymbol{\varOmega }}_0{\boldsymbol{\cdot }} {\boldsymbol{K}}_0\vert }{K_0}=2\varOmega _0\vert \mu \vert , \end{align}

(3.28b) \begin{align} \omega _a&={\boldsymbol{B}}_0{\boldsymbol{\cdot }} {\boldsymbol{K}}_0=B_0k_0=B_0 K_0\vert \mu \vert , \end{align}

(3.28c) \begin{align} \omega _g&=N\frac {\Vert \hat {\boldsymbol{z}}\times {\boldsymbol{K}}_0\Vert }{K_0}=N\sqrt {1-\mu ^2}. \end{align}

The eigenvalues $\sigma _j$ $( j = 1, 2, 3, 4)$ of the constant matrix $\boldsymbol{D}_0$ such that

(3.29) \begin{equation} \left \vert \boldsymbol{D}_0-\sigma \boldsymbol{I}_4\right \vert =\left \vert \boldsymbol{D}_0-\textrm{i} \omega \boldsymbol{I}_4\right \vert =0, \end{equation}

where $\boldsymbol{I}_4$ is the four-dimensional unit matrix, take the form

(3.30a) \begin{align} -\varOmega _0^2\sigma _{1,2}^2&=\omega ^2_{1,2}=\frac {1}{2}\left [2\omega _a^2+\omega _r^2+\omega _g^2+\sqrt {\left (\omega _r^2+\omega _g^2\right )^2+4\omega _r^2\omega _a^2}\ \right ]\!, \end{align}

(3.30b) \begin{align} -\varOmega _0^2\sigma _{3,4}^2&=\omega _{3,4}^2=\frac {1}{2}\left [2\omega _a^2+\omega _r^2+\omega _g^2-\sqrt {\left (\omega _r^2+\omega _g^2\right )^2+4\omega _r^2\omega _a^2}\ \right ]\!, \end{align}

or equivalently,

(3.31a) \begin{align} \omega _1=-\omega _{2}&=\frac {\varOmega _0}{\sqrt {2}}\sqrt {\left (4+2{\mathcal B}^2-{\mathcal N}^2\right )\mu ^2+{\mathcal N}^2+\sqrt {\left [\left (4-{\mathcal N}^2\right )\mu ^2+{\mathcal N}^2\right ]^2+16{\mathcal B}^2\mu ^4}}, \end{align}

(3.31b) \begin{align} \omega _3=-\omega _{4}&=\frac {\varOmega _0}{\sqrt {2}}\sqrt {\left (4+2{\mathcal B}^2-{\mathcal N}^2\right )\mu ^2+{\mathcal N}^2-\sqrt {\left [\left (4-{\mathcal N}^2\right )\mu ^2+{\mathcal N}^2\right ]^2+16{\mathcal B}^2\mu ^4}}. \end{align}

Because $\omega _1\geqslant \omega _3,$ fast waves have frequency $\omega _1,$ while slow waves have frequency $\omega _3.$ In the non-magnetized case, fast modes reduce to the purely hydrodynamic modes (i.e. inertia-gravity waves). In counterpart, slow waves are purely magnetic modes because $\omega _3$ reduces to zero for ${\mathcal B}=0$ (see LZ04). The expression of the eigenvectors associated with the eigenvalues $\sigma _1,\sigma _2,\sigma _3$ and $\sigma _4$ is reported in Appendix B.

The resonant cases of MIG waves are those parameter values $(\mu , {\mathcal N},{\mathcal B},)$ such that

(3.32) \begin{equation} \omega _i-\omega _j= n \varOmega _0,\quad (i\ne j). \end{equation}

As will be shown later, some resonant cases of order $n=1$ are destabilizing, whereas in the case of an elliptical flow the only resonant cases that can lead to instability are those for which the integer $n$ is even (see SC23). Indeed, in the precessing flow case, the components of the wavevector vary in time as $\sim \exp (\pm \textrm{i} \varOmega _{0} t)$ , reflecting the $\varOmega _{0}$ -frequency of the strain field, while in the elliptical flow case they vary as $\sim \exp (\pm 2 \textrm{i} \varOmega _{0} t)$ , corresponding to a doubled frequency. This difference directly accounts for the fact that precessing flows can exhibit destabilizing resonances already for $n = 1$ , whereas elliptical flows only admit even- $n$ resonances.

Note that the condition (3.32) for resonance readily extends the one by Bayly (Reference Bayly1986) for basic elliptical flow instability ( ${\mathcal B}=0$ and ${\mathcal N}=0),$ so that, $4 \varOmega \mu = n \varOmega _0,$ that immediately yielded $\mu = n/4=0.5$ for the elliptical flow ( $n=2)$ and $n/4=0.25$ for the precessing flow ( $n=1)$ (see, Kerswell Reference Kerswell1993; Salhi   Cambon Reference Salhi and Cambon2009) giving the origin of time-dependent instability tongues at vanishing $\varepsilon$ . We also note that the analysis of triad instability principle by Waleffe (Reference Waleffe1992) shows that the elliptical instability corresponds to a forward-interaction: the two modes with eigenfrequency $\omega _1=\omega _r$ and $\omega _2=-\omega _r$ have opposite polarities and are coupled with the mean flow, which is associated with a zero frequency, the unstable modes are thereby two resonant inertial waves associated with the uniform background rotation.

4.1. Expansion in Taylor series of the Floquet multiplier matrix

We denote by ${\boldsymbol \varPhi }(\tau ,\varepsilon ,\mu ,{\mathcal B}, {\mathcal N})$ the fundamental matrix solution,

(4.1) \begin{equation} d_\tau {\boldsymbol \varPhi }=\boldsymbol{D}{\boldsymbol{\cdot }} {\boldsymbol \varPhi },\quad {\boldsymbol \varPhi }(\tau =0)=\boldsymbol{I}_4 \end{equation}

and by $\boldsymbol{M}= {\boldsymbol \varPhi }(2\pi ,\varepsilon ,\mu ,{\mathcal B},{\mathcal N})$ the Floquet multiplier matrix. According to Floquet–Lyapunov theorem, $\boldsymbol \varPhi$ is expressible in the form (see, Kuchment Reference Kuchment1993),

(4.2) \begin{equation} \boldsymbol \varPhi (\tau )=\boldsymbol{F}\!(\tau )\exp \!\left (\boldsymbol{K}\!\tau \right )\! , \end{equation}

where $\boldsymbol{F}\!(\tau )$ is a non-singular continuous $2\pi \hbox{-}$ periodic $4\times 4$ matrix-function (whose derivative is an integrable piecewise-continuous function) and $\boldsymbol{K}$ is a constant matrix. The determinant of $\boldsymbol \varPhi$ is unity at $\tau =2\pi ,$ $\vert {\boldsymbol \varPhi (2\pi )}\vert =1$ since

(4.3) \begin{equation} \text{trace}\ \boldsymbol{D}=\sum _{j=1}^4 D_{\textit{jj}}=-4\varepsilon \frac {k_2k_3}{k^2}=-\frac {1}{k^2}{d_\tau k^2}. \end{equation}

It follows that whenever $\lambda$ is an eigenvalue of the monodromy matrix, $\boldsymbol{M}=\boldsymbol \varPhi (2\pi ),$ so also are its inverse $\lambda ^{-1}$ and its complex conjugate $\overline {\lambda }$ (see LZ04). Consequently, in the stable case, eigenvalues of $\boldsymbol{M}$ lie on the unit circle. If any eigenvalue $\lambda$ of $\boldsymbol{M}$ has modulus exceeding one, this implies that there is indeed an exponentially growing solution. The growth rates are then given by

(4.4) \begin{equation} \sigma =\frac {1}{2\pi }\log \left (\lambda \right )\! . \end{equation}

It follows that whenever $\lambda$ is an eigenvalue of $\boldsymbol{M},$ so also are its inverse $\lambda ^{-1}$ and its complex conjugate $\overline {\lambda }.$ We denote by

(4.5) \begin{equation} {p}(\lambda ,\varepsilon )=\vert \boldsymbol{M}-\lambda \boldsymbol{I}_{\!4}\vert \end{equation}

the characteristic polynomial of the Floquet multiplier matrix $\boldsymbol{M}$ and by $\varLambda _1,\varLambda _2,\varLambda _3$ and $\varLambda _4$ its roots. A necessary condition for stability is that each root lie on the unit circle (see LZ04).

We expand the Floquet multiplier matrix $\boldsymbol{M}\!(\varepsilon ,\mu , {\mathcal B}, {\mathcal N})$ in Taylor series in the neighbourhood of $(\varepsilon ,\mu )=(0,\mu _0),$ holding $({\mathcal B},{\mathcal N})$ constant,

(4.6) \begin{equation} \boldsymbol{M}=\boldsymbol{M}_0(\mu _0,0)+\varepsilon \boldsymbol{M}_\varepsilon (\mu _0,0)+ (\mu -\mu _0)\boldsymbol{M}_{\!\mu} (\mu _0,0)+{\cdots} \end{equation}

where $\boldsymbol{M}_\varepsilon = (\partial \boldsymbol{M}/\partial \varepsilon ),$ $\boldsymbol{M}_{\!\mu} = (\partial \boldsymbol{M}/\partial \mu )$ and the dots indicate higher-order terms in $\varepsilon$ and $\mu -\mu _0.$ Generally, at sufficiently small $\varepsilon ,$ the region in the $(\varepsilon , \mu )\hbox{-}$ plane where instability occurs is typically a wedge with apex at a point $(\varepsilon , \mu )) = (0, \mu _0)$ and boundaries

(4.7) \begin{equation} \mu =\mu _0+{\eta _\pm }\varepsilon, \end{equation}

where the slopes $\eta _+$ and $\eta _{-}$ are to be found. The instability (if it exists) has a bandwidth $ (\eta _+-\eta _- )\varepsilon$ that is, for given $\varepsilon ,$ $\mathcal B$ and ${\mathcal N},$ the length of the $\mu \hbox{-}$ interval for which the wavenumbers are unstable. Therefore, (4.6) can be rewritten as

(4.8a) \begin{align} \boldsymbol{M}&=\boldsymbol{M}_0+\varepsilon \boldsymbol{M}_1+{ O}(\varepsilon ^2), \end{align}

(4.8b) \begin{align} \boldsymbol{M}_1&=\boldsymbol{M}_\varepsilon +{\eta } \boldsymbol{M}\mu . \end{align}

Accordingly, we no longer need the designation $\mu _0$ and, hereinafter, use the symbol $\mu$ in its place.

To determine matrices $\boldsymbol{M}_0,$ $\boldsymbol{M}_\varepsilon$ and $\boldsymbol{M}_{\!\mu} ,$ we expand, for a given $\tau \in [0, 2\pi ],$ $\boldsymbol \varPhi$ and $\boldsymbol{D},$

(4.9a) \begin{align} {\boldsymbol \varPhi }(\tau ,\varepsilon )&={\boldsymbol \varPhi }_0(\tau ,\mu ,0)+\varepsilon {\boldsymbol \varPhi }_1(\tau ,\mu ,0)+{ O}(\varepsilon ^2), \end{align}

(4.9b) \begin{align} \boldsymbol{D}(\tau ,\varepsilon )&=\boldsymbol{D}_0+\varepsilon \boldsymbol{D}_\varepsilon (\tau ,0)+{ O}(\varepsilon ^2), \end{align}

where ${\boldsymbol \varPhi }_0(\tau =0)=\boldsymbol{I}_4$ and ${\boldsymbol \varPhi }_1(\tau =0)=\boldsymbol{0}.$ We notice that we find the same expression for the matrix $\boldsymbol{D}_\varepsilon (\tau , \varepsilon )$ for both KBF and MBF cases (see Appendix B2). Thus, to leading order in $\varepsilon ,$ the stability problem for these two basic flows is the same.

Substituting (4.9) into (4.1), we obtain

(4.10) \begin{align} d_\tau {\boldsymbol \varPhi }_0=\boldsymbol{D}_0{\boldsymbol{\cdot }} {\boldsymbol \varPhi }_0,\quad d_\tau {\boldsymbol \varPhi }_1=\boldsymbol{D}_0{\boldsymbol{\cdot }} {\boldsymbol \varPhi }_1+\boldsymbol{D}_\varepsilon {\boldsymbol{\cdot }} {\boldsymbol \varPhi }_0, \end{align}

with solution ${\boldsymbol \varPhi }_0(\tau )=\text{e}^{\tau \boldsymbol{D}_0}$ and

(4.11) \begin{equation} {\boldsymbol \varPhi }_1(\tau )={\boldsymbol \varPhi }_0(\tau ){\boldsymbol{\cdot }}\left (\int _0^\tau {\boldsymbol \varPhi }_0^{-1}(s){\boldsymbol{\cdot }}\boldsymbol{D}_\varepsilon (s) {\boldsymbol{\cdot }} {\boldsymbol \varPhi }_0(s){\rm d}s\right )\! . \end{equation}

Because the characteristic polynomial $\text{p}(\lambda ,\varepsilon )$ is the same in any coordinate system and the four eigenvalues of the matrix $\boldsymbol{D}_0,$ given by (3.30) are distinct as long as ${\mathcal B}\ne 0$ and $0\lt \mu ^2\lt 1,$ we transform the solution in the basis diagonalizing $\boldsymbol{D}_0,$ so that

(4.12a) \begin{align} \tilde {\boldsymbol{M}}_0&={\boldsymbol{T}}^{-1}{\boldsymbol{\cdot }} {\boldsymbol{M}}_0{\boldsymbol{\cdot }} \boldsymbol{T}=\textbf{diag}\big ({\text{e}}^{2\pi \sigma _1}, {\text{e}}^{-2\pi \sigma _1}, {\text{e}}^{2\pi \sigma _3}, \text{e}^{-2\pi \sigma _3}\big )= \textbf{diag}\left (\lambda _1, \lambda _2, \lambda _3, \lambda _4\right )\! , \end{align}

(4.12b) \begin{align} \tilde {\boldsymbol{M}}_\varepsilon &={\boldsymbol{T}}^{-1}{\boldsymbol{\cdot }} {\boldsymbol{M}}_\varepsilon {\boldsymbol{\cdot }} \boldsymbol{T}=\tilde {\boldsymbol{M}}_0{\boldsymbol{\cdot }} \tilde {\boldsymbol{J}}, \end{align}

(4.12c) \begin{align} \tilde {J}_{\textit{ij}}&=\left ({\boldsymbol{T}}^{-1}\right )_{im}{ T}_{lj}\int _0^{2\pi } {\text{e}}^{(\sigma _j-\sigma _i)\tau }\left ({\boldsymbol{D}}_\varepsilon \right )_{ml}(\tau ){\rm d}\tau , \end{align}

where the columns of the matrix $\boldsymbol{T}$ are the eigenvectors associated with the eigenvalues $\sigma _j$ ( $j=1,2,3,4)$ given by (B1) and $\boldsymbol{T}^{-1}$ is its inverse given by (B2).

To complete the construction of the matrix $\tilde {\boldsymbol{M}}_1,$ which appears in (4.8), we need the derivative of $\tilde {\boldsymbol{M}}_0(\mu )=\tilde {\boldsymbol{M}}(\mu ,0)$ with respect to $\mu ,$

(4.13) \begin{equation} {\eta }\tilde {\boldsymbol{M}}_\mu ={\eta }\frac {\partial \tilde {\boldsymbol{M}}_0}{\partial \mu }=\frac {2\textrm{i}\pi {\eta }}{\varOmega }\ \textbf{diag}\left ( \frac {\partial \omega _1}{\partial \mu }\lambda _1, \frac {\partial \omega _2}{\partial \mu }\lambda _2, \frac {\partial \omega _3}{\partial \mu }\lambda _3, \frac {\partial \omega _4}{\partial \mu }\lambda _4 \right )\! , \end{equation}

where $\omega _1,$ $\omega _2,$ $\omega _3$ and $\omega _4$ are given by (3.31).

When there is a subharmonic instability associated with the resonant case $\omega _\ell -\omega _j=\pm 1$ ( $j\ne \ell ),$ its maximum growth rate $\sigma _{m}$ is given by

(4.14) \begin{equation} \sigma _{m}=\frac {\varepsilon }{2\pi }\Re \sqrt {\tilde J_{\ell j}\tilde J_{j\ell }}\quad (\text{with}\ \ell \ne j), \end{equation}

recalling that the frequencies $\omega _j$ ( $j=1,2,3,4)$ are normalized by $\varOmega _0.$ In the following sections, we present the analytical results of the maximum growth rates associated with the three unstable cases, namely the fast–fast, slow–slow and fast–slow destabilizing resonances (calculation details are given in the Supplementary material).

4.2. Case of unstable resonance between two fast modes or between two slow modes

4.2.1. Domains of the resonant cases in the $({\mathcal B},{\mathcal N})\hbox{-}$ plane

The resonance condition between two fast modes is obtained by substituting the expression of $\omega _1$ given by (3.31) into the following identity:

(4.15) \begin{equation} \omega _1-\omega _2=2\omega _1=1\ \implies \ 16\omega _1^4=1, \end{equation}

or equivalently,

(4.16) \begin{equation} 4{\mathcal B}^2\big ({\mathcal B}^2-{\mathcal N}^2\big)\mu ^4+\big (4{\mathcal B}^2{\mathcal N}^2-2{\mathcal B}^2-4+{\mathcal N}^2\big )\mu ^2+\left (\frac {1}{4}-{\mathcal N}^2\right )=0, \end{equation}

with $0\lt \mu ^2\lt 1.$ Because the replacement $\mu \rightarrow -\mu$ in (4.16) result in the same condition, we may therefore assume without loss of generality that $\mu \gt 0.$ We note that the resonant condition between two slow modes, $\omega _3-\omega _4=1,$ is also described by (4.16), and the real roots of (4.16), such that $0\lt \mu ^2\lt 1,$ if they exist, are ordered. The smaller of the two roots characterizes the resonant cases between two fast modes where

(4.17) \begin{equation} \left ({\mathcal B},{\mathcal N}\right )\in \left [0,+\infty \right )\times \left [0,{1}/{2}\right ]\cup {\mathcal D}_f,\quad \lim _{\substack {{\mathcal B} \to +\infty \\ 0\leqslant {\mathcal N}\leqslant \frac {1}{2}}} {\mathcal B}^2\mu ^2=\frac {1}{4}-{\mathcal N}^2,\quad \lim _{\substack {{\mathcal N} \to +\infty \\ \frac {1}{2}\lt {\mathcal B}\leqslant \frac {\sqrt {5}}{2}}} \mu ^2=\frac {1}{4{\mathcal B}^2}. \end{equation}

The other root of (4.16), such that $0\lt \mu ^2\lt 1,$ characterizes the resonant cases between two slow modes where

(4.18) \begin{equation} \left ({\mathcal B},{\mathcal N}\right )\in \big [{\sqrt {5}}/{2},+\infty \big )\times \left [0,+\infty \right )\cup {\mathcal D}_f,\quad \lim _{\substack {{\mathcal B} \to +\infty \\ 0\leqslant {\mathcal N}\lt +\infty }} {\mathcal B}^2\mu ^2=\frac {1}{4},\quad \lim _{\substack {{\mathcal N} \to +\infty \\ \frac {1}{2}\lt {\mathcal B}\lt +\infty }} \mu ^2=1. \end{equation}

Here, the domain ${\mathcal D}_f$ is defined as follows (see also figure 2):

(4.19) \begin{equation} \forall \left ({\mathcal B},{\mathcal N}\right )\in {\mathcal D}_f \Leftrightarrow \ \text{for given}\ {\mathcal N}\in \big ({{\sqrt {6}}/{2}},+\infty \big )\quad \implies \quad \frac {1}{2}\lt f({\mathcal N})\leqslant {\mathcal B}\lt \frac {\sqrt {5}}{2}, \end{equation}

where $f : (\sqrt {6}/2, +\infty )\rightarrow (1/2,\sqrt {5}/2 )$ is a continuous decreasing function

(4.20) \begin{align} \forall {\mathcal N}\!\in \big ({\sqrt {6}}/{2},+\infty \big ), f({\mathcal N})=\frac {1}{2{\mathcal N}^2} \sqrt {\left ({\mathcal N}^2+2\right )^2-6+2\sqrt {(4{\mathcal N}^2-1)({\mathcal N}^4-1)}}\!\underset {{\mathcal N}\rightarrow +\infty }{\longrightarrow } \!\frac {1}{2}. \end{align}

Figure 2. Resonant cases of MIG waves. (a) Resonant cases (of order $n=1$ ) occur for $({\mathcal B},{\mathcal N})$ belonging to the coloured domains. (b) Resonance between two fast modes ( $\omega _1-\omega _2=\varOmega _0)$ where (I) $\equiv [0,+\infty )\times [0,1]$ and (II) $\equiv {\mathcal D}_f$ (see (4.19)). (c) Resonance between two slow modes ( $\omega _3-\omega _4=\varOmega _0)$ where (III) $\equiv [\sqrt {5}/2,+\infty )\times [0,+\infty ).$ (d) Resonance between a fast mode and a slow mode ( $\omega _1-\omega _4=\varOmega _0$ or $\omega _1-\omega _3=\varOmega _0)$ where (IV) $\equiv [0,+\infty )\times [0,1]$ and (V) $\equiv {\mathcal D}_m$ (see (4.30)).

For the pure hydrodynamic case ( ${\mathcal B}=0,{\mathcal N}=0),$ one has $\mu = {1}/{4}.$ For the resonant cases between to fast modes and for given ${{\mathcal N}\in [0,{1}/{2}]},$ the parameter $\mu$ decreases from $\mu =(\sqrt {1-4{\mathcal N}^2})/(2\sqrt {4-{\mathcal N}^2})$ as $\mathcal B$ increases, approaching zero when ${\mathcal B}\rightarrow +\infty .$ Likewise, for given ${1}/{2}\lt {\mathcal B}={\mathcal B}_f\leqslant {\sqrt {5}}/{2}$ such that $ ({\mathcal N},{\mathcal B} )\in {\mathcal D}_f,$ $\mu$ increases from $\mu (f^{-1}({\mathcal B}_f),{\mathcal B}_f )$ approaching $1$ when ${\mathcal N}\rightarrow +\infty$ (see also figure 3 a).

Figure 3. Resonant cases of MIG waves. (a) Variation of the parameter $\mu =k_0/K_0$ versus ${\mathcal N}=N/\varOmega _0$ for ${\mathcal B}= B_0K_0/\varOmega _0=0.75,\ 1.0,\ \sqrt {5}/2$ such that $({\mathcal B},{\mathcal N})\in {\mathcal D}_f$ (see (4.19)) associated with the resonant cases between two fast modes ( $\omega _1-\omega _2=\varOmega _0)$ or between two slow modes ( $\omega _3-\omega _4=\varOmega _0).$ (b) Variation of the parameter $\mu =k_0/K_0$ versus $\mathcal B$ for ${\mathcal N}=1.2,\ 2.0,\ 3.0$ such that $({\mathcal B},{\mathcal N})\in {\mathcal D}_m$ (see (4.30)) associated with the resonant cases between a fast mode and a slow mode ( $\omega _1-\omega _4=\varOmega _0$ or $\omega _1-\omega _3=\varOmega _0).$ Here, $K_0=\sqrt {k_0^2+k_{\!p}^2}.$

For the resonant cases between two slow modes, which exist only in the magnetized case, and for given ${\mathcal B}\in [{\sqrt {5}}/{2},+\infty ),$ $\mu$ decreases from $\mu =1/(2(\sqrt {1+{\mathcal B}^2}-1))$ as $\mathcal N$ increases approaching $1/(2{\mathcal B})$ when ${\mathcal N}\rightarrow +\infty .$ Likewise, for given ${1}/{2}\lt {\mathcal B}\leqslant {\sqrt {5}}/{2},$ so that $ ({\mathcal N},{\mathcal B} )\in {\mathcal D}_f,$ $\mu$ decreases from $\mu (f^{-1}({\mathcal B}),{\mathcal B} )$ approaching $1/(2{\mathcal B})$ when ${\mathcal N}\rightarrow +\infty$ (see also figure 3 a).

4.2.2. Maximal growth rate

If a subharmonic instability resulting from the resonance between two fast (respectively, slow) modes exists, its maximal growth rate, denoted by $\sigma _{\textit{mf}}$ ( $\sigma _{\textit{ms}}),$ is given by

(4.21) \begin{equation} \sigma _{\textit{mf}}=\frac {\varepsilon }{2\pi }\Re \sqrt {\tilde J_{12}\tilde J_{21}}, \quad \sigma _{\textit{ms}}=\frac {\varepsilon }{2\pi }\Re \sqrt {\tilde J_{34}\tilde J_{43}}. \end{equation}

After some analytical developments given in the Supplementary material, we find

(4.22) \begin{align}& \frac {\sigma _{\textit{mf}}}{\varepsilon },\ \frac {\sigma _{\textit{ms}}}{\varepsilon }=\left \vert \frac {\left (4{\mathcal B}^2\mu ^2+1\right )\mu \sqrt {1-\mu ^2}} {8\left [2\!\left (4+2{\mathcal B}^2-{\mathcal N}^2\right )\mu ^2+2{\mathcal N}^2-1\right ]}\right \vert \nonumber \\&\quad \times \big \vert \big (\big (4{\mathcal B}^2\mu ^2-1\big )\big ({\mathcal B}^2\big (1-4{\mathcal N}^2-4\big ({\mathcal B}^2-{\mathcal N}^2\big )\mu ^2\big )+4\big )+2\big (4{\mathcal B}^2\mu ^2-3\big )\big )\big \vert , \end{align}

in which $0\lt \mu ^2\lt 1$ is the smaller of the two roots of (4.16) when it is a resonance between two fast modes or the other root when it is a resonance between two slow modes.

Some results reported in previous studies of stratified or magnetized precessing flows can be recovered from (4.22),

(4.23a) \begin{align} \text{for}&\ {\mathcal B}={\mathcal N}=0,\quad \mu =\frac {1}{4},\quad \frac {\sigma _{\textit{mf}}}{\varepsilon }=\frac {5\sqrt {15}}{32}, \end{align}

(4.23b) \begin{align} \text{for}&\ {\mathcal B}=0,\quad \mu =\frac {1}{2}\sqrt {\frac {1-4{\mathcal N}^2}{4-{\mathcal N}^2}},\quad \frac {\sigma _{\textit{mf}}}{\varepsilon }=\frac {5\sqrt {15}}{8}\frac {\sqrt {1-4{\mathcal N}^2}}{\left (4-{\mathcal N}^2\right )}, \end{align}

(4.23c) \begin{align} \text{for}&\ {\mathcal N}=0,\ \mu =\frac {1}{2\big (\sqrt {1+{\mathcal B}^2}+1\big ) },\ \frac {\sigma _{\textit{mf}}}{\varepsilon }= \frac {1}{8}\frac {\big (2\sqrt {{\mathcal B}^2+1}+3\big )\sqrt {4 \big (1+\sqrt {1+{\mathcal B}^2}\big )^2\!-1}}{\big (1+\sqrt {{\mathcal B}^2+1}\big )^2}. \end{align}

Equation (4.23b ) indicates that the maximal growth rate is zero at ${\mathcal N}={1}/{2}.$ Therefore, in the non-magnetized stratified precessing case ( ${\mathcal B}=0),$ the subharmonic instability is completely suppressed by stratification when $\mathcal N$ reaches $ {1}/{2}$ (see, Benkacem et al. Reference Benkacem, Salhi, Khlifi, Nasraoui and Cambon2022). For the non-stratified magnetized precessing case ( ${\mathcal N}=0),$ (4.23c ) indicates that $\sigma _{\textit{mf}}/\varepsilon$ decreases as $\mathcal B$ increases so as

(4.24) \begin{equation} \lim _{{\mathcal N}=0,\ {\mathcal B}\rightarrow +\infty } \frac {\sigma _{\textit{mf}}}{\varepsilon }=\frac {1}{2}. \end{equation}

Recall that, at small scales, i.e. $\sqrt {k_{\!p}^2+k_0^2}\gg 1,$ $\mathcal B$ is large even if $B_0/\varOmega _0\sim 1.$

In the stratified magnetized precessing case, the analysis of (4.22) is more subtle as shown in the following.

For $0\lt {\mathcal N}\lt 0.5,$ we verify using the second identity in (4.17) that $\sigma _{\textit{mf}}$ approaches zero for large ${\mathcal B}.$ Given (4.24), this clearly proves that the ${\mathcal N}\rightarrow 0$ limit is, in fact, singular (discontinuous).

For ${\mathcal N}=0.5,$ (4.16) gives $\mu =0$ (i.e. the smaller of the two roots) and hence, $\sigma _{\textit{mf}}=0,$ signifying that there is no parametric instability associated with the resonant case of order $n=1$ between two fast modes. It is worth mentioning that the case $\mu =0,$ so that $k_0=0$ corresponds to a 2-D three component (2D3C) flow (see the end of § 3.3.2).

For $ {1}/{2}\lt {\mathcal B}\leqslant {\sqrt {5}}/{2},$ i.e. $({\mathcal B},{\mathcal N})\in {\mathcal D}_f,$ we verify using the third identity in (4.17) that $\sigma _{\textit{mf}}$ also approaches zero for large ${\mathcal N}.$

Figure 4(a) shows $\sigma _{\textit{mf}}/\varepsilon$ versus $\mathcal B$ for ${\mathcal N}=0,\ 0.2,\ 0.3,\ 0.35$ and ${\mathcal N}=0.45,$ so that $({\mathcal B},{\mathcal N})\in [0,+\infty )\times [0,1/2],$ while figure 4(b) illustrates the variation of $\sigma _{\textit{mf}}/\varepsilon$ versus $\mathcal B$ for ${\mathcal N}=2,$ so that $({\mathcal B},{\mathcal N})\in {\mathcal D}_f.$ The numerical results for $\varepsilon =0.05,$ which will be discussed in § 5, are also reported .

Figure 4. Maximal growth rate of destabilizing resonant cases given by the asymptotic formulae. (a) Resonant cases between two fast modes for some values of $({\mathcal B},{\mathcal N})\in [0,+\infty )\times [0,1/2].$ (b) Resonant cases between two fast modes or between two slow modes for some values of $({\mathcal B},{\mathcal N})\in {\mathcal D}_f$ (see (4.19)). (c) Resonant cases between two slow modes for some values of $({\mathcal B},{\mathcal N})\in [\sqrt {5}/2,+\infty )\times [0,+\infty ).$ (d) Resonant cases between a fast mode and a slow mode for some values of $({\mathcal B},{\mathcal N})\in [0,+\infty )\times [0,1].$ The numerical results obtained for $\varepsilon =0.05$ (symbol) are also reported.

For the non-stratified magnetized case ( ${\mathcal N}=0),$ the expression of the maximal growth rate $\sigma _{\textit{ms}}$ of the destabilizing resonance between two slow mode deduced from (4.22) is of the form

(4.25) \begin{equation} \frac {\sigma _{\textit{ms}}}{\varepsilon }= \frac {\big(2+2{\mathcal B}^2-3\sqrt {{\mathcal B}^2+1}\big)\sqrt {4 \big(\sqrt {1+{\mathcal B}^2}-1\big)^2-1}}{8\sqrt {1+{\mathcal B}^2}\big(\sqrt {{\mathcal B}^2+1}-1\big)^2}\underset {{\mathcal B}\rightarrow +\infty }{\longrightarrow }\frac {1}{2} \end{equation}

where ${\sqrt {5}}/{2}\lt {\mathcal B}\lt +\infty ,$ in agreement with the previous results by Salhi et al. (Reference Salhi, Lehner and Cambon2010). In this case, ${\sigma _{\textit{ms}}}/{\varepsilon }$ increases from $0$ at ${\mathcal B}={\sqrt {5}}/{2}$ to ${1}/{2}$ as ${\mathcal B}\rightarrow +\infty .$

For the stratified magnetized case such that $ ({\mathcal B},{\mathcal N} ) \in [ {\sqrt {5}}/{2},+\infty )\times [0,+\infty )\cup {\mathcal D}_f,$ we use the second and third identities in (4.18) to deduce from (4.22) that

(4.26) \begin{equation} \lim _{{\mathcal N}\gt 0,\ {\mathcal B}\rightarrow +\infty } \frac {\sigma _{\textit{ms}}}{\varepsilon }=0, \quad \lim _{{\mathcal B}\gt \frac {1}{2},\ {\mathcal N}\rightarrow +\infty } \frac {\sigma _{\textit{ms}}}{\varepsilon }=0. \end{equation}

Figure 4(c) shows $\sigma _{\textit{ms}}/\varepsilon$ versus $\mathcal B$ for ${\mathcal N}=0,\ 0.5,\ 1$ and ${\mathcal N}=2,$ so that $({\mathcal B},{\mathcal N})\in [0,+\infty )\times [\sqrt {5}/2,+\infty ),$ while figure 4(b) illustrates the variation of $\sigma _{\textit{ms}}/\varepsilon$ versus $\mathcal B$ for ${\mathcal N}=2,$ so that ${({\mathcal B},{\mathcal N})\in {\mathcal D}_f}.$ These figures show the expected agreement between the numerical results and the asymptotic formulae.

Accordingly, we may conclude that, unlike the non-stratified magnetized case for which $\lim _{{\mathcal B}\rightarrow +\infty }\sigma _{\textit{ms}}/\varepsilon ={1}/{2},$ the destabilizing resonance between two slow modes is weakened in the presence of stable stratification in the sense that its growth rate is reduced and approaches zero for strong stratification. Given (4.25), this clearly proves that the ${\mathcal N}\rightarrow 0$ limit is, in fact, singular (discontinuous).

4.3. Destabilizing resonance between a slow mode and a fast mode

In this section, we will show that the resonant case associated with $\omega _1-\omega _4=1$ is destabilizing. In contrast, the resonant case associated with $\omega _1-\omega _3=1$ is stable.

By substituting the expression of $\omega _1^2$ and $\omega _3^2$ given by (3.31) into the following identity:

(4.27) \begin{equation} \omega _1-\omega _4=\omega _1+\omega _3=1\implies \omega _1^4+\omega _3^4-2\omega _1^2\omega _3^2-2\big(\omega _1^2+\omega _3^2\big)+1=0, \end{equation}

we obtain

(4.28) \begin{equation} \big (\big (4-{\mathcal N}^2\big )^2+16{\mathcal B}^2\big )\mu ^4-2\!\big (2{\mathcal B}^2+\big ({\mathcal N}^2-4\big )\big ({\mathcal N}^2-1\big )\big )\mu ^2+\big ({\mathcal N}^2-1\big )^2=0. \end{equation}

We note that the resonant condition $\omega _1-\omega _3=1,$ is also described by (4.16), and the real roots of (4.16), such that $0\lt \mu ^2\lt 1,$ if they exist, are ordered. The smaller of the two roots characterizes the resonant cases characterized by $\omega _1-\omega _4=1$ where

(4.29) \begin{equation} \left ({\mathcal B},{\mathcal N}\right )\in (0,+\infty )\times \left [0,1\right ]\cup {\mathcal D}_m,\quad \lim _{\substack {{\mathcal B} \to +\infty {\mathcal N}\in [0,+\infty [}} {\mathcal B}^2\mu ^2=\frac {\left (1-{\mathcal N}^2\right )^2}{4}. \end{equation}

Here, the domain ${\mathcal D}_m$ is defined as follows:

(4.30) \begin{equation} \forall \left ({\mathcal B},{\mathcal N}\right )\in {\mathcal D}_m \Leftrightarrow \ \text{for given}\ {\mathcal N}\in \left ]1,+\infty \right [\ \implies \ 0\lt \sqrt {{3{\mathcal N}^2\!\left ({\mathcal N}^2-1\right )}}\leqslant {\mathcal B}\ . \end{equation}

Accordingly, for given ${\mathcal N}\in [0,1[,$ $\mu ({\mathcal B},{\mathcal N})$ decreases from $\mu (0,{\mathcal N})=\sqrt {(1-{\mathcal N}^2)/(4-{\mathcal N}^2)}$ as $\mathcal B$ increases approaching zero when ${\mathcal B}\rightarrow +\infty .$ For ${\mathcal N}=1,$ one has $\mu =0,$ and hence there is no instability for this resonant case (see the end of § 3.3.2). For given ${\mathcal N}\gt 1,$ such that $ ({\mathcal B},{\mathcal N} )\in {\mathcal D}_m,$ $\mu ({\mathcal B},{\mathcal N} )$ decreases from $\mu ({\mathcal B}_m,{\mathcal N} )=\sqrt { ({\mathcal N}^2-1 )/ (7{\mathcal N}^2-4 )}$ , approaching zero when ${\mathcal B}\rightarrow +\infty ,$ where ${\mathcal B}_m^2={{3{\mathcal N}^2 ({\mathcal N}^2-1 )}}$ (see also figure 2 c).

The maximal growth rate of the destabilizing resonant case $\omega _1-\omega _4=1,$ is denoted by $\sigma _{\textit{mm}}$ ( $m$ for mixed) is then described by (4.14) which, in this case, reduces to

(4.31) \begin{equation} \sigma _{\textit{mm}}=\frac {\varepsilon }{2\pi }\left (\Re \alpha \right )_{\textit{max}}=\frac {\varepsilon }{2\pi }\Re \sqrt {\tilde J_{14}\tilde J_{41}} \end{equation}

where

(4.32) \begin{align} \frac {{\tilde J}_{14}}{2\textrm{i} \pi }&= \frac {\mu \sqrt {1-\mu ^2}}{\left (\omega _3^2-\omega _1^2\right )}\left [ \frac {\left (\mu ^2{\mathcal B}^2-\omega _3^2\right )\!\omega _1\omega _3}{\mu ^2{\mathcal B}^2}+\frac {\left (\mu ^2{\mathcal B}^2-\omega _3^2\right )^2}{4\mu ^4{\mathcal B}^2}\frac {\omega _1}{\omega _3}\big (\omega _3-\big (1-\mu ^2\big ){\mathcal N}^2\big )\right . \nonumber \\&\quad \left . +\big (4\mu ^2-1\big )\omega _3- {\mathcal N}^2\big (1-\mu ^2\big)\frac {\left (\mu ^2{\mathcal B}^2-\omega _3^2\right )}{\omega _3} \right ]\!, \end{align}

(4.33) \begin{align} \frac {{\tilde J}_{41}}{2\textrm{i} \pi }&= \frac {\mu \sqrt {1-\mu ^2}}{\left (\omega _3^2-\omega _1^2\right )}\left [ \frac {\left (\mu ^2{\mathcal B}^2-\omega _1^2\right )\!\omega _1\omega _3}{\mu ^2{\mathcal B}^2}+\frac {\left (\mu ^2{\mathcal B}^2-\omega _1^2\right )^2}{4\mu ^4{\mathcal B}^2}\frac {\omega _3}{\omega _1}\big (\omega _1-\big (1-\mu ^2\big ){\mathcal N}^2\big )\right . \nonumber \\&\quad \left . +\big (4\mu ^2-1\big )\omega _1- {\mathcal N}^2\big (1-\mu ^2\big )\frac {\left (\mu ^2{\mathcal B}^2-\omega _1^2\right )}{\omega _1} \right ]\!, \end{align}

where detail calculations are given in the Supplementary material.

In the non-magnetized case ( ${\mathcal N}=0)$ the resulting expression of $\sigma _{\textit{mm}}$ is given by

(4.34) \begin{equation} \frac {\sigma _{\textit{mm}}}{\varepsilon }=\frac {{\mathcal B}^2\sqrt {3+4{\mathcal B}^2}}{8\left (1+{\mathcal B}^2\right )^2}\underset {{\mathcal B}\rightarrow +\infty }{\longrightarrow } 0. \end{equation}

However, in the magnetized case, analytically obtaining a compact form for the product ${\tilde J}_{14}{\tilde J}_{41}$ that only involves the parameters $\mu ,$ $\mathcal N$ and $\mathcal B$ is quite complex. Nevertheless, some results characterizing the behaviour of $\sigma _{\textit{mm}}$ can be drawn by analysing the behaviour of ${\tilde J}_{14}{\tilde J}_{41}$ for large ${\mathcal B}.$ Recall that, at small scales, i.e. $\sqrt {k_{\!p}^2+k_0^2}\gg 1,$ $\mathcal B$ is large even if $B_0/\varOmega _0\sim 1.$

Using the second identity in (4.29) we deduce the following equivalent forms for the normalized frequencies $\omega _1$ and $\omega _3$ given by (3.31),

(4.35) \begin{equation} \omega _1\underset {{\mathcal B}\rightarrow +\infty }{\large {\sim }}\frac {1+{\mathcal N}^2}{2},\quad \omega _3\underset {{\mathcal B}\rightarrow +\infty }{\sim }\mu {\mathcal B}\underset {{\mathcal B}\rightarrow +\infty }{\sim }\frac {1-{\mathcal N}^2}{2}, \end{equation}

so that $\omega _1+\omega _3=\omega _1-\omega _4=1.$ The resulting equivalent form for the product ${\tilde J}_{14}{\tilde J}_{41}$ is given by

(4.36) \begin{equation} \frac {{\tilde J}_{14}{\tilde J}_{41}} {4\pi ^2} \underset {{\mathcal B}\rightarrow +\infty }{\sim } \frac {\omega _3^2}{4\mu ^2{\mathcal B}^2} \frac {\left (\mu ^2{\mathcal B}^2-\omega _1^2\right )^2}{\left (\omega _3^2-\omega _1^2\right )^2}\left (1-\frac {{\mathcal N}^2}{\omega _1}\right ) \underset {{\mathcal B}\rightarrow +\infty }{\sim }\frac {1}{4}\frac {\left (1-{\mathcal N}^2\right )}{\left (1+{\mathcal N}^2\right )}. \end{equation}

Thus, the maximal growth rate of the destabilizing resonance between a fast mode and a slow mode such that $\omega _1-\omega _4=1$ has the following behaviour:

(4.37) \begin{equation} \frac {\sigma _{\textit{mm}}}{\varepsilon }\underset {{\mathcal B}\rightarrow +\infty }{\sim } \frac {1}{2}\sqrt {\frac {1-{\mathcal N}^2} {1+{\mathcal N}^2}}. \end{equation}

Figure 4(d) shows the variation of $\sigma _{\textit{mm}}/\varepsilon$ versus $\mathcal N$ for ${\mathcal B}=10^2,\ 10^3,\ 10^4$ and ${\mathcal B}\rightarrow +\infty .$ The numerical results obtained for $\varepsilon =0.05$ and ${\mathcal B}=10^2$ are also displayed in figure 4(d). As it can be seen, there is a good agreement between the numerical results (see § 5) and the asymptotic formulae.

Because in the non-magnetized case ( ${\mathcal N}=0),$ $\sigma _{\textit{mm}}$ tends to zero as ${\mathcal B}\rightarrow +\infty$ (see 4.34), while in the magnetized case (4.37) gives

(4.38) \begin{equation} \lim _{\substack {{\mathcal B} \to +\infty {\mathcal N}\rightarrow 0^+}}\frac {\sigma _{\textit{mm}}}{\varepsilon }=\frac {1}{2}, \end{equation}

we can conclude that for the destabilizing resonance associated with $\omega _1-\omega _4=1,$ the limit ${\mathcal N}\rightarrow 0$ is also, in fact, singular (discontinuous). This signifies that, at large ${\mathcal B},$ the presence of stable stratification enhances the destabilizing resonance between a fast (hydrodynamic) mode and a slow (magnetic) mode, provided $0\lt {\mathcal N}\lt 1.$ This would be relevant for the study of the generation of mean electromotive force ( ${\mathcal E}=\langle {\boldsymbol{u}}\times {\boldsymbol{b}}\rangle ,$ (see, Moffatt Reference Moffatt1978; Moffatt   Dormy Reference Moffatt and Dormy2019) in strongly magnetized precessing flows under a weak (stable) axial stratification, especially since in the present analysis the electromotive force generation mechanism does not require diffusivity to operate (see, Mizerski, Bajer   Moffatt Reference Mizerski, Bajer and Moffatt2012) for the case of the magneto-elliptic instability of rotating systems).

5.1. Identification of the instability tongues

For fixed $({\mathcal B},{\mathcal N}),$ we use the resonance conditions (4.16), and (4.28) for the identification of instabilities of order $n=1$ (if they exist). Note that the (3.32) also allows us to identify the high-order destabilizing resonances ( $n\geqslant 2)$ which are excluded by the procedure leading to the asymptotic formulae. On the other hand, for the same value of the couple $({\mathcal B} , {\mathcal N} ),$ we consider several values of $\varepsilon$ uniformly distributed in the interval $ (0, 0.25 )$ and, for each of these values of $\varepsilon ,$ we integrate numerically the Floquet system (4.1) and determine the growth rate $\Re \sigma$ (corresponding to the maximum modulus eigenvalue) for $5000$ values of $\mu$ (evenly distributed) in $ (0, 1 ).$ Obviously, $\sigma (\mu )=0$ if there is no instability. Thus, in the plane $(\mu , \varepsilon +\Re \sigma ) ,$ the region of instability that emanates from the point of abscissa $\mu$ characterizes a subharmonic instability as illustrated by figure 5 obtained for ${\mathcal B}=0.75$ and ${\mathcal N}=0.25.$ For these values, the destabilizing resonances detected by the computations are that between two fast modes which correspond to the tongue emanating from the point ( $\mu = 0.196, \varepsilon +\Re \sigma =0),$ and that between a fast mode and a slow mode which corresponds to the tongue emanating from the point ( $\mu = 0.380, \varepsilon +\Re \sigma =0$ ). The third instability tongue detected by the calculations emanates from the point ( $\mu = 0.434,\varepsilon +\Re \sigma =0$ ) and characterizes the destabilizing resonance of order $n=2$ between two fast modes. We note that this instability is more developed in the MBF case than in the KBF case, recalling that the matrix $\boldsymbol{D}$ in (3.23) does not have the same expression for these two base flows (see (3.24) and (A13)), except at order ${ O}(\varepsilon )$ (see (B4)).

Figure 5. Magneto-gravity-precessional instabilities. The figure shows, for fixed $\varepsilon ,$ $\varepsilon +\Re \sigma$ versus $\mu$ for ${\mathcal B}=0.75$ and ${\mathcal N}=0.25$ and $100$ values of $\varepsilon$ evenly distributed in the interval $[0,0.25].$ (a) The case of KBF; (b) the case of MBF. The instability region emanating from the point ( $\mu = 0.196, \varepsilon +\Re \sigma =0)$ is associated with the destabilizing resonance (of order $n=1)$ between two fast modes, while that emanating from the point ( $\mu = 0.380, \varepsilon +\Re \sigma =0$ ) is associated with the destabilizing resonance (of order $n=1)$ between a fast mode and a slow mode. The third region, which emanates from the point ( $\mu = 0.434,\varepsilon +\Re \sigma =0$ ), characterizes the destabilizing resonance of order $n=2$ between two fast modes.

For both background flows (KBF and MBF), the agreement between the asymptotic formulae and the numerical results is quite good for sufficiently small $\varepsilon (\leqslant 0.05),$ as illustrated in figures 4, 6 and 7. On the other hand, when $\varepsilon$ is not small enough, say $0.05\lt \varepsilon \lt 0.25$ , the subharmonic instabilities are strongest only in the KBF case, so the agreement between the asymptotic formulae and the numerical results is also quite good (see figures 7 a and 8 a).

Figure 6. Magneto-gravity-precessional instabilities. Maximal growth rate of dominant instability normalized by $\sigma _0=(5\sqrt {15}/32)\varepsilon$ plotted as a function of $0\leqslant {\mathcal B}=K_0B_0/\varOmega _0\leqslant 8$ and $0\leqslant {\mathcal N}=N/\varOmega _0\leqslant 2$ : (a) numerical results for $\varepsilon =0.05$ ; (b) asymptotic analysis results.

Figure 7. Magneto-gravity-precessional instabilities. Maximal growth rate of dominant instability normalized by $\sigma _0=(5\sqrt {15}/32)\varepsilon$ plotted as a function of $0\leqslant {\mathcal B}=K_0B_0/\varOmega _0\leqslant 8$ and $0\leqslant {\mathcal N}=N/\varOmega _0\leqslant 2$ for $\varepsilon =0.2.$ Here (a) KBF; (b) MBF.

Figure 8. Magneto-gravity-precessional instabilities. Maximal growth rate of dominant instability normalized by $\sigma _0=(5\sqrt {15}/32)\varepsilon$ plotted as a function of ${\mathcal N}=N/\varOmega _0$ for ${\mathcal B}=k_0B_0/\varOmega _0=5$ and $\varepsilon =0.05,\ 0.10,\ 0.20.$ Here (a) KBF; (b) MBF.

In figures 6 and 7, we show the continuous variation of the maximal growth rate $\sigma _m$ (maximum $\sigma$ over $0\leqslant \mu \leqslant 1)$ of the dominant instability, normalized by $\sigma _0=(5\sqrt {15}/32)\varepsilon ,$ plotted as a function of $0\leqslant {\mathcal B}\leqslant 8$ and $0\leqslant {\mathcal N}\leqslant 2.$ Figure 6(a) displays the numerical results obtained for $\varepsilon =0.05$ in the KBF case, while figure 6(b) shows the analytical results $\sigma _m/\sigma _0$ where

(5.1) \begin{equation} \sigma _m=\max (\sigma _{\textit{mf}},\sigma _{\textit{ms}},\sigma _{\textit{mm}}), \end{equation}

recalling that $\sigma _{\textit{mf}},$ $\sigma _{\textit{ms}}$ and $\sigma _{\textit{mm}}$ are given by (4.22) and (4.31), respectively. The numerical results for $\varepsilon =0.2$ are shown in figure 7 (figure 7 a, KBF case; figure 7 b, MBF case). Figure 8 shows $\sigma _m/\varepsilon$ as a function of $\mathcal N$ for ${\mathcal B}=5,$ $\varepsilon =0.05,\ 0.1$ and $\varepsilon =0.2.$ The analytical results are also reported in figure 8.

In the case of MBF and when $\varepsilon$ is not sufficiently small, it is rather the destabilizing resonance of order $n=2$ between two slow modes which is the strongest provided that ${\mathcal N}\gt 1$ and ${\mathcal B}\gt 2,$ as illustrated by figures 7(b) and 8(b). Recall that high-order destabilizing resonances ( $n\geqslant 2)$ are excluded by the procedure leading to the asymptotic formulae.

5.2. Diffusive effects in the case where $\nu =\eta =\kappa$

In this section, we consider the case where $\nu =\eta =\kappa ,$ so that the magnetic and thermal numbers are equal to one. In this case and given the definition (3.29), we show that the frequency of MIG waves is given by

(5.2) \begin{equation} \omega _j^{(v)}=\omega _j+\textrm{i} \nu k^2\quad (j=1,2,3,4), \end{equation}

where the superscript $(v)$ characterizes the viscous case and $\omega _j$ is the frequency of MIG waves for a non-diffusive fluid which is given by (3.30). Therefore, the condition of resonance

(5.3) \begin{equation} \omega _j^{(v)}-\omega _\ell ^{(v)}=\omega _j-\omega _\ell =n\varOmega _0,\quad (\ell \ne j) \end{equation}

is not modified when $\nu =\eta =\kappa .$ As for the destabilizing resonances (if they exist), they are affected by viscosity as shown in the following.

We formally rewrite the system (3.20) as

(5.4) \begin{equation} d_\tau \hat c_j={\mathcal L}_{\textit{jm}}\hat c_m \quad (j,m=1,{\cdots} ,5). \end{equation}

By taking into account diffusive terms such that $\nu =\eta =\kappa ,$ the above system is transformed as

(5.5) \begin{equation} d_\tau \hat c_j^{(v)}={\mathcal L}_{\textit{jm}}\hat c_m^{(v)} -\varOmega _0^{-1}\nu k^2\hat c_j^{(v)}\quad (j,m=1,{\cdots} ,5), \end{equation}

where its solution can be expressed in terms of that of (5.4) as

(5.6) \begin{equation} \hat c_j^{(v)}(\tau )=\hat c_j(\tau )\exp \!\left (-\frac {\nu }{\varOmega _0}\int _0^\tau k^2(s){\rm d}s\right )\! . \end{equation}

Thus, the linear part of the spectral counterpart of MIPS is zero because, for purposes of studying stability, we have set $\hat \varpi _m=0$ (see (3.22)),

(5.7) \begin{equation} k_0\varOmega _0^{-1}\hat \varpi ^{(v)}_m(\tau )=-\big ({\mathcal N}\hat c_4^{(v)}+\textrm{i} {\mathcal B}\hat \vartheta ^{(v)}\big )=k_0\varOmega _0^{-1}\hat \varpi _m \exp \!\left (-\frac {\nu }{\varOmega _0}\int _0^\tau k^2(s){\rm d}s\right )=0. \end{equation}

Accordingly, for a diffusive fluid such that $\nu =\eta =\kappa ,$ the growth rate of the destabilizing resonance (if it exists) takes the form (see (4.1) and (4.4))

(5.8) \begin{equation} \sigma ^{(v)}=\sigma -\frac {\nu }{2\pi \varOmega _0}\int _0^{2\pi } k^2(s){\rm d}s=\sigma -Re^{-1}L_0^2\!\left (k_0^2+k_{\!p}^2\right )\left (1+2\varepsilon ^2\!\left (1-\mu ^2\right )\right )\! , \end{equation}

where $Re=\varOmega _0L_0^2/\nu$ is the Reynolds number and $L_0$ a characteristic length scale. For the case of MBF (see (2.19) and Appendix A), we obtain

(5.9) \begin{equation} \sigma ^{(v)}=\sigma -Re^{-1}L_0^2\!\left (k_0^2+k_{\!p}^2\right )\left (1+4\varepsilon ^2\mu ^2\right )\! . \end{equation}

Thus, at order ${ O}(\varepsilon ),$ the effect of viscosity is the same for both basic flows. At ${L_0\sqrt {k_0^2+k_{\!p}^2}\sim 1},$ the dominant subharmonic instability survives the diffusive effects if $Re \gt Re_c=\sigma _m^{-1}.$