# Local magneto-shear instability in Newtonian gravity

Thomas Celora *  |,1 Ian Hawke *  |,1 Nils Andersson *  |1 and Gregory L.
Comer *  |2
1 School of Mathematics and STAG Research Centre, University of Southampton,
Southampton SO17 1BJ, UK
2 Department of Physics, Saint Louis University, St Louis, MO 63156-0907, USA
E-mail<EMAIL_ADDRESS><EMAIL_ADDRESS>
###### Abstract

The magneto-rotational instability (MRI)—which is due to an interplay between
a sheared background and the magnetic field—is commonly considered a key
ingredient for developing and sustaining turbulence in the outer envelope of
binary neutron star merger remnants. To assess whether (or not) the
instability is active and resolved, criteria originally derived in the
accretion disk literature—thus exploiting the symmetries of such systems—are
often used. In this paper we discuss the magneto-shear instability as a truly
local phenomenon, relaxing common symmetry assumptions on the background on
top of which the instability grows. This makes the discussion well-suited for
highly dynamical environments such as binary mergers. We find that—although
this is somewhat hidden in the usual derivation of the MRI dispersion
relation—the instability crucially depends on the assumed symmetries. Relaxing
the symmetry assumptions on the background we find that the role of the
magnetic field is significantly diminished, as it affects the modes’ growth
but does not drive it. We conclude by making contact with a suitable filtering
operation, as this is key to separating background and fluctuations in highly
dynamical systems.

###### keywords:

instabilities – stars: neutron – stars: oscillations – stars: rotation

††pubyear: 2023

##  1 Introduction

The magneto-rotational instability was discovered by Balbus and Hawley in the
early 1990s (Balbus & Hawley, 1991; Hawley & Balbus, 1991; Balbus & Hawley,
1998) (linking to earlier ideas from, for example, Chandrasekhar (1960) and
Velikhov (1959)). Due to the fast instability growth rate, this mechanism is
considered the most promising candidate for developing/sustaining magneto-
hydrodynamic turbulence in accretion disks as well as explaining enhanced
angular momentum transfer (Balbus & Hawley, 1998; Shakura & Sunyaev, 1973).
The instability is due to an interplay between a weak magnetic field and a
sheared background flow. With few exceptions (see for example Mahajan &
Krishan (2008) and Shakura et al. (2022)) and due to its “local” nature, the
magneto-rotational instability is commonly discussed in the so-called
“shearing box approximation” (Goldreich & Lynden-Bell, 1965; Hill, 1878). That
is, the instability is established in a frame that corotates with a fiducial
point in the mid-plane of the undisturbed disk (see also Goodman & Xu, 1994).
This is convenient for analytical studies as well as numerical analysis since
local simulations can reach much higher resolutions than global ones (see, for
example, Hawley et al., 1995; Zier & Springel, 2022, and references therein).
Shearing box simulations confirmed the predictions of the linear theory, and
allowed studies of non-linear features, such as the formation of the so-called
channel modes, eventually leading to a turbulent regime as these are destroyed
by a parasitic instability (Goodman & Xu, 1994). As turbulence decays, the
magneto-rotational instability may be revived and the process can start again
in cycles. Most importantly, while it has recently become possible to perform
fully kinetic simulations (see, e.g. Hoshino, 2015; Inchingolo et al., 2018),
shearing box simulations are necessary to explore the parameter space in a
variety of different setups (Sharma et al., 2006; Kempski et al., 2019; Guilet
et al., 2022; Held & Mamatsashvili, 2022).

Although originally discussed in the context of accretion disks, the magneto-
rotational instability is also thought to play a role in neutron-star mergers
(Duez et al., 2006; Siegel et al., 2013; Palenzuela et al., 2022; Margalit et
al., 2022; Hayashi et al., 2022; Kiuchi et al., 2022), especially for
sustaining a magneto-turbulent state in the outer envelope of the remnant,
where the Kelvin-Helmholtz instability is less significant or, indeed, not
active (Kiuchi et al., 2018). To assess whether or not the magneto-rotational
instability is active and resolved in merger simulations, criteria
discussed/established in the context of accretion disks (Hawley et al., 2011,
2013; Shibata, 2015) are often used. However, because binary neutron star
mergers are highly dynamical environments, framing a discussion of the
magneto-rotational instability using criteria that exploits restrictive
symmetry conditions might be misleading. Motivated by this, we aim to explore
the impact of relaxing common assumptions—well-motivated in the accretion
disks scenario, like an axisymmetric and circular background flow, but less so
for mergers—on the magneto-shear instability.

The paper is laid out as follows: we start in section 2 by introducing the
WKB-type expansion together with the concept of fast vs. slow background
gradients, as needed to derive dispersion relations associated with
perturbations propagating on top of a non-homogeneous background flow. In
section 3 we discuss how we can, in highly dynamical environments such as
mergers, still refer to an unspecified background—with no symmetries stated
from the outset—and consider perturbations rapidly evolving on top of that. We
continue in section 4 by deriving the non-inertial induction equation, and
making contact with the concept of a local frame. We do so as we will study
the instability from the perspective of an observer moving around with the
background flow, which cannot be expected to be inertial on a general basis.
We discuss our main findings in sections 5 and 6. These are based on the
dispersion relations, which are derived in detail in appendices C and D in
order to be able to focus on the physics in the main body of the paper. In
particular, we discuss the magneto-shear instability paying careful attention
to the differences introduced by the addition of a magnetic field. We draw our
conclusions, and make contact with the magneto-rotational instability (and the
associated criteria) in section 7.

##  2 Background gradients and plane-wave expansion

Let us begin by noting that the magneto-rotational instability is, in some
sense, a “global instability analyzed with local tools”. The local nature is
evident since the instability is established by means of a dispersion relation
(hence involves a plane-wave expansion and, by assumption, a short-wavelength
approximation). At the same time, one may appreciate the “global nature” of
the instability by recalling the key aspects of the instability: the addition
of a weak magnetic field turns axisymmetric modes (which would otherwise be
hydrodynamically stable) unstable. The global axisymmetry of the background,
then, plays a crucial role as the relevant hydrodynamic stability
criterion—the Rayleigh criterion (Rayleigh, 1917)—applies to axisymmetric
modes only. Although the standard derivation of the instability does not
highlight this subtlety, this aspect becomes apparent if we formulate the
problem using a co-rotating local frame (cf. the discussion in appendices A
and B).

With these points in mind, let us spell out how we intend to discuss the
magneto-shear instability without referring to a given axisymmetric and
circular background. Consistent with the shearing box idea (Goldreich &
Lynden-Bell, 1965; Hill, 1878), the strategy is to zoom in on a small region
of fluid—small enough for the analysis to be local but large enough to allow
for a meaningful hydrodynamic description. We then set up a local Cartesian
frame co-moving with the background flow—so that the background velocity
vanishes at the origin of the local box. As this frame moves with the flow—and
hence cannot be expected to be inertial—we need to consider the (at this point
Newtonian) ideal magneto-hydrodynamics equations in a non-inertial frame. This
step is commonly left out of the discussion, so we first of all have to fill
this gap. The non-inertial equations will then be perturbed—retaining
gradients in the background quantities as explained below—and a local WKB-type
dispersion relation will be derived and studied. This way we can account for
the effects of a sheared background and its interplay with the magnetic field
in a general setting.

Strictly speaking, the plane-wave expansion only makes sense for a homogeneous
background—that is, the plane-wave amplitude is assumed to vary on the same
scales as the background. At the same time, we know that a sheared background
is key to the magneto-rotational instability. Therefore, given any
quantity/field $a$, we first write it as a sum of background plus
perturbations

$a=A+\delta A\;,$ (1)

and then introduce a WKB-type expansion of the form (Thorne & Blandford, 2017;
Anile, 1989)

$\delta
A=\left(\sum_{q=0}\epsilon^{q}\bar{A}_{q}\right)e^{i\theta/{\epsilon}}\;\,\bar{\delta},$
(2)

with book-keeping parameters $\bar{\delta}$ and ${\epsilon}$ (see also
Palapanidis, 2018). The former $\left(\bar{\delta}\right)$ is introduced to
measure the relative magnitude of background vs. perturbations, while the
latter ($\epsilon$) is given by ${\epsilon}\approx\lambda/L$ where $\lambda$
is the typical wavelength of the waves and $L$ is the typical lengthscale over
which the wave amplitude, polarization and wavelength vary. Having split the
perturbations into amplitude and phase, we follow the standard convention
(Misner et al., 2017) and stick all “post-geometric optics” corrections into
the amplitude $\bar{A}_{q}$. With this Ansatz, the background equations are
obtained by collecting all terms of order
${\mathcal{O}}(\bar{\delta}^{0},{\epsilon}^{0})$, while the perturbation
equations are obtained collecting terms of order
${\mathcal{O}}(\bar{\delta}^{1},{\epsilon}^{0})$. Terms of higher order in
${\epsilon}$ correspond to post-geometric optics, while those of higher order
in $\bar{\delta}$ represent non-linear perturbations.

Along with this WKB-type Ansatz, we need to introduce the concept of fast and
slowly varying quantities. Given a specific choice of coordinates, a quantity
is slow in the variable $x$ if $A=A(X)$ where $X={\epsilon}x$ while it is fast
if $A=A(x)$. Deciding which quantities are fast or slow corresponds to
specifying (in a qualitative manner) the background configuration. As an
illustration, consider the simple toy problem

$a(\partial_{x}b+\partial_{x}c)=0\;,$ (3)

together with the Ansatz from eq. 2. Let us first assume that both $B$ and $C$
are fast, so that
$\partial_{x}B\approx{\mathcal{O}}(\bar{\delta}^{0},{\epsilon}^{0})$ and
similarly for $C$. The background equation is then

$A\left(\partial_{x}B+\partial_{x}C\right)=0\;.$ (4)

If we instead assume that, say, $B$ is fast while $C$ is slow, then
$\partial_{x}B\approx{\mathcal{O}}(\bar{\delta}^{0},{\epsilon}^{0})$ while
$\partial_{x}C\approx{\mathcal{O}}(\bar{\delta}^{0},{\epsilon})$ and the
background equation becomes

$A\partial_{x}B=0\;.$ (5)

Clearly, the two problems are different already at the background level.

Let us now turn to the linear perturbations. Because we have explicitly
introduced the book-keeping parameter ${\epsilon}$ in eq. 2, we take all
amplitude terms as well as the phase to be slowly varying. Then, to order
${\mathcal{O}}(\bar{\delta},{\epsilon}^{0})$ we have

$\left(A+\bar{\delta}\bar{A}_{0}e^{i\theta/{\epsilon}}\right)\partial_{x}\left[B+\bar{\delta}\bar{B}_{0}e^{i\theta/{\epsilon}}+C+\bar{\delta}\bar{C}_{0}e^{i\theta/{\epsilon}}\right]=0\;.$
(6)

Assuming again that the background quantity $B$ is fast, while $C$ is slow,
the perturbation equation becomes

$\bar{A}_{0}\left(\partial_{x}B\right)e^{i\theta/{\epsilon}}+A\left(\bar{B}_{0}\partial_{x}e^{i\theta/{\epsilon}}+\bar{C}_{0}\partial_{x}e^{i\theta/{\epsilon}}\right)=0\;.$
(7)

Next, Taylor expanding the phase—which is slowly varying—we get

$\frac{\theta(x)}{{\epsilon}}\approx\frac{\theta(0)}{{\epsilon}}+\partialderivative{\theta}{X}\Big{|}_{X=0}x+\dots=\theta(0)/{\epsilon}+k_{x}x+{\mathcal{O}}({\epsilon})\;,$
(8)

where we define the wave-vector $k_{x}=\partial\theta/\partial X$ from the
first order term in the expansion, while the overall constant can be
neglected. Then, introducing an analogous expansion for the fast background
gradients $\partial_{x}B(x)=\partial_{x}B(0)+{\mathcal{O}}({\epsilon})$ we end
up with

$\bar{A}_{0}(\partial_{x}B)+A\left(ik_{x}\bar{B}_{0}+ik_{x}\bar{C}_{0}\right)=0\;,$
(9)

where both $\partial_{x}B$ and $A$ are evaluated at a point (conveniently
chosen as the origin of the coordinate system). Therefore, if all background
quantities are “slow”, we get back the dispersion relation we would have
obtained ignoring all background gradients. This is quite intuitive. However,
the strategy also allows us to account for the impact that “fast” background
gradients have on the dispersion relation. In short, as long as these terms
are treated as constants, we may retain them and work out a dispersion
relation in the usual way.

##  3 The slowly evolving background

The starting point for any hydrodynamic perturbation analysis is the
choice/identification of a stationary background flow configuration, which is
then perturbed in order to establish stability (or not). Here, we want to
frame the analysis of the magneto-shear instability without considering a
specific background configuration with constraining symmetries stated from the
outset. Nonetheless, we need to clarify how we can refer to a suitable
“background” in highly dynamical environments like binary neutron star
mergers. We will first analyse the problem analytically, and return to discuss
the link to numerical simulations in section 7. Given real numerical
simulation data, this discussion will inevitably involve some kind of
filtering operation. Anticipating that this can be done in a meaningful way,
we consider perturbations evolving rapidly with respect to the evolution time-
scale of an unspecified “background” flow.

To make this statement more precise, let us consider the inertial ideal MHD
equations and introduce reference values for each quantity (indicated with an
“$r$” subscript) such as $\rho=\rho_{r}\tilde{\rho}$. We introduce the
(dimensionless) Strouhal, Mach, Froude and magnetic interaction numbers as

${\varepsilon}_{\text{St}}=\frac{l_{r}}{t_{r}v_{r}}\;,\quad{\varepsilon}_{\text{Ma}}=\frac{v_{r}}{c_{r}}\;,\quad{\varepsilon}_{\text{Fr}}=\frac{v_{r}}{\sqrt{\Phi_{r}}}\;,\quad{\varepsilon}_{B}=\frac{B_{r}^{2}}{\mu_{0}\rho_{r}v_{r}^{2}}\;,$
(10)

where $l_{r},\,t_{r},\,v_{r}$ are characteristic lengthscale, timescale and
velocity (respectively) while $B_{r},\,\Phi_{r},\,\rho_{r}$ are reference
values for the magnetic field, gravitational potential and density and $c_{r}$
is the (adiabatic) speed of sound. This way, the non-dimensional inertial
ideal MHD equations read (now dropping the tildes for notational clarity)

${\varepsilon}_{\text{St}}\,\partial_{t}{\rho}=-\rho\nabla_{i}v^{i}-v^{i}\nabla_{i}\rho\;,$
(11a)
${\varepsilon}_{\text{St}}\,\partial_{t}{B^{i}}=-v^{j}\nabla_{j}B^{i}+B^{j}\nabla_{j}v^{i}-B^{i}\nabla_{j}v^{j}\;,$
(11b)
${\varepsilon}_{\text{St}}\,\partial_{t}{v^{i}}=-v^{j}\nabla_{j}v^{i}-\frac{1}{{\epsilon}_{\text{Ma}}^{2}}\frac{1}{\rho}\nabla^{i}\rho-\frac{1}{{\epsilon}_{\text{Fr}}^{2}}\nabla^{i}\Phi-\\\
{\epsilon}_{B}\frac{1}{\rho}\left[B^{j}\nabla_{j}B^{i}-\nabla^{i}\left(\frac{B^{2}}{2}\right)\right]\;.$
(11c)

From this we see that a generic flow configuration can be considered slowly
evolving (in time) as long as the corresponding Strouhal number is small. In
practice, given a characteristic lengthscale $l_{r}$ and velocity $v_{r}$ of a
generic flow, we consider disturbances evolving on timescales $t_{r}$ such
that ${\epsilon}_{\text{St}}\ll 1$—over which the background can be
effectively taken as stationary. In turn, this determines the time-scales over
which we expect the following results to be reliable.

###  3.1 Velocity gradient decomposition

In the following we will consider the impact that gradients in the background
flow velocity have on the time evolution of perturbations. It is then
convenient to introduce the standard decomposition of the velocity gradient
into expansion, shear and vorticity. That is,

$\nabla_{i}v_{j}=\frac{1}{3}\theta g_{ij}+\sigma_{ij}+{\omega}_{ij}\;,$ (12)

where111We follow the standard notation and use (straight) round brackets to
indicate (anti-)symmetrization of the indices enclosed.

$\displaystyle\theta=\nabla_{i}v^{i}\;,$ (13a)
$\displaystyle\sigma_{ij}=\nabla_{(i}v_{j)}-\frac{1}{3}\theta
g_{ij}=\frac{1}{2}\left(\nabla_{i}v_{j}+\nabla_{j}v_{i}\right)-\frac{1}{3}\theta
g_{ij}\;,$ (13b)
$\displaystyle\omega_{ij}=\nabla_{[i}v_{j]}=\frac{1}{2}\left(\nabla_{i}v_{j}-\nabla_{j}v_{i}\right)\;.$
(13c)

In order to bring out the magneto-shear nature of the instability, we will
consider the impact of having a background with non-negligible shear and
vorticity separately. We will, however, not consider the impact of a
background expansion rate as exact non-linear results are sufficient to
predict this. In fact, due to the Alfvén theorem, we know that the magnetic
intensity must grow in a (ideal magneto-)fluid undergoing compression as the
field lines are squeezed together. Similarly, the field will get weaker in an
expanding fluid. In essence, we expect—and have verified explicitly—this non-
linear prediction to emerge in the analysis as a generic “instability”. The
background magnetic field cannot grow in time as it is assumed to be slowly
evolving by construction, so the required growth must be represented by
perturbations.

Before we move on to derive the non-inertial equations, it is useful to take a
brief detour and consider a realization of a flow with only non-negligible
shear. Because we are interested in flows that are slowly evolving we can
start by assuming that $v^{i}=v^{i}({\epsilon}t,x,y,z)$ and suppress the time
dependence in the following. We then take the velocity vector as mainly two-
dimensional, specifically in the $x-y$ plane of a set of local Cartesian
coordinates

$\mathbf{v}=v^{x}(x,y,z)\hat{x}+v^{y}(x,y,z)\hat{y}+{\mathcal{O}}({\epsilon})\;,$
(14)

where, in order to make sure the expansion is small, we take
$\partial_{x}{v^{x}}=-\partial_{y}{v^{y}}+{\mathcal{O}}({\epsilon})$. The
shear matrix is then given by

$\boldsymbol{\sigma}=\begin{pmatrix}\partial_{x}{v^{x}}&\frac{1}{2}\left(\partial_{x}{v^{y}}+\partial_{y}{v^{x}}\right)&\frac{1}{2}\partial_{z}{v^{x}}\\\
\frac{1}{2}\left(\partial_{x}{v^{y}}+\partial_{y}{v^{x}}\right)&-\partial_{x}{v^{x}}&\frac{1}{2}\partial_{z}{v^{y}}\\\
\frac{1}{2}\partial_{z}{v^{x}}&\frac{1}{2}\partial_{z}{v^{y}}&0\end{pmatrix}+{\mathcal{O}}({\epsilon})\;,$
(15)

while the curl of $\mathbf{v}$ becomes

$\nabla\times\mathbf{v}=-\left(\partial_{z}{v^{y}}\right)\hat{x}+\left(\partial_{z}{v^{x}}\right)\hat{y}+\left(\partial_{x}{v^{y}}-\partial_{y}{v^{x}}\right)\hat{z}\;.$
(16)

This has to be of ${\mathcal{O}}({\epsilon})$ for background flows with only
non-negligibile shear, in which case

$\partial_{z}{v^{x}}=\partial_{z}{v^{y}}={\mathcal{O}}({\epsilon})\;,\qquad\partial_{y}{v^{x}}=\partial_{x}{v^{y}}+{\mathcal{O}}({\epsilon}),$
(17)

and as result the determinant of the shear matrix vanishes (more precisely, is
of order ${\mathcal{O}}({\epsilon})$). This is equivalent to saying that two
eigenvalues of the shear matrix are opposite and the third is zero (to order
${\mathcal{O}}({\epsilon})$). In essence, a mainly two-dimensional flow with
negligible expansion and vorticity is characterized by a shear matrix with one
zero eigenvalue, and hence a vanishing determinant. This will turn out to be a
useful observation when we discuss the dispersion relations in sections 5 and
6. We also note that having negligible expansion (although relevant for the
present analysis) is not strictly necessary to the argument.

##  4 Non-inertial equations and the local frame

The ideal magneto-hydrodynamics equations above hold in an inertial frame.
However, as an observer locally co-moving with the fluid cannot (in general)
be expected to be inertial, we need to consider the equations according to a
non-inertial observer. As such, let us now derive the Newtonian induction
equation in a general (non-inertial) frame starting from the Maxwell equations
formulated in covariant form. This is the natural starting point as
electromagnetism is inherently a relativistic theory. Moreover, we stress that
the electric and magnetic field are observer dependent quantities, and the
standard textbook form of the Maxwell equations assume an inertial
observer/frame from the outset.

We then take as our starting point the covariant Maxwell equations (with
indices $a,b,\ldots$ representing space-time components, in contrast to the
$i,j,k,\ldots$ components from before)

$\nabla_{a}F^{ba}=\mu_{0}j^{b}\;,\qquad\nabla_{[a}F_{bc]}=0\;,$ (18)

where $F^{ab}$ is the Faraday tensor while $j^{b}$ is the four-current. From
these relations we derive the non-inertial (relativistic) induction equation
and, finally, consider the Newtonian limit. Because we are interested in the
non-inertial equations associated with an observer locally co-moving with the
fluid, it is natural to decompose the covariant Maxwell equations in the
“fluid frame” (Andersson & Comer, 2021)—leading to a (fibration) formulation,
commonly used in cosmology (Barrow et al., 2007). Hence, we introduce a four-
velocity $U^{a}$ associated with a generic observer and decompose the Faraday
tensor and the charge current as

$F_{ab}=2U_{[a}e_{b]}+{\varepsilon}_{abc}b^{c}\;,\qquad j^{a}=\sigma
U^{a}+J^{a}\;,$ (19)

where

$e_{a}=F_{ab}U^{b}\;,\;\;b_{a}=\frac{1}{2}{\varepsilon}_{abc}F^{bc}\;,\text{
and }{\varepsilon}_{abc}={\varepsilon}_{dabc}U^{d}\;.$ (20)

With these definitions we can rewrite the Maxwell equations as (see e.g.
Ellis, 1973; Andersson, 2012; Andersson & Comer, 2021)

$\perp^{a}_{b}\nabla_{a}e^{b}-\mu_{0}\sigma=2W^{a}b_{a}\;,$ (21a)
$\perp^{a}_{b}\nabla_{a}b^{b}=-2W^{a}e_{a}\;,$ (21b)
$\perp_{ab}\dot{e}^{b}-{\varepsilon}_{abc}\nabla^{b}b^{c}+\mu_{0}J_{a}=\\\
e^{b}\left(\sigma_{ba}+\omega_{ba}-\frac{2}{3}\theta\perp_{ba}\right)+{\varepsilon}_{abc}a^{b}b^{c}\;,$
(21c) $\perp_{ab}\dot{b}^{b}+{\varepsilon}_{abc}\nabla^{b}e^{c}=\\\
b^{b}\left(\sigma_{ba}+\omega_{ba}-\frac{2}{3}\theta\perp_{ba}\right)-{\varepsilon}_{abc}a^{b}e^{c}\;,$
(21d)

where dots stand for co-moving time derivatives $U^{a}\nabla_{a}$ and
$\perp_{ab}=g_{ab}+U_{a}U_{b}$ is the projection orthogonal to the observer
four-velocity ($g_{ab}$ is the spacetime metric). In eq. 21, the terms on the
left-hand side should be familiar, while those on the right-hand side are
associated with gradients of the observer four-velocity222The shear, vorticity
and expansion are defined as the four dimensional version of eq. 13—that is
with additional projections to ensure they are flow orthogonal—while the
observer four-acceleration is $a^{b}=U^{a}\nabla_{a}U^{b}$ and
$W^{a}=\frac{1}{2}{\varepsilon}^{abc}{\omega}_{bc}$.. As such, they vanish
identically for an inertial observer, and hence do not appear in most textbook
discussions.

To derive the induction equation, we follow the usual logic (see Andersson et
al., 2021, for a recent relativistic discussion) and “massage” the Faraday
equation

$\underbrace{\perp_{ab}\dot{b}^{b}-b^{b}\left(\sigma_{ba}+\omega_{ba}-\frac{2}{3}\theta\perp_{ba}\right)}_{\sim
b/T}\\\
+\underbrace{{\varepsilon}_{abc}\left(\nabla^{b}e^{c}+a^{b}e^{c}\right)}_{\sim
e/L}=0\;,$ (22)

to see that $e\sim Lb/T\sim Vb$, with $L$ and $T$ typical length- and time-
scales and $V$ the associated velocity333In order to avoid confusion, let us
state clearly that the scales $L,\,T$ (and the associated velocity $V$) in
this argument are, in principle, not related to the ones introduced above
$l_{r},\,t_{r}$ and $v_{r}$. The aim here is to reproduce the usual argument
for dropping the displacement current—and derive the induction equation—not to
argue in which sense we can consider perturbations rapidly evolving on a non-
stationary (but slowly evolving) background. When we perturb the induction
equation later on, however, we are in practice assuming the displacement
current to be small compared to perturbations.. As long as the electric and
magnetic fields are slowly evolving, a similar dimensional analysis then leads
us to neglecting terms involving the electric field (i.e. the displacement
current) in the Ampère law, so that

$J^{a}=\frac{1}{\mu_{0}}\left({\varepsilon}_{abc}\nabla^{b}b^{c}+{\varepsilon}_{abc}a^{b}b^{c}\right)\;.$
(23)

As the Ampère law has now been demoted to a constraint on the charge current
(see Andersson et al., 2021), we need to introduce a closure relation for the
electric field. In the ideal case, this can be obtained by looking at the
electric field measured by an observer locally co-moving with the fluid. In
fact, because the local fluid four-velocity $u^{a}$ is linked to the generic
observer $U^{a}$ via

$u^{a}=W\left(U^{a}+v^{a}\right)\;,\quad U^{a}v_{a}=0\;,\quad
W=\left(1-v^{a}v_{a}\right)^{-1/2}\;,$ (24)

where $v^{b}$ is the spatial fluid velocity as measured by $U^{a}$, the
electric field measured by the fluid is (cf. eq. 20)

$F_{ab}u^{b}=W\left[e_{a}+{\varepsilon}_{abc}v^{b}b^{c}+U_{a}(v^{b}b_{b})\right]\;.$
(25)

In a perfect conductor, where charges easily flow, one would expect the
electric field to “short out” as the matter becomes locally charge neutral, so
that

$e_{a}+{\varepsilon}_{abc}v^{b}b^{c}=0\;.$ (26)

With this constraint, we can derive the induction equation from Faraday’s law.
To do so, note that

${\varepsilon}_{abc}{\varepsilon}^{cde}=U^{f}U_{g}{\varepsilon}_{fabc}{\varepsilon}^{gdec}=-3!U^{f}U_{g}\delta^{[g}_{f}\delta^{d}_{a}\delta^{e]}_{b}\\\
=\left(\delta^{d}_{a}\delta^{e}_{b}-\delta^{e}_{a}\delta^{d}_{b}\right)-\left(\mathrel{/\mkern-5.0mu/}^{d}_{a}\delta^{e}_{b}-\mathrel{/\mkern-5.0mu/}^{e}_{a}\delta^{d}_{b}\right)-\left(\delta^{d}_{a}\mathrel{/\mkern-5.0mu/}^{e}_{b}-\delta^{e}_{a}\mathrel{/\mkern-5.0mu/}^{d}_{b}\right)\;,$
(27)

where we introduced the parallel projection
$\mathrel{/\mkern-5.0mu/}^{a}_{b}=-U^{a}U_{b}$. When this is contracted with a
spatial tensor (with respect to $U^{a}$) the last two terms in eq. 27 can be
dropped. It follows that

${\varepsilon}_{abc}{\varepsilon}^{cde}\left(a^{b}v_{d}b_{e}\right)=\left(a^{b}b_{b}\right)v_{a}-\left(a^{b}v_{b}\right)b_{a}\;.$
(28)

We also need to take care of the curl of $e^{a}$ term in the Faraday equation.
This can be written

$-{\varepsilon}_{abc}\nabla^{b}\left({\varepsilon}^{cde}v_{d}b_{e}\right)=\\\
\left[-{\varepsilon}_{abc}\left(\nabla^{b}{\varepsilon}^{cde}\right)v_{d}b_{e}\right]-\left[{\varepsilon}_{abc}{\varepsilon}^{cde}\nabla^{b}\left(v_{d}b_{e}\right)\right]\;,$
(29)

where it is convenient to consider the two terms separately. We start from the
second term and, even if $U^{a}$ is not necessarily surface forming (i.e. has
non-vanishing vorticity), we introduce a “spatial” covariant derivative $D$ in
the usual way (projecting each index in the sub-space orthogonal to $U^{a}$).
Then, it is easy to see that

${\varepsilon}_{abc}{\varepsilon}^{cde}D^{b}\left(v_{d}b_{e}\right)={\varepsilon}_{abc}{\varepsilon}^{cde}\left(\perp^{b}_{f}\perp^{g}_{d}\perp^{h}_{e}\right)\nabla^{f}\left(v_{g}b_{h}\right)\\\
={\varepsilon}_{abc}{\varepsilon}^{cde}\nabla^{b}\left(v_{d}b_{e}\right)\;,$
(30)

and hence

$-{\varepsilon}_{abc}{\varepsilon}^{cde}\nabla^{b}\left(v_{d}b_{e}\right)=D^{b}\left(v_{b}b_{a}\right)-D^{b}\left(v_{a}b_{b}\right)\;.$
(31)

As for the other term, writing it as

$-{\varepsilon}_{abc}\left(\nabla^{b}{\varepsilon}^{cde}\right)v_{d}b_{e}=-U^{g}{\varepsilon}_{gabc}{\varepsilon}^{fcde}\left(\nabla^{b}U_{f}\right)v_{d}b_{e}\\\
=-U^{g}\delta^{[f}_{g}\delta^{d}_{a}\delta^{e]}_{b}\,g_{fh}\left(-U^{b}a^{h}+{\omega}^{bh}+\sigma^{bh}+\frac{1}{3}\theta\perp^{bh}\right)v_{d}b_{e}\;,$
(32)

we see that—given the anti-symmetrization—it vanishes identically.

In summary, the (relativistic) induction equation according to a generic
observer444The worldlines of the generic observer $U^{a}$ constitute a
fibration of the spacetime, hence we may call this the ideal induction
equation in the fibration framework, as opposed to the corresponding 3+1 form
derived by, for example, Andersson et al. (2021). can be written as

$\perp^{ab}\dot{b}_{b}+D_{b}(v^{b}b^{a})-D_{b}(v^{a}b^{b})=\\\
\left(\sigma^{ab}-\omega^{ab}-\frac{2}{3}\theta\perp^{ab}\right)b_{b}+v^{a}(a_{b}b^{b})-b^{a}(a_{b}v^{b})\;.$
(33)

The terms on the left should be familiar, while those on the right vanish for
an inertial observer.

Next, let us show how the derived equation further simplifies in the Newtonian
limit. On dimensional grounds, we observe that the last two terms on the
right-hand side of eq. 33 contain an extra factor of $1/c^{2}$ (with respect
to the others, where $c$ is the speed of light), and will as a result be
negligible in the non-relativistic limit ($c^{2}\to\infty$). Similarly, let us
also consider the absence of monopoles constraint. From eq. 21 and eq. 26 we
immediately obtain

$\perp^{a}_{b}\nabla_{a}b^{b}=2W^{a}{\varepsilon}_{abc}v^{b}b^{c}\;,$ (34)

and we observe the term on the left hand side is $\sim b/L$ while that on the
right is $\sim bL/T^{2}$. Dimensional consistency implies the term on the
right-hand side contains an extra factor of $1/c^{2}$ and should be neglected
in the Newtonian limit. In essence, non-inertial effects do not affect the
absence of monopoles constraint at the Newtonian level. When it comes to the
Lorentz force, we expect it not to change at the Newtonian level, but let us
nonetheless check this for consistency. The Lorentz four-force can be written

$-j_{b}F^{ba}=-\left(\sigma
U_{b}+J_{b}\right)\left(U^{b}e^{a}+U^{a}e^{b}+{\varepsilon}^{bacd}U_{c}b_{d}\right)\\\
=-U^{a}\left(J_{b}{\varepsilon}^{bcd}v_{c}b_{d}\right)+{\varepsilon}^{abc}\left(J_{b}-\sigma
v_{b}\right)b_{c}\;,$ (35)

where we used the ideal magneto-hydrodynamics relation (26) in the second
step. The Lorentz three force corresponds to the second term, where the charge
density is measured by the observer, hence does not (in general) vanish.
However, if we insist on the local charge density to be zero—consistently with
(26)—then we have (cf. eq. 19)

$-U^{a}j_{a}=W(\sigma-v_{a}J^{a})=0\Longrightarrow J_{b}-\sigma
v_{b}=\left(g^{a}_{b}-v^{a}v_{b}\right)J_{a}\;.$ (36)

Re-inserting the factor of $1/c^{2}$ we see that the second term is negligible
with respect to the first. As also the second term in eq. 23 is negligible in
the Newtonian limit, we see that the Lorentz force in the Euler equation is
unchanged (as expected).

###  4.1 The local frame of an observer

Having derived the relativistic induction equation according to a generic
(non-inertial) observer and recalling that we are interested in a local
analysis, we now make contact with the concept of local frame associated with
an observer (Gourgoulhon, 2013; Misner et al., 2017). Given an observer
worldline with tangent $U^{a}$, the local frame is constructed by considering
three spatial unit vectors that complete $U^{a}$ to an orthonormal basis on
the tangent space at any given point. These unit vectors (including the
observer four-velocity)—the components of which are indicated by hats—are then
transported along the worldline according to

$U^{\hat{c}}\nabla_{\hat{c}}e_{\hat{a}}={\Omega}^{\hat{b}}_{\;\hat{a}}e_{\hat{b}}\;,\quad{\Omega}_{\hat{a}\hat{b}}=U_{\hat{a}}a_{\hat{b}}-a_{\hat{a}}U_{\hat{b}}-{\epsilon}_{\hat{a}\hat{b}\hat{c}\hat{d}}U^{\hat{c}}W^{\hat{d}}\;,$
(37)

where $a^{\hat{a}}$ is the four-acceleration of $U^{\hat{a}}$ (an intrinsic
property of the worldline) and $W^{\hat{d}}$ is the _arbitrary_ four-rotation
of the local frame. Focusing on the first term in the relativistic induction
equation (33), and using eq. 37,

$\dot{b}_{\hat{b}}=\left(U^{\hat{c}}\nabla_{\hat{c}}b\right)_{\hat{b}}=U^{\hat{c}}\partial_{\hat{c}}b_{\hat{b}}+\left(a^{\hat{c}}b_{\hat{c}}\right)U_{\hat{b}}+{\varepsilon}_{\hat{b}\hat{e}\hat{c}}W^{\hat{e}}b^{\hat{c}}\;.$
(38)

The second term vanishes due to the orthogonal projection, while555Identifying
the vorticity of the fibration observer with the four-rotation of the chosen
local frame.

$\perp^{\hat{a}\hat{b}}\dot{b}_{\hat{b}}+{\omega}^{\hat{a}\hat{b}}b_{\hat{b}}=\perp^{\hat{a}\hat{b}}\left(U^{\hat{c}}\partial_{\hat{c}}b_{\hat{b}}\right)\;.$
(39)

In practice, the term involving the four-rotation of the frame drops out of
the induction equation. We also note that, because we are now considering the
non-inertial equations in the local frame of a single observer, there is no
shear or expansion (associated with the fibration of spacetime induced by the
observer). Given this, the induction equation in the Newtonian limit
simplifies to

$\partial_{t}B^{\hat{i}}+\nabla_{\hat{j}}\left(v^{\hat{j}}B^{\hat{i}}-v^{\hat{i}}B^{\hat{j}}\right)=0\;.$
(40)

At the Newtonian level then, the induction equation in the local frame of a
_generic_ observer retains the same form as for an inertial one. This is
similar to the case of the Lorentz force (entering the Euler equation) and the
Ampère law. However, this is only true in the Newtonian case. Additional terms
involving the four-acceleration of the observer will appear at the special
relativistic level so some care will be required in order to extend our
results in that direction.

When it comes to the non-inertial terms in the Euler equations, these are
obviously well known: we have to account for fictitious acceleration. We refer
to Gourgoulhon (2013) for a rigorous derivation in special relativity, showing
how additional terms involving the observer four-acceleration also enter the
relativistic expressions. Let us also stress that working with a rotating or
non-rotating local frame is entirely a matter of choice (Misner et al., 2017).
At the local Newtonian level, we can always get rid of the non-inertial terms
associated with the frame rotation and effectively work with the inertial
equations.

We conclude this section by noting that, as previously anticipated, some kind
of filtering operation is key to separate between background and fluctuations
in a highly dynamical environment. Postponing a discussion of this to section
7, let us simply note at this point that the notion of local frame discussed
here is closely linked to the covariant filtering procedure discussed in
Celora et al. (2021).

##  5 Going back to hydrodynamics

As briefly hinted at in section 2, the magneto-rotational instability relies
on the _hydrodynamic_ stability of axisymmetric modes. The generic instability
problem is more involved. If we relax the symmetry assumptions on the
background, we first of all need to consider the fact that hydrodynamic shear
flows tend to be unstable (Drazin & Reid, 2004). That is, we expect to find
instabilities appearing already at the hydrodynamic level. Clearly, such
instabilities would be affected by a magnetic field but not caused by it in
the first place. This is an important distiction seeing as the magneto-
rotational instability is _due to_ the presence of the magnetic field.

With this observation in mind, let us first consider the fluid problem. This
is important for two reasons: first, it will allow us to get a better grasp on
the magnetic field impact on the instability. Second, it will allow us to make
contact with the Rayleigh criterion (and ultimately the magneto-rotational
instability). As the fluid problem is much simpler than the magneto-fluid one,
we will study the case where both shear and vorticity gradients are retained,
and also discuss the impact of shear viscosity—either of microphysical origin
or due to filtering666Any filtering operation will introduce additional
_residual_ terms every time it acts upon a non-linear term (M. Lesieur, 2005;
Berselli et al., 2006; McDonough, 2014; Schmidt, 2015). These terms are akin
to (but not quite the same as) dissipative terms (Celora et al., 2021), and
are meant to capture transport to/from unresolved scales.. Shear viscosity is
introduced in the usual way (see Landau & Lifshitz, 1959), and the shear
viscosity coefficient $\eta$ will be considered constant, consistent with the
local analysis.

Because the calculation is a bit tedious and we want to focus on the
implications for the physics, the derivation of the relevant fluid dispersion
relation(s) is provided in appendix C. There we also express the coefficients
of the resulting characteristic polynomial(s) in terms of scalars built from
the background quantities. This allows us to keep the discussion as general as
possible, without having to refer to a specific background configuration. It
is, however, worth stressing that, in many situations of interest the relevant
dynamics is either sub- or supersonic. Given this, for these problems it is
worth considering models that filter out modes that are either faster or
slower than the sound waves. This can be done starting from a fully
compressible dispersion relation and taking either of two limits: either we
assume the speed of sound to be very large, in which case the model becomes
sound-proof (we point to Vasil et al., 2013, for more details), or very small.
In the following, we typically work in the sound-proof limit, noting that the
MRI is typically discussed within the so-called Boussinesq approximation
(Barletta, 2022), thus removing fast magneto-sonic waves777They do, however,
retain perturbations in the fluid pressure in the Euler equation as they
consider a non-barotropic equation of state and the impact of stratification.
(Balbus & Hawley, 1991).

Starting from the continuity equation, perturbing it and introducing the
plane-wave expansion we readily obtain

$\partial_{t}\delta\rho+\delta\rho\nabla_{i}v^{i}+v^{i}\nabla_{i}\delta\rho+\rho\nabla_{i}\delta
v^{i}=0\Longrightarrow-i\omega\delta\rho+i\rho k_{i}\delta v^{i}=0\;,$ (41)

where ${\omega}$ and $k_{i}$ are defined as in section 2. Note that we set
$v^{i}=0$ as we evaluate the relation at the centre of the local box, and
assume that the background expansion rate $\nabla_{i}v^{i}$ can be neglected.
Similarly, we write the perturbed Euler equation including a shear-viscous
term as

$\partial_{t}\delta v_{i}+\delta
v^{j}\nabla_{j}v_{i}+\frac{1}{\rho}\nabla_{i}\delta
P-\delta\left(\eta\nabla^{j}\tau_{ji}\right)=0\\\
\Longrightarrow-i\omega\delta
v_{i}+i\frac{c_{s}^{2}}{\rho}k_{i}\delta\rho+\sigma_{ij}\delta
v^{j}+{\epsilon}_{ijk}W^{j}\delta
v^{k}-\delta\left(\eta\nabla^{j}\tau_{ji}\right)=0\;,$ (42)

where $\tau_{ji}$ is the rate-of-strain tensor and
$W^{i}=1/2{\epsilon}^{ijk}{\omega}_{jk}$. Working this out we have retained
gradients in the background flow only, used the velocity gradient
decomposition (section 3.1), introduced the adiabatic speed of sound
$c_{s}^{2}=\partial P/\partial\rho$, and considered the gravitational
potential to be externally sourced (hence neglecting its perturbations).

As a first example, consider a background with negligible vorticity, set the
shear viscosity to zero and take the sound-proof limit. Then consider the case
$\text{det}(\boldsymbol{\sigma})=0$, and first of all look for modes such that
$\sigma_{ij}k^{j}=0$. Recalling that, as discussed in section 3.1, a mainly
two dimensional flow with negligible vorticity is characterized by having a
shear matrix with vanishing determinant—that is
$\text{det}(\boldsymbol{\sigma})\sim{\mathcal{O}}({\epsilon})$—and noting that
we can choose the orientation of the local axes in such a way that the
background flow is, say, along the $\hat{x},\,\hat{y}$ directions, we can
always consider the determinant to be zero. This means that there always
exists a wave-vector living in the eigen-space corresponding to the zero
eigenvalue888We note that these modes correspond to those propagating in the
$z$-direction in terms of the adapted Cartesian coordinates introduced in
section 3.1. We also stress that the same condition is satisfied by the
fastest growing unstable modes, which also propagate vertically, although with
respect to the “global” cylindrical system.. Taking these steps, we end up
with999Note that the same dispersion relation applies in the opposite limit
where $c_{s}$ is small. (cf. eq. 87)

${\omega}^{2}=-\frac{1}{2}\text{Tr}(\boldsymbol{\sigma}^{2})\Longrightarrow{\omega}=\pm
i\sqrt{\frac{1}{2}\text{Tr}(\boldsymbol{\sigma}^{2})}\;.$ (43)

These modes are non-propagating, and half of them are unstable with a growth
rate independent of the wave-vector. Next, we consider wave-numbers such
that101010Using the adapted coordinates of section 3.1 it is possible to see
that these wave-vectors have non-zero $x$ and $y$ components, the ratio
between the two being equal to that of the independent entries of the shear
matrix. More importantly, it can be seen that the same condition is satisfied
by the axisymmetric modes usually considered in the derivation of the magneto-
rotational instability. $\sigma_{ij}k^{i}k^{j}=0$ (but $\sigma_{ij}k^{j}\neq
0$), noting that such modes will always exist—even in the more general case
(considered below) where $\text{det}(\boldsymbol{\sigma})\neq 0$. It follows
that for such modes

${\omega}^{2}=-\frac{1}{6}\text{Tr}(\boldsymbol{\sigma}^{2})\Longrightarrow{\omega}=\pm
i\sqrt{\frac{1}{6}\text{Tr}(\boldsymbol{\sigma}^{2})}\;.$ (44)

These modes are also non-propagating, and half of them are unstable with a
(constant) growth rate about a factor of 2 smaller than in the previous case.
As the dispersion relation is quadratic (in the sound-proof limit), we can
explicitly solve it and confirm the expectation (and well known fact) that
shearing flows are generically unstable. Given the role played by
$\text{Tr}(\boldsymbol{\sigma}^{2})$ here and in what follows, let us briefly
comment on how one can intuitively see why this makes sense. First, let us
note that the same term enters the growth of entropy in the generalized second
law equation. Second, we also note that the same term when evaluated on the
usual background for the magneto-rotational instability would lead to
$\text{Tr}(\boldsymbol{\sigma}^{2})=4(\text{d}\Omega/\text{dln}R)^{2}$. As
such a term enters the usual instability dispersion relation and criterion, it
is not surprising it plays a pivotal role also for the present discussion.

Let us now build on this and discuss how vorticity and shear viscosity impact
on the generic instability of shearing flows. First consider the case where
the background has negligible vorticity but non-vanishing shear viscosity,
observing that the corresponding modes in a homogeneous background are stable
provided $\eta>0$. If viscosity is of microphysical origin, then $\eta>0$
follows from the second law of thermodynamics (see Landau & Lifshitz, 1959;
Andersson & Comer, 2021). If the viscosity is instead due to filtering, a
positive value of $\eta$ corresponds to an eddy-type model where energy is
cascading to smaller/unresolved scales111111We also note that it is in
principle possible to have $\eta<0$ when this is not of micro-physical origin
but an effective viscosity instead. Negative values of $\eta$ would
correspond, in the simplest eddy-viscosity-type models to net energy source at
the resolved scales, and hence an inverse cascade from the unresolved ones
(see, e.g. section VIB in Celora et al. (2021)). (M. Lesieur, 2005; Berselli
et al., 2006; McDonough, 2014; Schmidt, 2015). Going back to the case with
both shear and viscosity (in the sound-proof limit), the dispersion relation
of modes such that $\sigma_{ij}k^{j}=0$ is (cf. eq. 89)

${\omega}^{2}+i\eta
k^{2}{\omega}-\frac{1}{4}\left[\eta^{2}(k^{2})^{2}-2\text{Tr}(\boldsymbol{\sigma}^{2})\right]=0\;,$
(45)

where $k=|\mathbf{k}|$. Assuming $\eta>0$, stability corresponds to

$\eta^{2}(k^{2})^{2}-2\text{Tr}(\boldsymbol{\sigma}^{2})>0\;.$ (46)

In essence, comparing this to eq. 43 we see that viscosity tends to stabilize
shear-unstable modes, with a larger impact at smaller scales. This makes
intuitive sense. Next, consider modes such that $\sigma_{ij}k^{i}k^{j}=0$,
which are solutions of

${\omega}^{2}+i\eta
k^{2}{\omega}-\frac{1}{12}\left[3\eta^{2}(k^{2})^{2}-2\text{Tr}(\boldsymbol{\sigma}^{2})\right]=0\;.$
(47)

As before, these modes—to be compared with their counterparts in eq. 44—are
also stable (assuming $\eta>0$) provided the last term in eq. 47 is negative.
This will be true when the wave-number is sufficiently large, and we have
verified that the same trend is true for generic wave-vectors. In essence, we
learn (as one may have expected) that shear viscosity generically slows the
growth rate of unstable shear modes, and stabilises modes with small enough
wavelengths.

Turning to the case where the background has non-negligible vorticity and
shear, and taking the sound-proof limit as before, we first observe that the
fastest growing modes encountered before, namely those characterized by
$\sigma_{ij}k^{j}=0$, are not guaranteed to exist anymore, as the determinant
of the shear matrix cannot in general be assumed to be negligible (see section
3.1). Should these modes exist, though, their dispersion relation would be
(cf. eq. 91)

${\omega}^{2}=-\frac{1}{2}\text{Tr}(\boldsymbol{\sigma}^{2})+(\hat{k}\cdot\mathbf{W})^{2}\\\
\Longrightarrow{\omega}=\pm
i\sqrt{\frac{1}{2}\text{Tr}(\boldsymbol{\sigma}^{2})-(\hat{k}\cdot\mathbf{W})^{2}}\;,\quad\hat{k}=\mathbf{k}/k$
(48)

and we see, comparing this to eq. 43, that vorticity tends to stabilize them.
We also observe that—in contrast to shear viscosity—vorticity affects all such
modes by reducing their growth rate in a way that does not depend on their
wave-number (although the direction of propagation is important). Next—and
also because the modes we just looked at may not exist—we consider modes such
that $\sigma_{ij}k^{i}k^{j}=0$, with dispersion relation

${\omega}^{2}=-\frac{1}{6}\text{Tr}(\boldsymbol{\sigma}^{2})+(\hat{k}\cdot\mathbf{W})^{2}\;.$
(49)

Comparing this to eq. 44, we observe again that vorticity tends to stabilize
such modes in a way that does not depend on the wave-number. We have verified
that the same trend is also true for generic wave-vectors. As a final point,
it is easy to verify that the case with only background vorticity is generally
stable (not only in the sound-proof limit).

In summary, a sheared background flow is generically unstable already at the
hydrodynamic level, which is a well-known fact. However, we have considered
the impact that shear viscosity and/or vorticity have on the instability of
the possible hydrodynamic modes. The results show that shear viscosity tends
to weaken the instability in general, with larger effects for larger wave-
numbers. Meanwhile, vorticity has a stabilizing effect which does not depend
on the wave-number. Finally, let us also point to appendix B where we show
that the general dispersion relation derived in appendix C (discussed here) is
shown to encompass the classic Rayleigh stability criterion.

##  6 Magneto-shear instability in the local frame

Having explored the hydrodynamic case, let us perturb the corresponding MHD
equations and study the impact of the magnetic field on the generic shear
instabilities we encountered. We consider a barotropic equation of state and
retain gradients in the background velocity only, as we want to focus on the
magneto-shear nature of the instability (cf. Hawley & Balbus, 1991; Shibata,
2015). The continuity equation is obviously unchanged, while the perturbed
Euler equation becomes

$\partial_{t}\delta v_{i}+\delta
v^{j}\nabla_{j}v_{i}+\frac{1}{\rho}\nabla_{i}\delta
P+\frac{1}{\mu_{0}\rho}\left[B_{j}\nabla_{i}\delta B^{j}-B^{j}\nabla_{j}\delta
B_{i}\right]=0\\\ \Longrightarrow-i\omega\delta
v_{i}+i\frac{c_{s}^{2}}{\rho}k_{i}\delta\rho+\frac{i}{\mu_{0}\rho}\left[(B_{j}\delta
B^{j})k_{i}-(B^{j}k_{j})\delta B_{i}\right]\\\ +\sigma_{ij}\delta
v^{j}+{\epsilon}_{ijk}W^{j}\delta v^{k}=0\;.$ (50)

Finally, the perturbed induction equation is

$\partial_{t}\delta B^{i}+B^{i}\nabla_{j}\delta v^{j}-B^{j}\nabla_{j}\delta
v^{i}-\delta B^{j}\nabla_{j}v^{i}+\delta B^{i}\nabla_{j}v^{j}=0\\\
\Longrightarrow-i\omega\delta B^{i}+iB^{i}(k_{j}\delta
v^{j})-i(B^{j}k_{j})\delta v^{i}\\\ -\sigma^{ij}\delta
B_{j}-{\epsilon}^{ijk}W_{k}\delta B_{j}+\frac{2}{3}\theta\delta B^{i}=0\;.$
(51)

We will now discuss the linearized system that follows from these equations,
as before focussing on the results and providing more detailed steps in
appendix D.

Let us first recap the mode analysis for the homogeneous case. In order to
derive the fully compressible dispersion relation, we first re-scale the
magnetic field as

$\mathbf{v}_{A}\doteq\frac{\mathbf{B}}{\sqrt{\mu_{0}\rho}}\;,\quad\delta\mathbf{v}_{A}\doteq\frac{\delta\mathbf{B}}{\sqrt{\mu_{0}\rho}}\;.$
(52)

Then, computing the dispersion relation from the linearized equations above
(neglecting background shear and vorticity) we obtain (cf. eq. 99)

$-{\omega}\left({\omega}^{2}-\left(\mathbf{v}_{A}\cdot\mathbf{k}\right)^{2}\right)\left[{\omega}^{4}-\left(v_{A}^{2}+c_{s}^{2}\right)k^{2}{\omega}^{2}+c_{s}^{2}k^{2}\left(\mathbf{v}_{A}\cdot\mathbf{k}\right)^{2}\right]=0\;,$
(53)

where the roots of the second factor correspond to Alfvén waves, while those
of the quartic polynomial in square brackets describe fast and slow magneto-
sonic waves. Before moving on to discuss the impact of shear and vorticity,
let us briefly note what happens to such modes when we take the sound-proof
limit—in which the speed of sound is large. From eq. 53 we see that in the
first case fast magneto-sonic waves are filtered out, while the slow ones
reduce to Alfvén waves. In the opposite limit—when disturbances are much
faster than the sound waves—the dispersion relation describes Alfvén waves and
the low-$c_{s}$ limit of fast magneto-sonic waves. This limit corresponds to
ignoring fluid pressure perturbations while retaining variations in the
magnetic pressure.

Let us now consider the impact of a non-negligible background shear. Again, we
focus on the results, with the details of the derivation provided in section
D.2. We start from the fully compressible dispersion relation and consider the
sound-proof limit, leading to the dispersion relation (103). As in the
hydrodynamic case considered earlier, we first consider the case
$\text{det}(\boldsymbol{\sigma})=0$, and look for modes such that
$\sigma^{ij}k_{j}=0$. It is then easy to see that the general dispersion
relation (103) simplifies to (ignoring a trivial root)

$\left[{\omega}^{2}-\left(\frac{1}{2}\text{Tr}(\boldsymbol{\sigma}^{2})-(\mathbf{v}_{A}\cdot\mathbf{k})^{2}\right)\right]^{2}=0\;.$
(54)

Comparing to the corresponding hydrodynamic modes in eq. 43, we immediately
see that the magnetic field tends to have a stabilizing effect (provided it is
not orthogonal to the wave-vector, in which case it has no effect whatsoever).

Next, we take (again, as before) $\text{det}(\boldsymbol{\sigma})=0$ and
consider modes such that $\sigma_{ij}k^{i}k^{j}=0$ (but $\sigma_{ij}k^{j}\neq
0$). The relevant dispersion relation can then be written (making use of eq.
86)

${\omega}^{4}+b_{2}{\omega}^{2}+b_{4}=0\;,$ (55)

with

$b_{2}=\frac{2}{3}\text{Tr}(\boldsymbol{\sigma}^{2})-2(\mathbf{v}_{A}\cdot\mathbf{k})^{2}\;,$
(56a) and
$b_{4}=\frac{1}{12}\text{Tr}(\boldsymbol{\sigma}^{2})^{2}-\frac{2}{3}\text{Tr}(\boldsymbol{\sigma}^{2})(\mathbf{v}_{A}\cdot\mathbf{k})^{2}+(\mathbf{v}_{A}\cdot\mathbf{k})^{4}\;.$
(56b)

The stabilizing effect of the magnetic field is evident from fig. 1, where
both the frequency and $|\mathbf{v}_{A}\cdot\mathbf{k}|$ are plotted in units
of $\sqrt{\text{Tr}(\boldsymbol{\sigma}^{2})}$.

Figure 1:  Real and imaginary part of the solutions of eq. 55, with the
frequency and $|\mathbf{v}_{A}\cdot\mathbf{k}|$ in units of
$\sqrt{\text{Tr}(\boldsymbol{\sigma}^{2})}$. The solutions plotted correspond
to the fastest growing modes evolving on top of an MHD sheared background. We
see that the magnetic field has a stabilizing effect, as the growth rates are
reduced with respect to those of the corresponding hydrodynamic modes. The
stabilizing effect is all the more pronounced the more the wave-vector is
aligned with the magnetic field lines, and is switched off for modes
propagating in the directions perpendicular to the magnetic field lines. In
particular, modes corresponding to sufficiently large values of
$|\mathbf{v}_{A}\cdot\mathbf{k}|$ are turned stable.

The key point here is that, while the background shear is required for the
instability (the vanishing-shear modes are stable Alfvèn waves in the sound-
proof limit), the magnetic field is not the main driver. This is evident from
the results as the imaginary part of the unstable modes remains finite in the
limit $\mathbf{v}_{A}\to 0$, and the limiting value coincides with the
hydrodynamic result (from the previous section). This observation, possibly
unexpected at first sight, deserves a thorough discussion, and we will return
to this issue in section 7.1. Before we expand on that aspect, let us stress
that the results make intuitive sense. The magnetic field impacts on the
instability in that it breaks the hydrodynamic isotropy and dampens the growth
of unstable modes propagating along magnetic field lines. This also suggests
that shear-instability driven turbulence is isotropic in the hydrodynamic case
but inherently anisotropic for magnetized flows, consistently with the overall
picture (see e.g. Schekochihin, 2020; Beresnyak, 2019). Before moving on, it
is also worth noting that the background velocity profile considered by Balbus
& Hawley (Balbus & Hawley, 1991, 1998) is characterized by having a shear
matrix with vanishing determinant (and expansion rate), and also that for
axisymmetric modes $\sigma_{ij}k^{i}k^{j}=0$, while the fastest growing MRI
modes propagate vertically, with $\sigma_{ij}k^{j}=0$.

As a final comment before we make contact with the usual MRI and the Rayleigh
criterion, we have also considered the case with non-negligible background
vorticity only. The details are discussed in section D.3, but the crucial
result is that magnetized flows are generically stable in this case. As this
feature is unchanged from the corresponding fluid case, it is reasonable to
expect that the same trend we discussed for the purely hydrodynamical case
will also apply to the magnetized case with both shear and vorticity:
Vorticity tends to stabilize shear-unstable modes in a manner independent of
the wave number (although the orientation of the wave-vector with respect to
the vorticity is expected to have an impact.).

##  7 Concluding Remarks: The instability in perspective

We set out with the intention of discussing the magneto-rotational instability
in a general background, relaxing the symmetry constraints associated with the
standard analysis and possibly deriving an instability “criterion” relevant
for (highly) dynamical environments and nonlinear simulations. In fact, the
usual criterion is often used to compute a quality factor for the instability,
and hence measure how well the mechanism “could be resolved” in a given
simulation. However, having set up the analysis (and the required tools) in an
arguably sensible way, we arrived at results which were not in line with the
“naïve” expectations. Given this, it makes sense to comment on the
implications. Moreover, we need to highlight an important “missing ingredient”
in the discussion; the need to involve some suitable filtering operation to
make the discussion sensible in the first place. We will deal with each of
these questions in turn, starting with the implications of our results for the
MRI.

###  7.1 The MRI vs the Rayleigh criterion

A key aspect of the MRI is that adding a weak magnetic field on top of a
hydrodynamically stable shearing flow changes the nature of the problem and
makes it unstable. In discussing this problem, however, it is often
“forgotten” that the relevant hydrodynamic stability criterion (Rayleigh,
1917) guarantees stability only for axisymmetric modes (cf. the discussion in
appendix B). Adding a magnetic field renders such modes unstable—technically,
the non-axisymmetric ones are not (Balbus & Hawley, 1992). Thus it is clear
that the MRI is relevant only in situations where we can think of axisymmetric
modes as being “preferred” in some sense. An immediate example of this is an
accretion disk, which involves a globally axisymmetric background for the
perturbations. This then immediately tells us that applying the results to the
dynamical context of neutron star mergers is a much more subtle endeavour. In
fact, this exercise is problematic from the outset.

To back up this claim, we show in appendix A that we can reproduce the MRI
perturbation equations and dispersion relation through the local frame
construction. However, for the specific MRI calculation there exists a
preferred local frame: the co-rotating frame. This local frame is set up
considering an observer that is co-rotating with the fluid along some orbit,
and the coordinate axes rotate in such a way that one of them always points in
the radial direction of a global cylindrical coordinate system. Another
coordinate axis always points in the azimuthal direction. This local frame is
“preferred” as the axes are (by construction) tied to those of the most
natural global coordinate system. In a sense, we could set up different local
co-rotating observers and construct the global axes by stitching together the
local ones. In the case of a general and truly local analysis, however, this
additional piece of information is not available.

Moreover, we show in appendix B how one may set up (for the circular and
axisymmetric background flow) a local frame that is “co-moving but not co-
rotating” with the fluid. In doing so, we derived the corresponding dispersion
relation, confirmed that the result is consistent with the general formulae
(cf. appendix C), and showed how we can recover the usual Rayleigh criterion
(and hence also the MRI criterion) as long as we perform the conversion to the
relevant co-rotating frame frequency.

These arguments clarify the sense in which the MRI (and similarly the Rayleigh
stability) is a “global instability analyzed with local tools”. The local
analysis needs to be “augmented” by pieces of information that cannot be truly
local. The upshot of this is that, in a merger-like scenario (where
assumptions regarding the global properties of the flow are debatable) we
should probably not expect the standard instability criteria to provide a
faithful indication/diagnostic of what is actually going on. The standard
argument will apply, but only if there is a meaningful sense of (Rayleigh
stable) flow on a scale larger than that at which the plane-wave analysis is
carried out. This would make a discussion much more difficult for any given
numerical simulation, but so be it.

###  7.2 The missing ingredient: Filtering

Throughout the discussion we have focussed on the analytical development,
sweeping issues associated with actual numerical data “under the carpet”. The
key issue here is that we ignored the question of how one would, in practice,
construct the background suitable for the perturbation analysis given
nonlinear simulation dynamics. If the numerical data shows some degree of
axisymmetry, then the split could be achieved via an azimuthal averaging (and
the usual criteria will apply). If this is not the case, however, we need to
apply some suitable filtering operation to remove small scale fluctuations
from a gradually varying “background”. With such a filtering procedure at hand
we could start from real numerical data and perform the split into background
plus fluctuations. The background data could then be used as numerical input
in the dispersion relations we have discussed to assess whether or not the
instability is active and how rapidly it grows. In a nonlinear setting this
split is (obviously) not guaranteed to make sense. Suppose that the
instability we are trying to uncover acts on some characteristic scale $L$,
say. Then we need a background that varies on a larger scale than this,
otherwise the notion of a shear flow that becomes unstable due to smaller
scale waves makes no sense. This argument relies on an explicit filtering
step, separating the instability scale $L$ from the variation of the
background. The construction of such a filter should be possible, at least in
principle, in many situations (see, for example, Celora et al. (2021)). Of
course, the scale separation may not apply in actual problems of interest. In
addition, the analysis in this work (like any linear stability analysis)
relies on the expansion detailed in section 2 and the neglect of non-linear
fluctuations. The intention is to describe the early stages of the development
of the instability. Crucially, because the split into background and
fluctuations depends on the filtering procedure, the requirement that
fluctuations are indeed small will further constrain the procedure (e.g. the
filter width).

Further complicating the discussion is the unavoidable implicit filtering
associated with the finite numerical resolution. We know from the large body
of work on turbulence simulations that sub-grid dynamics may play an important
role in a robust description of the dynamics. This typically involves a
suitable large-eddy scheme to represent the subgrid dynamics. Hence, the
analysis involves elements of choice (effectively, the closure relations).
Crucially, the effective field theory that is/should be simulated is not that
of the ideal theory. All current models—both the ones discussed in Carrasco et
al. (2020) and Radice (2020) as well as the covariant scheme of Celora et al.
(2021)—modify the principal part of the equations of motion. Therefore the
analysis of the model “that is actually solved” is fundamentally changed, even
when the closure terms are small. In essence, an instability analysis of
numerical simulation data needs to consider the impact of an effective
viscosity/resistivity. Given the presently available tools, we do not have a
particularly good handle on this issue. We are forced to conclude that we also
need to make progress on the development of robust large-eddy models before we
can make a sensible attempt to demonstrate the presence of the MRI in a highly
dynamical environment.

## Acknowledgements

NA and IH gratefully acknowledge support from Science and Technology Facility
Council (STFC) via grant numbers ST/R00045X/1 and ST/V000551/1.

## Data availability

Supplementary data underlying this article are available at the following
links: https://doi.org/10.5281/zenodo.7612469 and
https://doi.org/10.5281/zenodo.7612479.

## References

  * Andersson (2012) Andersson N., 2012, Phys. Rev. D, 86, 043002
  * Andersson & Comer (2021) Andersson N., Comer G. L., 2021, Living Rev. Relativ., 24
  * Andersson et al. (2021) Andersson N., Hawke I., Celora T., Comer G. L., 2021, Mon. Not. R. Astron. Soc., 509, 3737
  * Anile (1989) Anile A. M., 1989, Relativistic fluids and magneto-fluids. Cambridge University Press
  * Balbus & Hawley (1991) Balbus S. A., Hawley J. F., 1991, Astrophys. J., 376, 214
  * Balbus & Hawley (1992) Balbus S. A., Hawley J. F., 1992, Astrophys. J., 400, 610
  * Balbus & Hawley (1998) Balbus S. A., Hawley J. F., 1998, Rev. Mod. Phys., 70, 1
  * Barletta (2022) Barletta A., 2022, Mech. Res. Commun., 124, 103939
  * Barrow et al. (2007) Barrow J. D., Maartens R., Tsagas C. G., 2007, Phys. Rept., 449, 131
  * Beresnyak (2019) Beresnyak A., 2019, Living Rev. Comput. Astrophys., 5
  * Berselli et al. (2006) Berselli L., Iliescu T., Layton W., 2006, Mathematics of Large Eddy Simulation of Turbulent Flows. Springer, Heidelberg
  * Carrasco et al. (2020) Carrasco F., Viganò D., Palenzuela C., 2020, Phys. Rev. D, 101, 063003
  * Carroll (2019) Carroll S. M., 2019, Spacetime and Geometry: An Introduction to General Relativity. Cambridge University Press
  * Celora et al. (2021) Celora T., Andersson N., Hawke I., Comer G. L., 2021, Phys. Rev. D, 104, 084090
  * Chandrasekhar (1960) Chandrasekhar S., 1960, Proc. Natl. Acad. Sci., 46, 253
  * Drazin & Reid (2004) Drazin P., Reid W., 2004, Hydrodynamic Stability. Cambridge Mathematical Library, Cambridge University Press
  * Duez et al. (2006) Duez M. D., Liu Y. T., Shapiro S. L., Shibata M., Stephens B. C., 2006, Phys. Rev. Lett., 96, 031101
  * Ellis (1973) Ellis G. F. R., 1973, Cargèse Lect. Phys., 6, 1
  * Goldreich & Lynden-Bell (1965) Goldreich P., Lynden-Bell D., 1965, Mon. Not. R. Astron. Soc., 130, 125
  * Goodman & Xu (1994) Goodman J., Xu G., 1994, Astrophys. J., 432, 213
  * Gourgoulhon (2013) Gourgoulhon E., 2013, Special Relativity in General Frames. Springer-Verlag Berlin Heidelberg
  * Guilet et al. (2022) Guilet J., Reboul-Salze A., Raynaud R., Bugli M., Gallet B., 2022, Mon. Not. R. Astron. Soc., 516, 4346
  * Hawley & Balbus (1991) Hawley J. F., Balbus S. A., 1991, Astrophys. J., 376, 223
  * Hawley et al. (1995) Hawley J. F., Gammie C. F., Balbus S. A., 1995, Astrophys. J., 440, 742
  * Hawley et al. (2011) Hawley J. F., Guan X., Krolik J. H., 2011, Astrophys. J., 738, 84
  * Hawley et al. (2013) Hawley J. F., Richers S. A., Guan X., Krolik J. H., 2013, Astrophys. J., 772, 102
  * Hayashi et al. (2022) Hayashi K., Kiuchi K., Kyutoku K., Sekiguchi Y., Shibata M., 2022, arXiv:2211.07158
  * Held & Mamatsashvili (2022) Held L. E., Mamatsashvili G., 2022, Mon. Not. R. Astron. Soc., 517, 2309
  * Hill (1878) Hill G. W., 1878, Am. J. Math., 1, 5
  * Hoshino (2015) Hoshino M., 2015, Phys. Rev. Lett., 114, 061101
  * Inchingolo et al. (2018) Inchingolo G., Grismayer T., Loureiro N. F., Fonseca R. A., Silva L. O., 2018, Astrophys. J., 859, 149
  * Kempski et al. (2019) Kempski P., Quataert E., Squire J., Kunz M. W., 2019, Mon. Not. R. Astron. Soc., 486, 4013
  * Kiuchi et al. (2018) Kiuchi K., Kyutoku K., Sekiguchi Y., Shibata M., 2018, Phys. Rev. D, 97, 124039
  * Kiuchi et al. (2022) Kiuchi K., Fujibayashi S., Hayashi K., Kyutoku K., Sekiguchi Y., Shibata M., 2022,  arXiv:2211.07637
  * Korn & Korn (2013) Korn G., Korn T., 2013, Mathematical Handbook for Scientists and Engineers. Dover Civil and Mechanical Engineering, Dover Publications
  * Landau & Lifshitz (1959) Landau L. D., Lifshitz E. M., 1959, Fluid mechanics. Pergamon Press
  * M. Lesieur (2005) M. Lesieur O. Métais P. C., 2005, Large-Eddy Simulations of Turbulence. Cambridge University Press
  * Mahajan & Krishan (2008) Mahajan S. M., Krishan V., 2008, Astrophys. J., 682, 602
  * Margalit et al. (2022) Margalit B., Jermyn A. S., Metzger B. D., Roberts L. F., Quataert E., 2022, Astrophys. J., 939, 51
  * McDonough (2014) McDonough J., 2014, Introductory Lectures on Turbulence. CreateSpace Independent Publishing Platform, Scotts Valley
  * Misner et al. (2017) Misner C., Thorne K., Wheeler J., Kaiser D., 2017, Gravitation. Princeton University Press
  * Palapanidis (2018) Palapanidis K., 2018, PhD thesis, University of Southampton
  * Palenzuela et al. (2022) Palenzuela C., Aguilera-Miret R., Carrasco F., Ciolfi R., Kalinani J. V., Kastaun W., Miñano B., Viganò D., 2022, Phys. Rev. D, 106, 023013
  * Radice (2020) Radice D., 2020, Symmetry, 12, 1249
  * Rayleigh (1917) Rayleigh L., 1917, Proc. R. Soc. Lond. A, 93, 148
  * Schekochihin (2020) Schekochihin A. A., 2020,  arXiv:2010.00699
  * Schmidt (2015) Schmidt W., 2015, Living Rev. Comput. Astrophys., 1, 2
  * Shakura & Sunyaev (1973) Shakura N. I., Sunyaev R. A., 1973, Astron. Astrophys., 24, 337
  * Shakura et al. (2022) Shakura N., Postnov K., Kolesnikov D., Lipunova G., 2022, arXiv:2210.15337
  * Sharma et al. (2006) Sharma P., Hammett G. W., Quataert E., Stone J. M., 2006, Astrophys. J., 637, 952
  * Shibata (2015) Shibata M., 2015, Numerical Relativity. World Scientific
  * Siegel et al. (2013) Siegel D. M., Ciolfi R., Harte A. I., Rezzolla L., 2013, Phys. Rev. D, 87, 121302
  * Thorne & Blandford (2017) Thorne K. S., Blandford R. D., 2017, Modern Classical Physics. Princeton University Press
  * Vasil et al. (2013) Vasil G. M., Lecoanet D., Brown B. P., Wood T. S., Zweibel E. G., 2013, Astrophys. J., 773, 169
  * Velikhov (1959) Velikhov E., 1959, Sov. Phys. JETP, 36, 995
  * Zier & Springel (2022) Zier O., Springel V., 2022, Mon. Not. R. Astron. Soc., 517, 2639

##  Appendix A Formulating the MRI in the local frame

In this appendix we derive the MRI using the local frame construction
discussed in the main text. We consider the circular velocity profile assumed
in Balbus & Hawley (1991), $\mathbf{v}=v^{\hat{\varphi}}\hat{\varphi}$ with
$v^{\hat{\varphi}}={\Omega}(R)R$ where we use cylindrical coordinates and an
orthonormal basis on the “tangent space” (as usual). We then distinguish
between indices with a “hat” corresponding to the orthonormal basis, and those
without that correspond to the coordinate basis. We then pick an orbit at some
radial distance $R_{0}$ and choose an observer that is co-rotating with
angular frequency identical to that of the background flow at $R_{0}$, that is
$\mathbf{v}_{obs}=\Omega_{0}R\hat{\varphi}$ where
${\Omega}_{0}={\Omega}(R_{0})$. The observer is then accelerated with
acceleration $\mathbf{a}=-{\Omega}_{0}^{2}R\hat{R}$, and the velocity of the
fluid with respect to such an observer then is
$\mathbf{v}^{\prime}=({\Omega}-{\Omega}_{0})R\hat{\varphi}$. We then set up
the axes of the observer’s local frame so that one is pointing in the radial
direction ($\hat{e}_{1}$), one is pointing in the azimuthal direction
($\hat{e}_{2}$) and the third one is aligned with the rotation axis
($\hat{e}_{3}$). Introducing coordinates associated with this observer, we can
then write the background fluid velocity as

$\mathbf{v}^{\prime}=\frac{\text{d}{\Omega}}{\text{dln}R}\bigg{|}_{R_{0}}x^{\prime}\hat{e}_{2}+{\mathcal{O}}(x^{\prime
2})\;.$ (57)

We have neglected terms of order ${\mathcal{O}}\left(x^{\prime 2}\right)$ as
we will only need the velocity and its gradients evaluated at the origin of
the frame—so that such terms will not enter the perturbation equations anyway.
Computing the gradients we then obtain

$\partial^{\prime}_{i}v^{\prime}_{j}=\begin{pmatrix}0&s_{0}&0\\\ 0&0&0\\\
0&0&0\end{pmatrix}\;,\qquad
s_{0}=\frac{\text{d}{\Omega}}{\text{dln}R}\bigg{|}_{R_{0}}\;.$ (58)

As the local frame of the observer is rotating with angular velocity
${\Omega}_{0}\hat{e}_{3}$, we need to include the Coriolis force into the
perturbation equations. We then write the perturbed Euler and continuity
equations (dropping the primes for clarity, and retaining only gradients in
the background velocity) as

$\partial_{t}\delta\rho+\rho\partial_{i}\delta v^{i}=0\;,$ (59a)
$\partial_{t}\delta v_{i}+2{\Omega}_{0}{\epsilon}_{i3k}\delta v^{k}+\delta
v^{j}\partial_{j}v_{i}+\frac{c_{s}^{2}}{\rho}\partial_{i}\delta\rho\\\
+\frac{1}{\mu_{0}\rho}\left[B_{j}\partial_{i}\delta
B^{j}-B^{j}\partial_{j}\delta B_{i}\right]=0\;,$ (59b)

and, introducing a WKB plane-wave expansion,

$-i{\omega}\delta\rho+i\rho k_{i}\delta v^{i}=0\;,$ (60a) $-i{\omega}\delta
v_{i}+2{\Omega}_{0}{\epsilon}_{i3k}\delta v^{k}+s_{0}\delta_{i2}\delta
v^{1}+i\frac{c_{s}^{2}}{\rho}k_{i}\delta\rho\\\
+\frac{i}{\mu_{0}\rho}\left[B_{j}k_{i}\delta B^{j}-B^{j}k_{j}\delta
B_{i}\right]=0\;,$ (60b)

Next, focus on the induction equation. As we have discussed above, the
induction equation in the co-rotating frame retains the inertial form. We then
have

$\delta\left[\partial_{j}\left(v^{j}B^{i}-v^{i}B^{j}\right)\right]=\delta
v^{j}\partial_{j}B^{i}+v^{j}\partial_{j}\delta B^{i}-B^{j}\partial_{j}\delta
v^{i}-\delta B^{j}\partial_{j}v^{i}\;,$ (61)

where we made use of the no-monopoles constraint, the vanishing expansion of
the background flow and got rid of the divergence of the perturbed velocity
consistently with the Boussinesq approximation. Introducing the WKB plane-wave
expansion and evaluating the background quantities at the origin of the local
frame we then end up with

$-i{\omega}\delta B^{i}-iB^{j}k_{j}\delta v^{i}-\delta
B^{1}s_{0}\delta^{i2}=0\;.$ (62)

In eqs. 60 and 62 we recognize the terms entering the perturbation equations
in Balbus & Hawley (1991) (with the exception of gradients in the background
pressure that we are neglecting here). We also note that we do not need to
formally neglect terms of the form $B/R$ as these terms do not appear in the
explicit local frame construction. We conclude by noting that, at the special
relativistic level, a uniformly rotating observer and the co-rotating one are
not the same as the latter is also accelerated (see Gourgoulhon, 2013, ch.
13). However, this difference is irrelevant at the level of the Newtonian
perturbation equations since i) pseudo-acceleration terms drop out of the
perturbed Euler equation ii) non-inertial terms in the induction equation
involving the four-acceleration are negligible in the Newtonian limit.

###  A.1 Another look at the non-inertial MHD equations

Before we move on to take a closer look at the Rayleigh criterion, let us show
how the terms involving the local frame rotation drop out of the induction
equation. Even though we have already argued this happens in general (cf. eq.
39), we here prove this for the specific case of a co-rotating observer. We do
so as this allows us to appreciate better why the cancellation comes about. We
will use a notation that is common in general relativity, that is the notion
of spin coefficients associated with a non-coordinate basis (Carroll, 2019).
The covariant derivative of a tensor $T^{\hat{a}}_{\hat{b}}$ is

$\nabla_{a}T^{\hat{a}}_{\hat{b}}=\partial_{a}T^{\hat{a}}_{\hat{b}}+{\omega}_{a\;\hat{c}}^{\;\hat{a}}T^{\hat{c}}_{\hat{b}}-{\omega}_{a\;\hat{b}}^{\;\hat{c}}T^{\hat{a}}_{\hat{c}}\;,$
(63)

with

${\omega}_{a\;\hat{c}}^{\;\hat{b}}=e^{\hat{c}}_{d}\,e^{e}_{\hat{b}}\,\Gamma^{d}_{ae}-e^{d}_{\hat{b}}\partial_{a}e^{\hat{c}}_{d}\;,$
(64)

where $\Gamma^{a}_{bc}$ are the connection coefficients associated with the
coordinates chosen while $e^{\hat{a}}_{b}$ is the matrix connecting the
coordinate basis to the orthonormal one. We also note here that for a local
frame we identify
$u^{a}{\omega}_{a\;\hat{b}}^{\;\hat{a}}={\Omega}^{\hat{a}}_{\;\hat{b}}$
introduced in eq. 37. We then introduce (Born) coordinates associated with a
uniformly rotating observer (axes suitably oriented so that the angular
velocity is ${\Omega}_{0}\hat{z}$)

$t^{\prime}=t\;,\quad z^{\prime}=z\;,\quad
x^{\prime}=R\text{cos}({\Omega}_{0}t+\varphi)\;,\quad
y^{\prime}=R\text{sin}({\Omega}_{0}t+{\varphi})\;,$ (65)

where primed coordinates are Cartesian (i.e. non-rotating). Computing the spin
coefficients (starting from a flat metric) we then obtain

${\omega}_{\varphi\;\hat{\varphi}}^{\;\hat{R}}=-1\;,\quad{\omega}_{\varphi\;\hat{R}}^{\;\hat{\varphi}}=+1\;,\quad{\omega}_{t\;\hat{\varphi}}^{\;\hat{R}}=-{\Omega}_{0}\;,\quad{\omega}_{t\;\hat{R}}^{\;\hat{\varphi}}={\Omega}_{0}\;,$
(66)

showing that, as the coordinates “mix space and time” we need to introduce a
covariant derivative in the time-direction as well. We then write the non-
inertial induction equation as (cf. eq. 33)

$\nabla_{t}B^{\hat{i}}+\nabla_{\hat{j}}\left(v^{\hat{j}}B^{\hat{i}}-v^{\hat{i}}B^{\hat{j}}\right)+{\epsilon}^{\hat{i}\hat{j}\hat{k}}B_{\hat{j}}{\Omega}^{0}_{\hat{k}}=0\;.$
(67)

It is then easy to verify, by means of eq. 66 that

$\nabla_{t}B^{\hat{i}}+{\epsilon}^{\hat{i}\hat{j}\hat{k}}B_{\hat{j}}{\Omega}^{0}_{\hat{k}}=\partial_{t}B^{\hat{i}}\;,$
(68)

thus confirming the result in eq. 39 and the use of the inertial induction
equations. We stress that the co-rotating frame rotation vector is
$\Omega_{0}\hat{e}_{3}$, and is the same (by construction) as the vorticity of
the observer. This is why we see the same cancellation as in eq. 39, where we
assumed they are equal.

##  Appendix B A closer look at the Rayleigh stability criterion

The key point of the magneto-rotational instability is that the circular
velocity background is stable against axisymmetric hydrodynamic perturbations,
while adding a (however weak) magnetic field changes completely the nature of
the system and makes it unstable to such perturbations. We now revisit the
Rayleigh criterion in order to highlight the key role played by the co-
rotating observer in deriving the criterion. These important aspects have to
be kept in mind when looking at the general results derived below (cf.
appendices C and D) and discussed in sections 6 and 5.

Starting from eq. 60, and ignoring terms associated with the magnetic field,
we write the coefficient matrix (ordering the perturbed quantities as
$\\{\delta\rho/\rho,\delta v^{1},\delta v^{2},\delta v^{3}\\}$)

$\begin{pmatrix}-{\omega}&k_{1}&k_{2}&k_{3}\\\
c_{s}^{2}k_{1}&-{\omega}&2i{\Omega}_{0}&0\\\
c_{s}^{2}k_{2}&-i\frac{\kappa^{2}}{2{\Omega}_{0}}&-{\omega}&0\\\
c_{s}^{2}k_{3}&0&0&-{\omega}\end{pmatrix}\;,$ (69)

where

$\frac{\kappa^{2}}{2{\Omega}_{0}}=2{\Omega}_{0}+s_{0}=2{\Omega}_{0}+\frac{\text{d}{\Omega}}{\text{dln}R}\bigg{|}_{R_{0}}\;,$
(70)

and the dispersion relation reads

${\omega}^{4}-\left(c_{s}^{2}k^{2}-\kappa^{2}\right){\omega}^{2}+ic_{s}^{2}s_{0}k_{1}k_{2}{\omega}+c_{s}^{2}\kappa^{2}(k_{3})^{2}=0\;.$
(71)

Taking the sound-proof limit we then end up with

$k^{2}{\omega}^{2}-is_{0}k_{1}k_{2}{\omega}-\kappa^{2}k_{3}^{2}=0\;.$ (72)

From this we easily see that, if we assume axisymmetric perturbations, namely
$k_{2}=0$, we obtain the usual Rayleigh stability criterion, that is
$\kappa^{2}>0$ (Rayleigh, 1917). We stress that, as is well-known, the
criterion does not guarantee that non-axisymmetric modes are stable. In fact,
rewriting the dispersion relation in terms of $\Delta=-i{\omega}$ and using
the Routh-Hurwitz criterion (Korn & Korn, 2013) we find that, on top of the
Rayleigh criterion we would also need

$s_{0}k_{1}k_{2}\leq 0\;.$ (73)

We also note that, the story changes if we take the opposite limit, in which
case the Rayleigh criterion is sufficient to guarantee stability of non-
axisymmetric perturbations as well. This would also be the case had we assumed
incompressibility from the start.

Having discussed the usual Rayleigh criterion using the co-rotating observer,
we now re-work it using an observer that is orbiting with the fluid at a given
orbital distance $R_{0}$ but whose (local frame) axes are non rotating. We do
this for two reasons. First, it will allow for a direct comparison with the
general results in appendix C. Second, we have argued that choosing to work
with a rotating or non-rotating observer is, in general, just a matter of
taste. We then pick up an orbit $R_{0}$ as before, and choose the observer to
be co-orbiting with the background flow at the specific orbit

$\mathbf{v}_{obs}=-{\Omega}_{0}y_{0}\hat{x}+{\Omega}_{0}x_{0}\hat{y}\;,$ (74)

where we used global Cartesian coordinates and
$x_{0}(t)=R_{0}\text{cos}({\Omega}_{0}t),y_{0}(t)=R_{0}\text{sin}({\Omega}_{0}t)$
describe the worldline of the observer (the origin of the axes is suitably
chosen so that $z_{0}(t)=0$). The background fluid velocity is then

$\mathbf{v}=-{\Omega}x\hat{y}+{\Omega}y\hat{x}\;,\qquad{\Omega}={\Omega}(\sqrt{x^{2}+y^{2}})\;,$
(75)

so that, considering the relative velocity
$\mathbf{v}^{\prime}=\mathbf{v}-\mathbf{v}_{obs}$ and expanding around
$(x_{0},y_{0})$ we obtain

$\mathbf{v}^{\prime}=-\left[s_{0}\frac{x_{0}y_{0}}{x_{0}^{2}+y_{0}^{2}}x^{\prime}+\left({\Omega}_{0}+s_{0}\frac{y_{0}^{2}}{x_{0}^{2}+y_{0}^{2}}\right)y^{\prime}\right]\hat{x}\\\
+\left[\left({\Omega}_{0}+s_{0}\frac{x_{0}^{2}}{x_{0}^{2}+y_{0}^{2}}\right)x^{\prime}+s_{0}\frac{x_{0}y_{0}}{x_{0}^{2}+y_{0}^{2}}x^{\prime}\right]\hat{y}\;,$
(76)

where $x^{\prime}=x-x_{0},\,y^{\prime}=y-y_{0}$. We can now choose a local
region around a specific point $(x_{0},y_{0},z_{0})$ on the orbit and choose
to re-orient the axes by a constant rotation so that the observer velocity is
moving only in the $y-$direction. We then set up the local frame in such a way
that the local axes are non-rotating and oriented like the global cartesian
ones. We can therefore write the gradients as

$\partial^{\prime}_{i}v^{\prime}_{j}=\begin{pmatrix}0&{\Omega}_{0}+s_{0}&0\\\
-{\Omega}_{0}&0&0\\\ 0&0&0\end{pmatrix}\;,$ (77)

and the coefficient matrix of the linearized Euler plus continuity system is
(cf. eq. 60 and ignore both magnetic field terms and the Coriolis force as the
axis are non-rotating)

$\begin{pmatrix}-{\omega}&k_{1}&k_{2}&k_{3}\\\
c_{s}^{2}k_{1}&-{\omega}&i{\Omega}_{0}&0\\\
c_{s}^{2}k_{2}&-i({\Omega}_{0}+s_{0})&-{\omega}&0\\\
c_{s}^{2}k_{3}&0&0&-{\omega}\end{pmatrix}\;.$ (78)

We can then compute the dispersion relation to find

${\omega}^{4}-\left[c_{s}^{2}k^{2}+{\Omega}_{0}({\Omega}_{0}+s_{0})\right]{\omega}^{2}+ic_{s}^{2}k_{1}k_{2}s_{0}\,{\omega}\\\
+c_{s}^{2}{\Omega}_{0}({\Omega}_{0}+s_{0})(k_{3})^{2}=0\;,$ (79)

and observe this is consistent with the general dispersion relation in eq. 91
when restricted to the shear and vorticity associated with eq. 77. However,
this is not quite the dispersion relation we obtained above. The reason for
this is that the two local observers we have considered measure different
frequencies, as the axes of the co-rotating observer rotate with angular
velocity ${\Omega}_{0}\hat{e}_{3}$ with respect to the other. To show why this
is the resolution to the apparent conflict, let us consider once again the
Born coordinates (cf. eqs. 65 and 66). Given any vector $a^{\hat{i}}$ we have

$\nabla_{t}a^{\hat{i}}=\partial_{t}\hat{a}^{i}+{\Omega}_{0}{\epsilon}^{\hat{i}\hat{3}\hat{k}}a_{\hat{k}}\;.$
(80)

This relation, when we introduce a plane-wave WKB expansion translates to

$-i{\omega}_{rot}\delta a^{\hat{i}}=-i{\omega}_{nr}\delta
a^{\hat{i}}+{\Omega}_{0}{\epsilon}^{\hat{i}\hat{3}\hat{k}}\delta
a_{\hat{k}}\;,$ (81)

where ${\omega}_{rot}$ is the frequency measured by the co-rotating observer,
while ${\omega}_{nr}$ is the frequency measured by an observer that has the
same worldline but uses non-rotating axes. Specifying eq. 81 to the perturbed
velocity (noting that it would not apply to the continuity equation as the
density is a scalar), and noting that the frequency in eq. 78 corresponds to
${\omega}_{nr}$, we can reconcile the results obtained from eq. 78 with those
from eq. 69. We also note here that the same logic applies when we consider
magnetized flows. That is, if we work with the inertial induction equation and
compute the background velocity gradients as in eq. 77, we also need to take
into account the relation in eq. 81 for magnetic field disturbances to get
back to eq. 62 and the MRI dispersion relation.

##  Appendix C Hydrodynamic dispersion relations

In this appendix we provide more details of the results discussed in section 5
of the main text. In order to derive the dispersion relation and study the
effects of a sheared background, it is convenient to choose a basis that is
adapted to it. Because the shear is a symmetric trace-free matrix there exists
a basis such that

$\sigma_{ij}=\text{diag}\left(\sigma_{1},\sigma_{2},-(\sigma_{1}+\sigma_{2})\right)\;.$
(82)

We will make use of this basis to write down the coefficients matrix of the
linearized system. Before doing so, however, it is reasonable to wonder
whether this change of basis has any impact on the perturbation equations. We
are, in fact always free to choose a basis in the tangent space that is not
associated with the coordinates chosen, but this (in general) introduces
additional terms in the covariant derivative. Let us spell out why this is not
the case here. Working with a non-coordinate basis, we need to account for
spin-coefficients—defined as in eq. 64—when a derivative acts on vectors and
tensors. The first term in eq. 64 then vanishes as we are here working with a
non-rotating Cartesian frame (so that the Christoffel symbols vanish), while
the second term is in general non vanishing, and accounts for the fact that
the change of basis (in the tangent space) needed to diagonalize the shear
matrix may change from point to point. In the context of this analysis,
however, we are looking at scales smaller than those over which background
quantities vary. In essence, also this second term vanishes as the shear
matrix is (by construction) constant over the local region of fluid we are
zooming on.

Working in the shear adapted basis (cf. eq. 82), we write the coefficient
matrix of the linearized system as (cf. eq. 42)

$\begin{pmatrix}-{\omega}&\rho k_{1}&\rho k_{2}&\rho k_{3}\\\
\frac{c_{s}^{2}}{\rho}k_{1}&-{\omega}-i\sigma_{1}-i\eta
L_{1}&iW^{3}-\frac{i}{6}\eta k_{1}k_{2}&-iW^{2}-\frac{i}{6}\eta k_{1}k_{3}\\\
\frac{c_{s}^{2}}{\rho}k_{2}&-iW^{3}-\frac{i}{6}\eta
k_{2}k_{1}&-{\omega}-i\sigma_{2}-i\eta L_{2}&iW^{1}-\frac{i}{6}\eta
k_{2}k_{3}\\\ \frac{c_{s}^{2}}{\rho}k_{3}&iW^{2}-\frac{i}{6}\eta
k_{3}k_{1}&-iW^{1}-\frac{i}{6}\eta k_{3}k_{2}&-{\omega}+i\sigma_{3}-i\eta
L_{3}\end{pmatrix}\;,$ (83)

where $\sigma_{3}=\sigma_{1}+\sigma_{2}$,
$L_{1}=\frac{2}{3}k_{1}^{2}+\frac{1}{2}k_{2}^{2}+\frac{1}{2}k_{3}^{2}$ and
$L_{2},\,L_{3}$ are similarly defined.

The dispersion relation is computed taking the determinant of this matrix and
equating it to zero. In order to keep the discussion as general as possible
(i.e. without having to refer to a specific background configuration) we will
decompose the coefficients of the characteristic polynomial in terms of
scalars built from background quantities. In the simplest cases this can be
done “by eye”, but the procedure can easily become quite messy. The logic is
nonetheless simple: we group the different terms in each coefficients
according to the power of the various background quantities, for example we
group all the terms quadratic in the shear and wave-vector components. We then
build all the possible scalars that are quadratic in the shear and wave-
vector, and look for the correct linear combination of them. This logic can be
easily implemented on a computer algebra program like Mathematica121212See
https://doi.org/10.5281/zenodo.7612469 for more details.. We now discuss the
dispersion relations obtained by retaining only shear terms, both shear and
viscous terms, and lastly shear and vorticity terms.

Before doing so, we first observe that the coefficients will involve scalars
constructed from the shear matrix only. As with any $3\crossproduct 3$ matrix,
the shear matrix $\boldsymbol{\sigma}$ has three invariants

$I_{1}=\text{Tr}(\boldsymbol{\sigma})\;,\quad
I_{2}=\frac{1}{2}\left[\text{Tr}(\boldsymbol{\sigma}^{2})-\left(\text{Tr}(\boldsymbol{\sigma})\right)^{2}\right]\;,\quad
I_{3}=\text{det}(\boldsymbol{\sigma})\;,$ (84)

related via the Cayley-Hamilton theorem as

$\boldsymbol{\sigma}^{3}-I_{1}\boldsymbol{\sigma}^{2}+I_{2}\boldsymbol{\sigma}-I_{3}\mathbb{I}=0\;,$
(85)

where $\mathbb{I}$ is the $3\crossproduct 3$ identity matrix. Because the
shear matrix is trace-free, we will write the coefficients in terms of
$\nicefrac{{1}}{{2}}\text{Tr}(\boldsymbol{\sigma}^{2})$ and
$\text{det}(\boldsymbol{\sigma})$. It is also useful to note that there will
always exist modes such that $\sigma_{ij}k^{i}k^{j}=0$. In the convenient
shear basis, these are characterized by $k^{1}=k^{2}=k^{3}$ if the determinant
is not vanishing (that is $s_{1}\neq s_{2}$), and $k^{1}=k^{2}$ when it does.
It follows that for such modes

$-\frac{1}{2}\text{Tr}(\boldsymbol{\sigma}^{2})+\sigma^{2}_{ij}\hat{k}^{i}\hat{k}^{j}=-\frac{1}{6}\text{Tr}(\boldsymbol{\sigma}^{2})\;,\quad\hat{k}=\mathbf{k}/|\mathbf{k}|=\mathbf{k}/k\;.$
(86)

The resulting dispersion relation for the case with negligible vorticity and
viscosity is then

${\omega}^{4}+a_{2}{\omega}^{2}+a_{1}{\omega}+a_{0}=0$ (87)

with

$\displaystyle a_{2}$
$\displaystyle=-c_{s}^{2}k^{2}+\frac{1}{2}\text{Tr}(\boldsymbol{\sigma}^{2})\;,$
(88a) $\displaystyle a_{1}$
$\displaystyle=i\left[c_{s}^{2}\sigma_{ij}k^{i}k^{j}-\text{det}(\sigma)\right]\;,$
(88b) $\displaystyle a_{0}$
$\displaystyle=c_{s}^{2}\left[\sigma^{2}_{ij}k^{i}k^{j}-\frac{1}{2}\text{Tr}(\boldsymbol{\sigma}^{2})k^{2}\right]\;,$
(88c)

whose solutions are discussed in the main text (cf. section 5).

Turning to the case with non-vanishing shear viscosity, the dispersion
relation is then

${\omega}^{4}+a_{3}{\omega}^{3}+a_{2}{\omega}^{2}+a_{1}{\omega}+a_{0}=0\;,$
(89)

with

$\displaystyle a_{3}$ $\displaystyle=\frac{5}{3}i\eta k^{2}\;,$ (90a)
$\displaystyle a_{2}$
$\displaystyle=-c_{s}^{2}k^{2}+\frac{1}{2}\text{Tr}(\boldsymbol{\sigma}^{2})-\frac{1}{12}\eta\left[11\eta(k^{2})^{2}-2\sigma_{ij}k^{i}k^{j}\right]\;,$
(90b) $\displaystyle a_{0}$
$\displaystyle=c_{s}^{2}\left[\sigma^{2}_{ij}k^{i}k^{j}-\frac{1}{2}\text{Tr}(\boldsymbol{\sigma}^{2})k^{2}-\frac{1}{2}\eta
k^{2}\sigma_{ij}k^{i}k^{j}+\frac{1}{4}\eta^{2}(k^{2})^{3}\right]\;,$ (90c) and
$a_{1}=ic_{s}^{2}\left[\sigma_{ij}k^{i}k^{j}-\eta(k^{2})^{2}\right]\\\
+i\Big{\\{}-\frac{1}{6}\left[\sigma^{2}_{ij}k^{i}k^{j}-2\text{Tr}(\boldsymbol{\sigma}^{2})k^{2}\right]\\\
+\frac{1}{12}\eta^{2}k^{2}(\sigma_{ij}k^{i}k^{j})-\frac{1}{6}\eta^{3}(k^{2})^{3}-\text{det}(\boldsymbol{\sigma})\Big{\\}}\;.$
(90d)

We can sanity check this dispersion relation by considering the homogeneous
background limit. It is then immediate to see that, in the sound-proof limit,
this would be stable provided $\eta>0$. It turns out that this condition
guarantees stability even outside of the sound-proof limit, as can be verified
by means of the Routh-Hurwitz criterion (Korn & Korn, 2013). As discussed in
the main text, this is the case both when shear viscosity is of micro-physical
origin (due to the second law of thermodynamics) but also when this is an
effective viscosity due to filtering provided this models a net energy
transfer to smaller unresolved scales—which is intuitive.

Finally, the dispersion relation in the case where the background has non-
negligible vorticity and shear is

${\omega}^{4}+a_{2}{\omega}^{2}+a_{1}{\omega}+a_{0}=0\;,$ (91)

with

$\displaystyle a_{2}$
$\displaystyle=-c_{s}^{2}k^{2}-W^{2}+\frac{1}{2}\text{Tr}(\boldsymbol{\sigma}^{2})\;,$
(92a) $\displaystyle a_{1}$
$\displaystyle=i\left[c_{s}^{2}\sigma_{ij}k^{i}k^{j}-\text{det}(\boldsymbol{\sigma})-\sigma_{ij}W^{i}W^{j}\right]\;,$
(92b) $\displaystyle a_{0}$
$\displaystyle=c_{s}^{2}\left[\sigma^{2}_{ij}k^{i}k^{j}-\frac{1}{2}\text{Tr}(\boldsymbol{\sigma}^{2})k^{2}+(\mathbf{k}\cdot\mathbf{W})^{2}\right]\;.$
(92c)

##  Appendix D MHD dispersion relations

In this appendix we provide more details on the results discussed in section 6
of the main text. In order to keep the presentation tidy, we will discuss
separately the homogeneous background case (cf. section D.1), the sheared case
(cf. section D.2) and the case with non-negligible background vorticity (cf.
section D.3).

###  D.1 Homogeneous background

In order to derive the fully compressible dispersion relation for the
homogeneous case, we first re-scale the magnetic field as in eq. 52, and
introduce a convenient basis $\\{\hat{v}_{A},\,\hat{q},\,\hat{s}\\}$ where
$\hat{v}_{A}=\mathbf{v}_{A}/|\mathbf{v}_{A}|$ while $\hat{q},\hat{s}$ complete
it to an orthonormal basis. For instance, assuming $\mathbf{v}_{A}$ is not
aligned with $\mathbf{k}$ we can construct it as

$\mathbf{q}=\mathbf{k}-\left(\mathbf{k}\cdot\hat{v}_{A}\right)\hat{v}_{A}\,,\qquad\hat{q}=\frac{\mathbf{q}}{|\mathbf{q}|}\,,\qquad\hat{s}=\hat{v}_{A}\crossproduct\hat{q}\;,$
(93)

so that131313If the wave-vector is along the background magnetic field we just
have to set $k^{q}=0$ in the following.

$\mathbf{k}=k^{v_{A}}\hat{v_{A}}+k^{q}\hat{q}\;.$ (94)

The coefficient matrix of the linearized system can then be written as (cf.
eqs. 41, 50 and 51 and ignore background vorticity and shear)

$\mathbf{M}=\begin{pmatrix}\mathbf{A}&\mathbf{C}\\\
\mathbf{C}^{\top}&\mathbf{D}\end{pmatrix}$ (95)

with

$\mathbf{A}=\begin{pmatrix}-{\omega}&\rho k_{v_{A}}&\rho k_{q}&0\\\
\frac{c_{s}^{2}}{\rho}k_{v_{A}}&-{\omega}&0&0\\\
\frac{c_{s}^{2}}{\rho}k_{q}&0&-{\omega}&0\\\ 0&0&0&-{\omega}\end{pmatrix}\;,$
(96a) $\mathbf{C}=\begin{pmatrix}0&0&0\\\ 0&0&0\\\ bk_{q}&-bk_{v_{A}}&0\\\
0&0&-bk_{v_{A}}\end{pmatrix}\;,$ (96b) and
$\mathbf{D}=\begin{pmatrix}-{\omega}&0&0\\\ 0&-{\omega}&0\\\
0&0&-{\omega}\end{pmatrix}\;.$ (96c)

As $\mathbf{D}$ is clearly invertible, we can reduce $\mathbf{M}$ into factors
via the Schur complement

$\begin{pmatrix}\mathbf{A}&\mathbf{C}\\\
\mathbf{C}^{\top}&\mathbf{D}\end{pmatrix}\\\
=\begin{pmatrix}\mathbf{I}_{4}&\mathbf{C}\mathbf{D}^{-1}\\\
\mathbf{0}_{4\times
3}&\mathbf{I}_{3}\end{pmatrix}\begin{pmatrix}\mathbf{A}-\mathbf{C}\mathbf{D}^{-1}\mathbf{C}^{\top}&\mathbf{0}_{4\times
3}\\\ \mathbf{0}_{3\times
4}&\mathbf{D}\end{pmatrix}\begin{pmatrix}\mathbf{I}_{4}&\mathbf{0}_{4\times
3}\\\ \mathbf{D}^{-1}\mathbf{C}^{\top}&\mathbf{I}_{3}\end{pmatrix}\;,$ (97)

and then compute the determinant as

$\text{det}(\mathbf{M})=\text{det}(\mathbf{D})\,\text{det}(\mathbf{A}-\mathbf{C}\mathbf{D}^{-1}\mathbf{C}^{\top})\;.$
(98)

The resulting dispersion relation is

$-{\omega}\left({\omega}^{2}-\left(\mathbf{v}_{A}\cdot\mathbf{k}\right)^{2}\right)\left[{\omega}^{4}-\left(v_{A}^{2}+c_{s}^{2}\right)k^{2}{\omega}^{2}+c_{s}^{2}k^{2}\left(\mathbf{v}_{A}\cdot\mathbf{k}\right)^{2}\right]=0\;.$
(99)

###  D.2 Sheared Background

Let us now consider the case where the background vorticity is negligible
while shear terms are not. Re-scaling the magnetic field as in eq. 52 and
decomposing eqs. 41, 50 and 51 (ignoring vorticity terms) as well as
$\delta\mathbf{v}$ and $\delta\mathbf{v}_{A}$ in the shear-adapted basis, the
coefficient matrix of the linearized system of equations reads

$\mathbf{M}=\begin{pmatrix}\mathbf{A}&\mathbf{C}\\\
\mathbf{C}^{\top}&\mathbf{D}\end{pmatrix}\;,$ (100)

where

$\mathbf{A}=\begin{pmatrix}-{\omega}&\rho k_{1}&\rho k_{2}&\rho k_{3}\\\
\frac{c_{s}^{2}}{\rho}k_{1}&-{\omega}-i\sigma_{1}&0&0\\\
\frac{c_{s}^{2}}{\rho}k_{2}&0&-{\omega}-i\sigma_{2}&0\\\
\frac{c_{s}^{2}}{\rho}k_{3}&0&0&-{\omega}+i(\sigma_{1}+\sigma_{2})\end{pmatrix}\;,$
(101a) $\mathbf{D}=\begin{pmatrix}-{\omega}+i\sigma_{1}&0&0\\\
0&-{\omega}+i\sigma_{2}&0\\\
0&0&-{\omega}-i(\sigma_{1}+\sigma_{2})\end{pmatrix}\;,$ (101b)
$\mathbf{C}=\begin{pmatrix}0&0&0\\\ I_{1}&v_{A}^{2}k_{1}&v_{A}^{3}k_{1}\\\
v_{A}^{1}k_{2}&I_{2}&v_{A}^{3}k_{2}\\\
v_{A}^{1}k_{3}&v_{A}^{2}k_{3}&I_{3}\end{pmatrix}\;,$ (101c)

while

$I_{1}=v_{A}^{1}k_{1}-\left(\mathbf{v}_{A}\cdot\mathbf{k}\right)\;,$ (102)

and $I_{2},I_{3}$ are defined similarly.

In a similar fashion as for the hydrodynamic case considered above, we will
decompose the coefficients of the characteristic polynomial in terms of
scalars built from background quantities. As we might have expected, the
resulting dispersion relation is a complicated seventh-degree
polynomial141414See https://doi.org/10.5281/zenodo.7612479 for the details of
the scalar decomposition. (and we sanity-checked it reduces to the homogeneous
case when we set to vanish the shear terms). In order to learn something
useful out of it, we then consider the sound-proof limit and retain only terms
proportional to the speed of sound. We end up with the following dispersion
relation

$a_{5}{\omega}^{5}+a_{4}{\omega}^{4}+a_{3}{\omega}^{3}+a_{2}{\omega}^{2}+a_{1}{\omega}+a_{0}=0\;,$
(103)

with

$a_{0}=-i\bigg{\\{}\text{det}(\boldsymbol{\sigma})\Big{[}\sigma^{2}_{ij}k^{i}k^{j}-\frac{1}{2}\text{Tr}(\boldsymbol{\sigma}^{2})\Big{]}\\\
+(\mathbf{v}_{A}\cdot\mathbf{k})^{2}\Big{[}\text{det}(\boldsymbol{\sigma})k^{2}-\frac{1}{2}(\sigma_{ij}k^{i}k^{j})\text{Tr}(\boldsymbol{\sigma}^{2})\Big{]}+(\mathbf{v}_{A}\cdot\mathbf{k})^{4}\sigma_{ij}k^{i}k^{j}\bigg{\\}}\;,$
(104a)
$a_{1}=\bigg{\\{}(\mathbf{v}_{A}\cdot\mathbf{k})^{4}k^{2}+\left(\mathbf{v}_{A}\cdot\mathbf{k}\right)^{2}\left[\sigma^{2}_{ij}k^{i}k^{j}-\text{Tr}(\boldsymbol{\sigma}^{2})k^{2}\right]+\\\
\text{det}(\boldsymbol{\sigma})\left(\sigma_{ij}k^{i}k^{j}\right)+\frac{1}{2}\text{Tr}(\boldsymbol{\sigma}^{2})\left[\frac{1}{2}\text{Tr}(\boldsymbol{\sigma}^{2})k^{2}-\sigma^{2}_{ij}k^{i}k^{j}\right]\bigg{\\}}\;,$
(104b)
$a_{2}=i\Bigg{\\{}-\frac{1}{2}(\sigma_{ij}k^{i}k^{j})\text{Tr}(\boldsymbol{\sigma}^{2})+\text{det}(\boldsymbol{\sigma})k^{2}+2(\mathbf{v}_{A}\cdot\mathbf{k})^{2}(\sigma_{ij}k^{i}k^{j})\Bigg{\\}}\;,$
(104c)
$a_{3}=\left[\text{Tr}(\boldsymbol{\sigma}^{2})k^{2}-2(\mathbf{v}_{A}\cdot\mathbf{k})^{2}k^{2}-\sigma^{2}_{ij}k^{i}k^{j}\right]\;,$
(104d) $a_{4}=-i(\sigma_{ij}k^{i}k^{j})\;,$ (104e) $a_{5}=k^{2}\;.$ (104f)

###  D.3 Background with vorticity

We now turn to the case with non-negligible background vorticity (but
negligible shear). We re-scale the magnetic field as in eq. 52 and introduce a
convenient basis $\\{\hat{W},\,\hat{q},\,\hat{s}\\}$, where
$\hat{W}=\mathbf{W}/|\mathbf{W}|$ while $\hat{q},\hat{s}$ complete it to an
orthonormal basis. For instance, assuming $\mathbf{v}_{A}$ is not aligned with
$\mathbf{W}$ we can construct it as

$\mathbf{q}=\mathbf{v}_{A}-\left(\mathbf{v}_{A}\cdot\hat{W}\right)\hat{W}\,,\qquad\hat{q}=\frac{\mathbf{q}}{|\mathbf{q}|}\,,\qquad\hat{s}=\hat{W}\crossproduct\hat{q}\;,$
(105)

and the magnetic field151515Note that the definition of $\hat{q}$ changes when
the background magnetic field is aligned with the vorticity, even though in
what follows we would simply have to set $v_{A}^{q}=0$.

$\mathbf{v}_{A}=v_{A}^{W}\hat{W}+v_{A}^{q}\hat{q}\;.$ (106)

The coefficient matrix of the linearized system then is (cf. eqs. 41, 50 and
51 and ignore shear terms)

$\mathbf{M}=\begin{pmatrix}\mathbf{A}&\mathbf{C}\\\
\mathbf{C}^{\top}&\mathbf{D}\end{pmatrix}\;,$ (107)

with

$\mathbf{C}=\begin{pmatrix}0&0&0\\\ -v_{A}^{q}k^{q}&v_{A}^{q}k^{W}&0\\\
v_{A}^{W}k^{q}&-v_{A}^{W}k^{W}&0\\\
v_{A}^{W}k^{s}&v_{A}^{q}k^{s}&-\left(v_{A}^{W}k^{W}+v_{A}^{q}k^{q}\right)\end{pmatrix}\;,$
(108a) $\mathbf{A}=\begin{pmatrix}-{\omega}&\rho k^{W}&\rho k^{q}&\rho
k^{s}\\\ \frac{c_{s}^{2}}{\rho}k^{W}&-{\omega}&0&0\\\
\frac{c_{s}^{2}}{\rho}k^{q}&0&-{\omega}&+iW\\\
\frac{c_{s}^{2}}{\rho}k^{s}&0&-iW&-{\omega}\end{pmatrix}\;,$ (108b) and
$\mathbf{D}=\begin{pmatrix}-{\omega}&0&0\\\ 0&-{\omega}&-iW\\\
0&+iW&-{\omega}\end{pmatrix}$ (108c)

As for the sheared case, we compute the dispersion relation by taking the
determinant of this matrix and decompose each coefficient as a sum of
scalars161616More details on the decomposition can be found at
https://doi.org/10.5281/zenodo.7612479.. Having sanity-checked the result by
contrasting it against the homogeneous background dispersion relation, we take
the sound proof limit. The sound-proof dispersion relation can then be written
as

${\omega}^{4}+b_{2}{\omega}^{2}+b_{4}=0\;,$ (109)

with

$\displaystyle b_{2}$
$\displaystyle=-\left[W^{2}+(\hat{k}\cdot\mathbf{W})^{2}+2(\mathbf{v}_{A}\cdot\mathbf{k})^{2}\right]\;,$
(110a) $\displaystyle b_{4}$
$\displaystyle=\left[(\mathbf{v}_{A}\cdot\mathbf{k})^{2}+W^{2}\right]\left[(\hat{k}\cdot\mathbf{W})^{2}+(\mathbf{v}_{A}\cdot\mathbf{k})^{2}\right]\;.$
(110b)

As this is a particularly simple quartic polynomial, we can study the
stability of its roots analytically. Considering eq. 109 as an equation for
${\omega}^{2}$ and computing the discriminant we obtain

$\left[W^{2}-(\hat{k}\cdot\mathbf{W})^{2}\right]^{2}\geq 0$ (111)

so that ${\omega}^{2}$-roots are real. As complex roots of a real algebraic
polynomials occur in pairs of complex conjugates, complex ${\omega}^{2}$-roots
would correspond to an instability. In order to have stable roots though, we
also need other conditions to be met. We, in fact need $b_{2}<0$ and $b_{4}>0$
to make sure that the ${\omega}^{2}$-roots are real and positive, so that
${\omega}$-roots are real as well. As this is evidently the case, we conclude
that background vorticity—just like in the purely hydrodynamic case—does not
lead to an instability.