Our philosophy throughout this work is to consider the simplest and most general modifications to MRI evolution on the largest (MHD) scales due to kinetic physics. We anticipate (but cannot prove) that these are the most important kinetic modifications to the MRI. We thus focus on the development of a gyrotropic pressure anisotropy – i.e. a pressure tensor that differs in the directions parallel (
$p_{\Vert }$
) and perpendicular (
$p_{\bot }$
) to the magnetic-field lines, but that is unchanged by rotations about the field line. This pressure anisotropy,
$\unicode[STIX]{x0394}p\equiv p_{\bot }-p_{\Vert }$
, causes an additional stress in the momentum equation, which can nonlinearly affect the MRI modes at much lower amplitudes than occurs in standard MHD. The gyrotropic approximation is generally valid when the magnetic field varies on spatial and temporal scales much larger than the ion gyroradius and inverse gyrofrequency.

### 2.1 Basic equations and closure models

Our equations are obtained as follows. A small patch of an accretion disc, co-orbiting with a fiducial point
$R_{0}$
in the mid-plane of the unperturbed disc at an angular velocity
$\unicode[STIX]{x1D734}=\unicode[STIX]{x1D6FA}\hat{\boldsymbol{z}}$
, is represented in Cartesian coordinates with the
$x$
and
$y$
directions corresponding to the radial and azimuthal directions, respectively. Differential rotation is accounted for by including the Coriolis force and by imposing a background linear shear flow,
$\boldsymbol{U}_{0}=-Sx\hat{\boldsymbol{y}}$
, where
$S\equiv -\text{d}\unicode[STIX]{x1D6FA}/\text{d}\ln R_{0}>0$
is the shear frequency; Keplerian rotation yields
$S=(3/2)\unicode[STIX]{x1D6FA}$
. The evolutionary equations for the first three moments of the plasma distribution function are then (Chew, Goldberger & Low 1956; Kulsrud 1983; Schekochihin *et al.*
2010),

where
$\text{d}/\text{d}t\equiv \unicode[STIX]{x2202}/\unicode[STIX]{x2202}t+\boldsymbol{u}\boldsymbol{\cdot }\unicode[STIX]{x1D735}$
is the convective derivative. The velocity
$\boldsymbol{u}$
in (2.1)–(2.5) includes both the background shear flow
$\boldsymbol{U}_{0}$
and perturbations on top of this
$\unicode[STIX]{x1D6FF}\boldsymbol{u}$
(e.g. the MRI). The other symbols have their usual meanings:
$\unicode[STIX]{x1D70C}$
is the mass density,
$\boldsymbol{B}$
is the magnetic field,
$\unicode[STIX]{x1D708}_{c}$
is the particle collision frequency, while
$B\equiv |\boldsymbol{B}|$
and
$\hat{\boldsymbol{b}}\equiv \boldsymbol{B}/B$
denote the magnetic-field strength and direction (note that
$\unicode[STIX]{x1D735}\boldsymbol{\cdot }(\hat{\boldsymbol{b}}\hat{\boldsymbol{b}}B^{2})$
is simply
$\boldsymbol{B}\boldsymbol{\cdot }\unicode[STIX]{x1D735}\boldsymbol{B}$
). The pressures perpendicular and parallel to
$\hat{\boldsymbol{b}}$
are
$p_{\bot }$
and
$p_{\Vert }$
respectively, while
$q_{\bot }$
and
$q_{\Vert }$
denote the fluxes of perpendicular and parallel heat in the direction parallel to
$\hat{\boldsymbol{b}}$
. Note that
$p_{\bot }$
and
$p_{\Vert }$
in (2.2) should in principle be summed over both particle species (with separate pressures for each species), while
$\unicode[STIX]{x1D70C}$
and
$\boldsymbol{u}$
in (2.1)–(2.3) are the ion density and flow velocity. For simplicity, in this work we solve only the ion pressure equations (i.e.
$p_{\bot }=p_{\bot ,i}$
,
$p_{\Vert }=p_{\Vert ,i}$
), which is justified in the limit of cold electrons (as expected in RIAFs, e.g. Sharma *et al.*
2007). We have also neglected non-ideal corrections to the induction equation (2.3) (e.g. the Hall term), which is appropriate given our neglect of finite Larmor radius (FLR) effects in (2.2) (we will, however, include a hyper-resistivity term in equation (2.3) for numerical reasons; see § 2.2). For convenience, we define the dimensionless anisotropy
$\unicode[STIX]{x1D6E5}\equiv \unicode[STIX]{x0394}p/p_{0}$
, where
$p_{0}=(p_{\bot }+2p_{\Vert })/3$
is the total thermal pressure, as well as the ratio of thermal to magnetic pressure
$\unicode[STIX]{x1D6FD}\equiv 8\unicode[STIX]{x03C0}p_{0}/B^{2}$
, the Alfvén speed
$v_{A}=B/\sqrt{4\unicode[STIX]{x03C0}\unicode[STIX]{x1D70C}}$
and its
$\hat{\boldsymbol{z}}$
component
$v_{Az}=B_{z}/\sqrt{4\unicode[STIX]{x03C0}\unicode[STIX]{x1D70C}}$
, the sound speed
$c_{s}=\sqrt{p_{0}/\unicode[STIX]{x1D70C}}$
and the parallel sound speed
$c_{s\Vert }=\sqrt{p_{\Vert }/\unicode[STIX]{x1D70C}}$
. The double-dot notation used in (2.4)–(2.5) and throughout this work is
$\hat{\boldsymbol{b}}\hat{\boldsymbol{b}}\,\boldsymbol{ : }\,\unicode[STIX]{x1D735}\boldsymbol{u}\equiv \sum _{i,j}b_{i}b_{j}\unicode[STIX]{x2202}_{i}u_{j}$
.

In their present form, equations (2.1)–(2.5) are not closed, due to the presence of the unspecified heat fluxes,
$q_{\bot }$
and
$q_{\Vert }$
. These must be either specified using a closure scheme, neglected or solved for using the full kinetic equations. In this work we consider three closures for
$q_{\bot }$
and
$q_{\Vert }$
(or equivalently, three approximations to (2.4)–(2.5)). These will be seen to lead to quite different behaviour in solutions of (2.1)–(2.5). They are:

*Collisionless Landau-fluid closure.* Landau-fluid (LF) closures have been used extensively in the fusion community (Hammett & Perkins 1990; Hammett, Dorland & Perkins 1992; Snyder, Hammett & Dorland 1997) and, to a lesser degree, for astrophysical applications (S06; Sharma *et al.*
2007). They are particularly well suited for modelling collisionless (
$\unicode[STIX]{x1D708}_{c}=0$
) plasmas. In the LF closure, the heat fluxes,

are constructed to replicate the effects of linear Landau damping. Here
$\unicode[STIX]{x1D735}_{\Vert }\equiv \hat{\boldsymbol{b}}\boldsymbol{\cdot }\unicode[STIX]{x1D735}$
denotes the gradient parallel to the field and
$|k_{\Vert }|$
denotes the wavenumber parallel to the field, which must be considered as an operator because it appears in the denominator of (2.6)–(2.7). The forms of the heat fluxes in (2.6)–(2.7) accurately reproduce the true kinetic growth rates and frequencies for a variety of large-scale (MHD) modes, including the MRI (Q02; Sharma *et al.*
2003). We refer the reader to Snyder *et al.* (1997) and S06 for more information.

*Weakly collisional ‘Braginskii’ closure.* In the Braginskii regime,
$|\unicode[STIX]{x1D735}\boldsymbol{u}|\ll \unicode[STIX]{x1D708}_{c}$
(Braginskii 1965), the pressure anisotropy is strongly influenced by collisional relaxation. We thus neglect
$\text{d}_{t}p_{\bot }$
and
$\text{d}_{t}p_{\Vert }$
in (2.4)–(2.5) and balance the double-adiabatic production of pressure anisotropy (the
$\hat{\boldsymbol{b}}\hat{\boldsymbol{b}}\,\boldsymbol{ : }\,\unicode[STIX]{x1D735}\boldsymbol{u}$
and
$\unicode[STIX]{x1D735}\boldsymbol{\cdot }\boldsymbol{u}$
terms) against its collisional relaxation (the
$\unicode[STIX]{x1D708}_{c}\unicode[STIX]{x0394}p$
terms) to find

where we have also used the fact that
$\unicode[STIX]{x1D708}_{c}/|\unicode[STIX]{x1D735}\boldsymbol{u}|\gg 1$
implies
$\unicode[STIX]{x0394}p\ll p_{0}$
. (For
$\unicode[STIX]{x1D6FD}\gg 1$
, the
$\unicode[STIX]{x1D735}\boldsymbol{\cdot }\boldsymbol{u}$
term can also be neglected.) When inserted into the momentum equation (2.2), equation (2.8) has the form of an anisotropic viscous stress, and is thus referred to as ‘Braginskii viscosity’ (or Braginskii MHD for the full set of equations). Note that we have neglected heat fluxes in arriving at (2.8), a simplification that is rigorously obtained if
$\unicode[STIX]{x1D708}_{c}/|\unicode[STIX]{x1D735}\boldsymbol{u}|\gg \unicode[STIX]{x1D6FD}^{1/2}$
(the ‘high-collisionality’ regime). On the other hand, if
$\unicode[STIX]{x1D708}_{c}/|\unicode[STIX]{x1D735}\boldsymbol{u}|\ll \unicode[STIX]{x1D6FD}^{1/2}$
(the ‘moderate-collisionality’ regime), the heat fluxes are strong over the time scales of the motion (Mikhailovskii & Tsypin 1971), and their contribution to the (ion) pressure anisotropy must be retained. (See § B.4 for further discussion.) In this case, there is no simple closure that can be devised (e.g. see appendix B of Squire, Schekochihin & Quataert 2017) and it is usually easier to consider the full LF system.

*Double-adiabatic closure.* The double-adiabatic, or Chew–Goldberger–Low (CGL), closure (Chew *et al.*
1956) simply involves setting
$q_{\bot }=q_{\Vert }=0$
. This approximation is far from justified for subsonic motions in the high-
$\unicode[STIX]{x1D6FD}$
plasmas considered here; however, the closure is useful for comparison with the LF closure by virtue of its relative simplicity. It has also been employed in a variety of previous computational studies (e.g. S06; Kowal, Falceta-Gonçalves & Lazarian 2011, Santos-Lima *et al.*
2014), and so it is worthwhile to diagnose the model’s successes and limitations.

An important caveat for each of these approximations to (2.4)–(2.5) relates to plasma microinstabilities. For our purposes, given the focus of this work on the high-
$\unicode[STIX]{x1D6FD}$
regime, the most significant of these are the firehose instability (Rosenbluth 1956; Chandrasekhar, Kaufman & Watson 1958; Parker 1958; Yoon, Wu & de Assis 1993), which is excited if

and the mirror instability (Hasegawa 1969; Southwood & Kivelson 1993; Hellinger 2007), which is excited if

(There are corrections to these
$\unicode[STIX]{x1D6FD}$
-dependent thresholds that arise from particle resonances and depend on the specific form of the distribution function; see, e.g. Klein & Howes 2015.) Important aspects of these instabilities (e.g. their regularization at small scales or particle scattering in their nonlinear evolution) are not captured by the closures we employ here, and kinetic calculations are needed to correctly understand their saturation. There have been a variety of recent works in this vein (Hellinger & Trávníček 2008; Schekochihin *et al.*
2008; Rosin *et al.*
2011; Kunz, Schekochihin & Stone 2014*a*
; Hellinger *et al.*
2015; Rincon, Schekochihin & Cowley 2015; Riquelme, Quataert & Verscharen 2015; Sironi & Narayan 2015; Melville, Schekochihin & Kunz 2016; Riquelme, Quataert & Verscharen 2016), which have shown that these microinstabilities generally act to pin the pressure anisotropy at the marginal stability limits. Interestingly, when
$\unicode[STIX]{x1D6E5}$
is driven beyond the stability boundaries, the microinstabilities achieve this in two stages: first, while the microscale fluctuations are growing secularly, by increasing (if
$\text{d}_{t}B<0$
) or decreasing (if
$\text{d}_{t}B>0$
) the magnetic-field strength
$B$
(in other words, the small-scale fluctuations contribute to
$B$
); second, as the instabilities saturate, by enhancing the scattering of particles and thus increasing
$\unicode[STIX]{x1D708}_{c}$
in (2.4)–(2.5). While the time for the firehose instability at moderate
$\unicode[STIX]{x1D6FD}$
to saturate is essentially set by gyro-scale physics, and so might be considered as instantaneous in a fluid model, the mirror instability saturates on a time scales comparable to the turnover time of the large-scale motions driving the anisotropy (see § 6.3 and Kunz *et al.*
2014*a*
, Rincon *et al.*
2015, Riquelme *et al.*
2015, Melville *et al.*
2016). We also note that there are also various other kinetic instabilities that could be important, for instance, the ion-cyclotron instability, or electron instabilities. We do not consider these in detail because the mirror and firehose instabilities are thought to be the most relevant to the high-
$\unicode[STIX]{x1D6FD}$
, ion-dominated regime that is the focus of this work (see § 5 for further discussion).

In practice, because the primary effect in both the secular and scattering regimes is to limit
$\unicode[STIX]{x1D6E5}$
at the threshold boundaries, we model these effects as a ‘hard-wall’ limit on
$\unicode[STIX]{x1D6E5}$
, following prior work (S06; Sharma *et al.*
2007, Santos-Lima *et al.*
2014, Chandra *et al.*
2015, Foucart *et al.*
2015). This simply limits
$\unicode[STIX]{x1D6E5}$
to
$1/\unicode[STIX]{x1D6FD}$
or
$-2/\unicode[STIX]{x1D6FD}$
if the dynamics drives
$\unicode[STIX]{x1D6E5}$
across these boundaries.
^{3}
One should, however, be careful with this simple ‘limiter’ method not to inadvertently remove interesting physics from the model. For instance, the parallel firehose instability is destabilized at the same point
$\unicode[STIX]{x1D6E5}=-2/\unicode[STIX]{x1D6FD}$
as where the magnetic tension is nullified by
$\unicode[STIX]{x0394}p$
(indeed, this is the cause of the instability), which can have a strong influence on the largest scales (Squire, Quataert & Schekochihin 2016; Squire *et al.*
2017) and is captured in even our simplest 1-D models. For this reason, in considering azimuthal-field KMRI modes, we have run calculations both with and without a firehose limiter, seeing very similar dynamics in each case. Finally, we note that with finite scale separations between
$\unicode[STIX]{x1D6FA}_{i}$
and
$S$
, as is the case in numerical simulations (Riquelme *et al.*
2012; Hoshino 2015; Kunz *et al.*
2016), there can be significant overshoot of the pressure anisotropy beyond the limits (2.9) and (2.10) (Kunz *et al.*
2014*a*
). This overshoot may be important for the large-scale evolution (see § 3) but is probably not representative of what happens in real systems, which usually have a very large dynamic range between
$\unicode[STIX]{x1D6FA}_{i}$
and
$S$
. These effects are discussed in detail in § 6.

### 2.2 Computational methods

A number of different numerical methods are used to solve (2.1)–(2.7). For investigating the 1-D evolution of a channel mode, we use a pseudo-spectral method, with standard dealiasing and hyper-diffusion operators used to remove the energy just above the grid scale. A very similar numerical method, albeit on a 3-D Fourier grid, is used for the studies of parasitic modes. For studies of the fully nonlinear 3-D evolution, we use a modified version of the finite-difference code Zeus, as described in S06. For simulations in the weakly collisional regime, we solve the full LF system (2.1)–(2.7), so as to correctly capture the effects of the heat fluxes in the moderate- and high-collisionality regimes (see § B.4).

A few words are needed regarding the numerical treatment of heat fluxes in the LF model (2.6)–(2.7). In particular, the
$1/|k_{\Vert }|$
operator is numerically awkward, because it is not diagonal in either Fourier space or real space. We thus use the prescription of S06 and replace this by a pre-chosen
$k_{L}$
for the Zeus implementation and the parasitic-mode studies, while for the 1-D collisionless calculations we use
$1/|k_{\Vert }|=1/|k_{z}|$
(once the mode reaches larger amplitude this may somewhat underestimate the heat fluxes, since
$k_{\Vert }<k_{z}$
if the field lines are not straight). Following S06, we have checked that varying the choice of
$k_{L}$
within a reasonable range, or using the choice
$|k_{\Vert }|=|k_{z}|$
, does not significantly affect the dynamics.

We use the methods detailed in appendix A3 of S06 to limit the pressure anisotropy:
$\unicode[STIX]{x1D708}_{c}$
is instantaneously increased in (2.4)–(2.5) whenever
$\unicode[STIX]{x1D6E5}$
passes the limits (2.9) or (2.10). We do not make any distinction between the ‘secular growth’ and ‘particle scattering’ phases of microinstability evolution with this method (see discussion above, around (2.10), and § 6.3), and more study is needed to better understand the successes and limitations of this simple limiter approach.

### 2.3 General comparison of kinetic models

Before commencing with our analysis of the KMRI, we highlight in this section some of the key similarities and differences between the LF, Braginskii and CGL models, as well as that of standard MHD.

First, irrespective of the closure details, the general effect of pressure anisotropy is to modify the Lorentz force, either by enhancing (
$\unicode[STIX]{x0394}p>0$
) or reducing (
$\unicode[STIX]{x0394}p<0$
) the effective magnetic tension (see (2.2)). The same relative pressure anisotropy
$\unicode[STIX]{x1D6E5}=\unicode[STIX]{x0394}p/p_{0}$
will thus have a greater dynamical effect as
$\unicode[STIX]{x1D6FD}$
increases, because
$p_{0}$
increases compared to
$B$
. In contrast, the generation of
$\unicode[STIX]{x1D6E5}$
in a changing magnetic field,

(or
$\unicode[STIX]{x1D6E5}\sim \unicode[STIX]{x1D708}_{c}^{-1}\hat{\boldsymbol{b}}\hat{\boldsymbol{b}}\,\boldsymbol{ : }\,\unicode[STIX]{x1D735}\boldsymbol{u}$
for Braginskii), does not depend on
$\unicode[STIX]{x1D6FD}$
. Thus, the dynamics of higher-
$\unicode[STIX]{x1D6FD}$
plasmas is in general more strongly influenced by self-generated pressure anisotropies than is the dynamics of lower-
$\unicode[STIX]{x1D6FD}$
plasmas.

Secondly, the spatial form of
$\unicode[STIX]{x0394}p$
generated in a given time-dependent
$B$
differs significantly between the LF, CGL and Braginskii models. This
$\unicode[STIX]{x0394}p$
is important for the nonlinear behaviour of the collisionless MRI. The influence of a spatially constant
$\unicode[STIX]{x1D6E5}$
can be considered from an essentially linear standpoint: it simply acts to enhance or reduce the Lorentz force by the factor
$(1+\unicode[STIX]{x1D6FD}\unicode[STIX]{x1D6E5}/2)$
. In contrast, when the spatial variation in
$\unicode[STIX]{x1D6E5}$
is similar to its magnitude, its effect is inherently nonlinear. For the same reason, detailed conclusions about the expected change in the spatial shape of an MRI mode due to pressure anisotropy do depend on the regime of interest (and model) – for example, collisionless versus Braginskii (Squire *et al.*
2016), or low versus high
$\unicode[STIX]{x1D6FD}$
– rather than being generic consequences of any self-generated pressure anisotropy. For our purposes, one can consider high-collisionality Braginskii MHD (equation (2.8)) as the limiting model that develops large spatial variation in
$\unicode[STIX]{x1D6E5}$
(since
$\unicode[STIX]{x1D6E5}$
is tied directly to
$\text{d}_{t}B$
), while the high-
$\unicode[STIX]{x1D6FD}$
limit of the LF model is the opposite, developing large
$\unicode[STIX]{x1D6E5}$
with very little spatial variation.
^{4}
This is because collisions generically act to reduce the magnitude of
$\unicode[STIX]{x1D6E5}$
without affecting the spatial variation in
$p_{\bot }$
and
$p_{\Vert }$
, whereas heat fluxes act to reduce spatial variation in
$p_{\bot }$
and
$p_{\Vert }$
without affecting the spatial average of
$\unicode[STIX]{x1D6E5}$
.

This last statement warrants further explanation, given the complexity of (2.4)–(2.7). The effect of the LF heat fluxes (2.6)–(2.7) can be clarified if we assume
$\unicode[STIX]{x1D708}_{c}=0$
and
$\unicode[STIX]{x0394}p\ll p_{0}$
(the latter is always valid at high
$\unicode[STIX]{x1D6FD}$
), so that

Then, assuming
$\hat{\boldsymbol{b}}\boldsymbol{\cdot }\unicode[STIX]{x1D735}q_{\bot ,\Vert }\gg q_{\bot ,\Vert }\unicode[STIX]{x1D735}\boldsymbol{\cdot }\hat{\boldsymbol{b}}$
, which is valid when the perturbation to the background field is small (i.e. when the field lines are nearly straight), the contribution to the
$p_{\bot }$
and
$p_{\Vert }$
evolution equations has the form
$\unicode[STIX]{x2202}_{t}p_{\bot }\sim -\unicode[STIX]{x1D70C}c_{s}|k_{\Vert }|(p_{\bot }/\unicode[STIX]{x1D70C})$
and
$\unicode[STIX]{x2202}_{t}p_{\Vert }\sim -\unicode[STIX]{x1D70C}c_{s}|k_{\Vert }|(p_{\Vert }/\unicode[STIX]{x1D70C})$
, respectively. These terms, which model the effect of Landau damping (Snyder *et al.*
1997), act as a ‘scale-independent’ diffusion of
$p_{\bot }$
and
$p_{\Vert }$
, damping inhomogeneities
^{5}
over the sound-wave time scale
$|k_{\Vert }|c_{s}$
. This is important because
$|k_{\Vert }|c_{s}\sim \unicode[STIX]{x1D6FD}^{1/2}\unicode[STIX]{x1D6FE}_{\text{MRI}}$
(where
$\unicode[STIX]{x1D6FE}_{\text{MRI}}$
is the MRI growth rate), showing that for
$\unicode[STIX]{x1D6FD}\gg 1$
, the heat fluxes will rapidly erase spatial variation in
$\unicode[STIX]{x0394}p$
on the time scale that the MRI grows. As a result, the spatial variation in
$\unicode[STIX]{x0394}p$
will be dwarfed by its mean; i.e.
$\unicode[STIX]{x0394}p$
will be nearly spatially constant. A more thorough discussion of these ideas is given in appendix B, where we solve explicitly for the
$\unicode[STIX]{x0394}p$
that arises in an exponentially growing, spatially varying magnetic perturbation. This shows that the double-adiabatic model for
$p_{\bot ,\Vert }$
generates a
$\unicode[STIX]{x0394}p$
with spatial variation of the order of its mean, while the addition of LF heat fluxes decreases the spatial variation of
$\unicode[STIX]{x0394}p$
by a factor
$\unicode[STIX]{x1D6FD}^{1/2}$
while leaving the mean
$\unicode[STIX]{x0394}p$
unchanged.