1 Introduction
Accretion of matter onto a central massive object is one of the most spectacular astronomical phenomena. A number of theoretical and numerical studies of accretion flows have been conducted over past decades (see Balbus & Hawley (Reference Balbus and Hawley1998), Lesur (Reference Lesur2021) and references therein), including the discovery of momentum transport due to turbulence driven by magnetorotational instability (MRI; Balbus & Hawley Reference Balbus and Hawley1991). On the observational front, the Event Horizon Telescope (EHT) successfully captured an image of a radiating disc at M87 (EHT Collaboration 2019), opening the door to direct comparisons between models and observations. However, there are many unsolved questions in MRI-driven turbulence that are crucial for interpreting such observations. In this study, we focus on the energy partition between Alfvénic and slow-mode-like compressive fluctuations in collisional MRI turbulence. This is important in deciding the partition of heating between ions and electrons at dissipation scales, where plasma is collisionless (Kawazura et al. Reference Kawazura, Schekochihin, Barnes, TenBarge, Tong, Klein and Dorland2020).
In order to calculate the partition of energy between Alfvénic and compressive fluctuations in a numerical simulation, one must access scales small enough that these fluctuations become energetically decoupled. It is known that such decoupling is established at the scale where the reduced magnetohydrodynamics (RMHD) approximation, 
$k_{\|}/k_{\perp} \ll 1$ and 
$|\delta {\boldsymbol {B}}|/B_0 \sim |\delta {\boldsymbol {u}}|/v_A \ll 1$, is satisfied (Schekochihin et al. Reference Schekochihin, Cowley, Dorland, Hammett, Howes, Quataert and Tatsuno2009). Here, 
$\boldsymbol{B}$ denotes the magnetic field, 
$\boldsymbol{u}$ denotes the flow velocity, 
$\boldsymbol{k}$ denotes the wavenumber, the subscript 
$\|$ (
${\perp}$) denotes the component parallel (perpendicular) to the ambient magnetic field, the prefix 
$\delta$ and the subscript 0 denote the fluctuation and equilibrium fields, respectively, and 
$v_A$ denotes the Alfvén speed. The RMHD approximation is expected to be satisfied at sufficiently small scales in the inertial range, because the large-scale magnetic field serves as an effective mean field for the fluctuations at the smaller scales (Kraichnan Reference Kraichnan1965). Once the cascade reaches the RMHD range, the partition of Alfvénic and compressive fluctuations is maintained down to the ion Larmor scale (Schekochihin et al. Reference Schekochihin, Cowley, Dorland, Hammett, Howes, Quataert and Tatsuno2009).
While the partition of Alfvénic and compressive fluctuations has been studied in externally forced magnetohydrodynamics (MHD) turbulence (Cho & Lazarian Reference Cho and Lazarian2002, Reference Cho and Lazarian2003; Makwana & Yan Reference Makwana and Yan2020), none of the previous studies of MRI turbulence have investigated this problem. Previous MRI turbulence simulations have suggested that, in order to reach the RMHD range, significantly higher numerical resolution is necessary for MRI turbulence than for externally forced MHD turbulence, because there is non-local energy transfer (Lesur & Longaretti Reference Lesur and Longaretti2011), meaning that the injection range is broad in the Fourier space.
Here, instead of carrying out brute-force high-resolution MHD simulations, we study the partition of Alfvénic and compressive fluctuations by reducing the MHD equations to a more tractable form that is valid only when there is a mean magnetic field in approximately azimuthal direction. The presence of the near-azimuthal mean magnetic field is a natural consequence of the differential rotation of the disc; even when the system is initialized with a purely vertical magnetic field, MRI creates a radial magnetic field which will then be twisted in the azimuthal direction. Indeed, a predominantly azimuthal magnetic field is quite often seen both in local and global simulations of MRI turbulence (e.g. Suzuki & Inutsuka Reference Suzuki and Inutsuka2009, Reference Suzuki and Inutsuka2014). A statistical analysis of MRI turbulence in incompressible MHD also supports the presence of a near-azimuthal mean field (Zhdankin et al. Reference Zhdankin, Walker, Boldyrev and Lesur2017). We show that the RMHD approximation captures the fastest-growing MRI modes in such a system. We then simulate this type of MRI turbulence numerically and show that the compressive fluctuations carry almost twice as much energy flux as Alfvénic fluctuations at the small scales, where the two kinds of fluctuations are decoupled.
2 Model
We consider a local shearing-box approximation (Goldreich & Lynden-Bell Reference Goldreich and Lynden-Bell1965) for a plasma in Cartesian coordinates 
$(X, Y, Z)$ located at a fixed radius 
$r = r_0$ and rotating with an angular velocity 
${\boldsymbol {\varOmega }} = \varOmega \hat {{\boldsymbol {Z}}}$, where 
$X$, 
$Y$ and 
$Z$ correspond to the radial, azimuthal and vertical (rotation-axis) directions. The MHD equations in these conditions are
$$\begin{gather} \frac{{\partial} \rho}{{\partial} t} + {\boldsymbol{u}}\boldsymbol{\cdot}\boldsymbol{\nabla}\rho + {\boldsymbol{u}}_0 \boldsymbol{\cdot}\boldsymbol{\nabla}\rho ={-}\rho(\boldsymbol{\nabla}\boldsymbol{\cdot}{\boldsymbol{u}}), \end{gather}$$
$$\begin{gather}\rho\left( \frac{{\partial} }{{\partial} t} + {\boldsymbol{u}}\boldsymbol{\cdot}\boldsymbol{\nabla} + {\boldsymbol{u}}_0\boldsymbol{\cdot}\boldsymbol{\nabla} \right){\boldsymbol{u}} ={-}\boldsymbol{\nabla}\left( p + \frac{B^2}{8{\rm \pi}} \right) + \frac{{\boldsymbol{B}}\boldsymbol{\cdot}\boldsymbol{\nabla}{\boldsymbol{B}}}{4{\rm \pi}} - 2\rho{\boldsymbol{\varOmega}}\times{\boldsymbol{u}} - \rho{\boldsymbol{u}}\boldsymbol{\cdot}\boldsymbol{\nabla}{\boldsymbol{u}}_0 , \end{gather}$$
$$\begin{gather}\frac{{\partial} {\boldsymbol{B}}}{{\partial} t} + {\boldsymbol{u}}\boldsymbol{\cdot}\boldsymbol{\nabla}{\boldsymbol{B}} + {\boldsymbol{u}}_0\boldsymbol{\cdot}\boldsymbol{\nabla}{\boldsymbol{B}} + {\boldsymbol{B}}(\boldsymbol{\nabla}\boldsymbol{\cdot}{\boldsymbol{u}}) = {\boldsymbol{B}}\boldsymbol{\cdot}\boldsymbol{\nabla}{\boldsymbol{u}} + {\boldsymbol{B}}\boldsymbol{\cdot}\boldsymbol{\nabla}{\boldsymbol{u}}_0, \end{gather}$$
$$\begin{gather}\frac{{\partial} p}{{\partial} t} + {\boldsymbol{u}}\boldsymbol{\cdot}\boldsymbol{\nabla} p + {\boldsymbol{u}}_0\boldsymbol{\cdot}\boldsymbol{\nabla} p + \varGamma p\boldsymbol{\nabla}\boldsymbol{\cdot}{\boldsymbol{u}} = 0, \end{gather}$$
where 
$\rho$ is the mass density, 
${\boldsymbol {u}}$ is the fluid velocity, 
${\boldsymbol {B}}$ is the magnetic field, 
$p$ is the thermal pressure, 
${\boldsymbol {u}}_0 \equiv q {\boldsymbol {X}}\times {\boldsymbol {\varOmega }}$ is the background shear flow, 
$q$ is a shear rate and 
$\varGamma = 5/3$ is the specific heat ratio. Hereafter, we only consider Keplerian rotation (
$q = 3/2$) and call (2.1a)–(2.1d) the ‘full-MHD’ equations.
Balbus & Hawley (Reference Balbus and Hawley1992b) showed that the fastest-growing MRI modes have 
$k_Z \to \infty$ when the ambient magnetic field 
${\boldsymbol {B}}_0$ approaches the azimuthal direction. These fastest-growing modes also satisfy 
$k_{\|} v_A/\varOmega \simeq 1$. For a near-azimuthal 
${\boldsymbol {B}}_0$, 
$\hat {{\boldsymbol {Z}}}$ is almost perpendicular to 
${\boldsymbol {B}}_0$, meaning that the fastest-growing modes satisfy 
$k_{\|}/k_{\perp} \ll 1$. Therefore, if the fastest-growing modes decide the nature of MRI turbulence at the smaller scales, we can ignore the scales that are outside of the 
$k_{\|}/k_{\perp} \ll 1$ approximationFootnote 1 . This idea motivates us to impose the RMHD approximation on the full-MHD equations (2.1a)–(2.1d). We assume also that the magnetic perturbations are separated from the time-invariant and spatially uniform mean fields as 
${\boldsymbol {B}} = {\boldsymbol {B}}_0 + \delta {\boldsymbol {B}}$, where 
${\boldsymbol {B}}_0$ is taken to have finite 
$Y$ (azimuthal) and 
$Z$ (vertical) components but no 
$X$ (radial) component. We assume that the density and pressure are also separated into a constant background and perturbations as 
$\rho = \rho _0 + \delta \rho$ and 
$p = p_0 + \delta p$. The angle between 
$\hat {{\boldsymbol {Y}}}$ and 
${\boldsymbol {B}}_0$ is denoted by 
$\theta \ [0 \le \theta \le {\rm \pi};$ the same definition as Quataert et al. (Reference Quataert, Dorland and Hammett2002)]. In this study, we focus on the near-azimuthal background field, i.e. 
$\sin \theta \ll 1$. Then, we introduce a ‘tilted’ coordinate set 
$(x, y, z)$ in which the 
$z$-axis is aligned with 
${\boldsymbol {B}}_0$, and the 
$x$-axis is aligned with the 
$X$-axis (figure 1), i.e. 
$(x,y,z)$ is a rotation of 
$(X,Y,Z)$ by 
${\rm \pi} /2 - \theta$ about the 
$\hat {{\boldsymbol {X}}}$ axis. When 
$\sin \theta \ll 1$, 
$\hat {{\boldsymbol {z}}}$ and 
$\hat {{\boldsymbol {y}}}$ almost align with 
$\hat {{\boldsymbol {Y}}}$ and 
$-\hat {{\boldsymbol {Z}}}$, respectively. This tilted coordinate set is more convenient than the standard coordinate set 
$(X, Y, Z)$ because 
$k_z \sim k_{\|} \ll k_{\perp}$ is a key criterion for the decoupling of Alfvénic and compressive fluctuations. In the standard coordinate set, however, 
$k_{\|}$ and 
$k_{\perp}$ are more difficult to separate, both being a mixture of 
$k_Y$ and 
$k_Z$.
Figure 1. Schematic of the conventional coordinate system 
$(X, Y, Z)$ and our tilted coordinate system 
$(x, y, z)$.
Thus, we impose the RMHD ordering 
$k_{\|} \ll k_{\perp},\ {\boldsymbol {u}} \ll v_A$ and 
$\delta {\boldsymbol {B}} \ll {\boldsymbol {B}}_0$ on (2.1a)–(2.1d). We also assume 
$ {\partial }/ {\partial } t \sim \varOmega \sim k_{\|} v_A$ and a near-azimuthal 
${\boldsymbol {B}}_0$, i.e. 
$\sin \theta \ll 1$, with 
$\sin \theta$ the same order of smallness as 
$k_{\|}/k_{\perp}$ and 
$\delta B/B_0$. The assumed anisotropy between 
$k_{\|}$ and 
$k_{\perp}$ is motivated by the critical balance conjecture (Goldreich & Sridhar Reference Goldreich and Sridhar1995, Reference Goldreich and Sridhar1997)
\begin{equation} k_{\|} v_A \sim k_{\perp} u_{\perp}, \end{equation}
which physically means that the time scales of linear wave propagation along 
${\boldsymbol {B}}_0$ and nonlinear cascade in the plane perpendicular to 
${\boldsymbol {B}}_0$ are of the same orderFootnote 2 . Then, we obtain the RMHD equations with differential rotation (see Appendix A for the detailed derivation)
$$\begin{gather} \frac{\mathrm{d}\varPsi}{{\rm d}t} = v_A\frac{{\partial} \varPhi}{{\partial} z}, \end{gather}$$
$$\begin{gather}\frac{\mathrm{d}}{{\rm d}t}\boldsymbol{\nabla}_{\perp}^2\varPhi = v_A\boldsymbol{\nabla}_{\|}\boldsymbol{\nabla}_{\perp}^2\varPsi - 2\varOmega \frac{{\partial} u_{\|}}{{\partial} y}, \end{gather}$$
$$\begin{gather}\frac{{\rm d}u_{\|}}{{\rm d}t} = v_A^2\boldsymbol{\nabla}_{\|}\frac{\delta B_{\|}}{B_0} + (2 - q)\varOmega \frac{{\partial} \varPhi}{{\partial} y}, \end{gather}$$
$$\begin{gather}\frac{\mathrm{d}}{{\rm d}t}\left( 1 + \frac{v_A^2}{c_S^2} \right)\frac{\delta B_{\|}}{B_0} = \boldsymbol{\nabla}_{\|} u_{\|} + \frac{q\varOmega }{v_A}\frac{{\partial} \varPsi}{{\partial} y}, \end{gather}$$
where 
$c_S$ is the sound speed, and 
$\varPhi$ and 
$\varPsi$ are the streamfunction and magnetic flux function defined by 
${\boldsymbol {u}}_{\perp} = \hat {{\boldsymbol {z}}}\times \boldsymbol {\nabla }_{\perp}\varPhi$ and 
$\delta {\boldsymbol {B}}_{\perp} = \sqrt {4{\rm \pi} \rho _0}\hat {{\boldsymbol {z}}}\times \boldsymbol {\nabla }_{\perp}\varPsi$, respectively. We have also defined 
$\mathrm {d}/\mathrm {d} t \equiv {\partial }/ {\partial } t + {\boldsymbol {u}}_\perp \boldsymbol {\cdot }\boldsymbol {\nabla }_\perp$ and 
$\boldsymbol {\nabla }_{\|} \equiv {\partial }/ {\partial } z + (\delta {\boldsymbol {B}}_{\perp}/B_0)\boldsymbol {\cdot }\boldsymbol {\nabla }_{\perp}$. Hereafter, we call these equations rotating RMHD (RRMHD)Footnote 3 . When 
$\varOmega = 0$, these become the standard RMHD equations (Kadomtsev & Pogutse Reference Kadomtsev and Pogutse1974; Strauss Reference Strauss1976), which is a long-wavelength limit of gyrokinetics and in which Alfvénic and compressive fluctuations are decoupled (Schekochihin et al. Reference Schekochihin, Cowley, Dorland, Hammett, Howes, Quataert and Tatsuno2009).
One may notice that (2.3a)–(2.3d) do not have the shearing effect that originates from 
${\boldsymbol {u}}_0\boldsymbol {\cdot }\boldsymbol {\nabla }$ terms in (2.1a)–(2.1d). This is due to 
$k_{\|}/k_{\perp} \ll 1$ and 
$\sin \theta \ll 1$; in a shearing box, the radial wavenumber depends on time as 
$k_x(t) = k_x(0) + q\varOmega t (k_y\sin \theta + k_{\|} \cos \theta )$ (see, e.g. the fourth term on the left-hand side of (A2)), and the time-dependent term on the right-hand side is ordered out because 
$\sin \theta \sim k_{\|}/k_x(0) \sim \epsilon$. However, when we consider a long-time evolution 
$\varOmega t \sim \epsilon ^{-1}$, the time dependence is not negligible. In that case, the non-modal growth of MRI (Squire & Bhattacharjee Reference Squire and Bhattacharjee2014a,Reference Squire and Bhattacharjeeb) becomes important. On the other hand, as we will show below, the eddy turnover time in RRMHD turbulence becomes shorter than the disc rotation time, i.e. 
$k_{\perp} u_{\perp}/\varOmega \gg 1$ immediately below the injection scale (see figure 7). Therefore, we do not need to consider a long-time evolution with 
$\varOmega t \sim \epsilon ^{-1}$.
As we shall see in the next section, when 
$\varOmega \ne 0$, this system can be MRI unstable. In the turbulent state, the magnitudes of the nonlinear terms in (2.3a)–(2.3d) increase as the cascade proceeds to smaller scales, and at some point, the linear terms that are proportional to 
$\varOmega$ become negligible. Below the scale at which this happens, the turbulence is governed by standard RMHD, and thus Alfvénic and compressive fluctuations are decoupled. As we will see below, this critical scale roughly corresponds to the scale at which the eddy turnover time becomes shorter than 
$\varOmega ^{-1}$. In other words, when an eddy's lifetime is much shorter than the orbital time of the disc, the effects of the disc's rotation are insignificant. Therefore, the transient growth effects (Balbus & Hawley Reference Balbus and Hawley1992a; Mamatsashvili et al. Reference Mamatsashvili, Chagelishvili, Bodo and Rossi2013) are absent. We also note that, with the normalizations 
$t \varOmega \to t$, 
$z \varOmega /v_A \to z$, 
$x/L_{\perp} \to x$, 
$\varPhi /L_{\perp}^2\varOmega \to \varPhi$, 
$\varPsi /L_{\perp}^2\varOmega \to \varPsi$, 
$u_{\|}/L_{\perp}\varOmega \to u_{\|}$ and 
$v_A \delta B_{\|}/B_0 L_{\perp}\varOmega \to \delta B_{\|}$, the rotation rate is no longer a free parameter, and the only remaining parameter is 
$c_S^2/v_A^2 = \varGamma \beta /2$, where 
$\beta = 8{\rm \pi} p_0/B_0^2$.
The nonlinear free-energy invariant of (2.3a)–(2.3d) consists of Alfvénic and compressive portions 
$W_\mathrm {tot} = W_\mathrm {AW} + W_\mathrm {compr}$, where
$$\begin{gather} W_\mathrm{AW} = \frac{1}{2}\int\mathrm{d}^3{\boldsymbol{r}} \left[ |\boldsymbol{\nabla}_{\perp}\varPhi|^2 + |\boldsymbol{\nabla}_{\perp}\varPsi|^2 \right] , \end{gather}$$
$$\begin{gather}W_\mathrm{compr} = \frac{1}{2}\int\mathrm{d}^3{\boldsymbol{r}} \left[ u_{\|}^2 + v_A^2\left( 1 + \frac{v_A^2}{c_S^2} \right)\frac{\delta B_{\|}^2}{B_0^2} \right]. \end{gather}$$
The time evolution of 
$W_\mathrm {AW}$ and 
$W_\mathrm {compr}$ is given by
$$\begin{gather} \frac{{\rm dW}_\mathrm{AW}}{{\rm d}t} ={-} 2\varOmega\int\mathrm{d}^3{\boldsymbol{r}} u_{\|}\frac{{\partial} \varPhi}{{\partial} y} \equiv I_\mathrm{AW}, \end{gather}$$
$$\begin{gather}\frac{{\rm dW}_\mathrm{compr}}{{\rm d}t} = q\varOmega\int\mathrm{d}^3{\boldsymbol{r}} \left[ v_A\frac{\delta B_{\|}}{B_0}\frac{{\partial} \varPsi}{{\partial} y} - u_{\|}\frac{{\partial} \varPhi}{{\partial} y} \right] + 2\varOmega\int\mathrm{d}^3{\boldsymbol{r}} u_{\|}\frac{{\partial} \varPhi}{{\partial} y} \equiv I_\mathrm{compr}, \end{gather}$$
where 
$I_\mathrm {AW}$ and 
$I_\mathrm {compr}$ are the energy-injection rates of Alfvénic and compressive fluctuations. Noticing that the second term of 
$I_\mathrm {compr}$ is identical to 
$-I_\mathrm {AW}$, we may write 
$I_\mathrm {compr} = I_\mathrm {MRI} - I_\mathrm {AW}$. Then the net injection rate by the MRI is 
$I_\mathrm {MRI}$, and it goes into compressive fluctuations, which then exchange energy with Alfvénic fluctuations at the rate 
$I_\mathrm {AW}$ via linear coupling.
3 Linear MRI of RRMHD
Next, we compare the linear MRI growth rate of full MHD and RRMHD to show that RRMHD can capture the MRI growth rate of full MHD when 
${\boldsymbol {B}}_0$ is nearly azimuthal, viz., 
$\sin \theta \ll 1$. Here, we focus on the modes that are symmetric with respect to the rotation axis 
$\hat {{\boldsymbol {Z}}}$, viz., 
$k_Y = 0$, which is equivalent to 
$k_y = -k_z/\tan \theta$. We focus on these modes because they are the fastest-growing ones. The linear dispersion relation of full MHD (Balbus & Hawley Reference Balbus and Hawley1998, (99)) is
\begin{align} & \left[\omega^2 -
(k_{\|}v_A)^2\right]\left\{ \omega^4 - \left[k_x^2 +
(k_{\|}/\sin\theta)^2\right](c_S^2 + v_A^2) \omega^2 +
(k_{\|} v_A)^2\left[k_x^2 +
(k_{\|}/\sin\theta)^2\right]c_S^2 \right\} \nonumber\\
&\quad = 2\left\{ (2 - q)\omega^4 - k_{\|}^2\left[ 2c_S^2 +
(2\cos^2\theta - q) (c_S^2 + v_A^2)/\sin^2\theta
\right]\omega^2 - qc_S^2 v_A^2 k_{\|}^4/\sin^2\theta
\right\}\varOmega^2. \end{align}
On the other hand, the dispersion relation of RRMHD, (2.3a)-(2.3d), is
\begin{align} & \left[ \omega^2 - (k_{\|} v_A)^2 \right]\left[ \left(1 + \frac{v_A^2}{c_S^2}\right) \omega^2 - (k_{\|} v_A)^2 \right] \nonumber\\ &\quad = 2\left[ (2 - q)\left(1 + \frac{v_A^2}{c_S^2}\right) \omega^2 + q(k_{\|} v_A)^2 \right]\frac{k_y^2}{k_{\perp}^2} \varOmega^2. \end{align}
One can show that, for both (3.1) and (3.2), 
$k_x = 0$ gives the fastest-growing mode. For the RRMHD dispersion relation (3.2), the growth rate does not depend on 
$k_y$ when 
$k_x = 0$. When 
$\varOmega = 0$, (3.1) recovers the Alfvén, slow, and fast modes, while (3.2) recovers the Alfvén and slow modes (the fast mode is eliminated in the RMHD ordering). The maximum growth rate of RRMHD is given by
\begin{equation} \frac{\gamma_\mathrm{max}}{\varOmega} = \sqrt{\frac{5}{18}\beta[20\beta + 15 - \sqrt{8(50\beta^2 + 75\beta + 18})]}, \end{equation}
where we have used 
$q = 3/2$ and 
$\varGamma = 5/3$. One finds that 
$\gamma _\mathrm {max}$ is an increasing function of 
$\beta$ ranging from 
$\gamma _\mathrm {max}/\varOmega \to 0$ for 
$\beta \to 0$Footnote 4 to 
$\gamma _\mathrm {max}/\varOmega \to 3/4$ for 
$\beta \to \infty$. Note that the high-
$\beta$ limit of the maximum growth rate in RRMHD is the same as in full MHD (Balbus & Hawley Reference Balbus and Hawley1998, (114)), and the stabilization of MRI at 
$\beta \to 0$ is consistent with the study by Kim & Ostriker (Reference Kim and Ostriker2000), who found that MRI in full MHD was stabilized when 
$\beta \to 0$ and 
$\theta < 30^\circ$.
In figure 2, we compare the solutions with (3.1) and (3.2). Figure 2(a) shows the growth rate obtained with RRMHD for different values of 
$\beta$; one finds that 
$\gamma _\mathrm {max}$ of RRMHD decreases as 
$\beta$ decreases as expected from (3.3). Figure 2(b–d) shows the growth rate obtained with full MHD for different values of 
$\beta$ and 
$\theta$. For full MHD, the growth rate does not depend on 
$\beta$ when 
$\theta = {\rm \pi}/2$; however, when 
$\sin \theta \ll 1$, the growth rate decreases as 
$\beta$ decreases. Clearly, the growth rates in RRMHD match those in full MHD with 
$\sin \theta \ll 1$, meaning that RRMHD captures the fastest-growing MRI modes when 
${\boldsymbol {B}}_0$ is nearly azimuthal.
Figure 2. The linear MRI growth rate of (a) RRMHD and (b–d) full MHD. The line colours correspond to the values of 
$\beta$ as given in the legend of (a). The line thickness for (b–d) corresponds to the value of 
$\theta$ as given in the legends of these panels. The horizontal dotted lines indicate that, independently of 
$\beta$, the maximum growth rates in RRMHD coincide with those in full MHD in the limit of 
$\theta \to 0$.
4 Simulation of MRI turbulence in RRMHD
Next, we carry out nonlinear simulations of the RRMHD equations to compute the energy partition between the Alfvénic and compressive fluctuations in the saturated state of MRI turbulence. We solve (2.3a)–(2.3d) using a three-dimensional pseudo-spectral code Calliope (Kawazura Reference Kawazura2022). In order to terminate the energy cascade at small scales, we add hyper-viscous and hyper-resistive terms proportional to 
$k_\perp ^8$ and 
$k_z^8$ to the right-hand sides of (2.3a)–(2.3d). As mentioned above, the Alfvénic and compressive fluctuations are expected to be decoupled below some critical scale where the nonlinear terms start to dominate the linear terms. We set the computational grids so that this critical scale is well resolved, which we confirm later in this section. Therefore, the dissipation caused by the hyper-viscosity and hyper-resistivity in (2.3a) and (2.3b) is a measure of the energy flux carried by the Alfvénic fluctuations. Likewise, we can measure the energy flux carried by compressive fluctuations via the hyper-dissipation in (2.3c) and (2.3d). We denote the dissipation rates of the Alfvénic and compressive fluctuations by 
$D_\mathrm {AW}$ and 
$D_\mathrm {compr}$, respectively. The power balance of the system is then
\begin{equation} \frac{{\rm d}W_\mathrm{tot}}{{\rm d}t} = I_\mathrm{AW} + I_\mathrm{compr} - D_\mathrm{AW} - D_\mathrm{compr}. \end{equation}
In a statistically stationary state, 
$\langle I_\mathrm {AW} \rangle + \langle I_\mathrm {compr} \rangle = \langle D_\mathrm {AW} \rangle + \langle D_\mathrm {compr}\rangle$, where 
$\langle \cdots \rangle$ denotes the time average. We set the box size of the simulations as 
$(L_x, L_y, L_z) = (8{\rm \pi} L_{\perp}, 2{\rm \pi} L_{\perp}, 8{\rm \pi} v_A/\varOmega )$ which is discretized by ‘low-resolution grids’ 
$(n_x, n_y, n_z) = (512, 128, 512)$, ‘medium-resolution grids’ 
$(n_x, n_y, n_z) = (1024, 256, 1024)$ and ‘high-resolution grids’ 
$(n_x, n_y, n_z) = (1024, 256, 2048)$. We choose 
$L_z$ so that the fastest-growing mode (
$k_{z} v_A/\varOmega \simeq 1$, as seen in figure 2) fits in the box. We investigate three cases: 
$\beta = 0.1, 1$ and 
$10$. For all of these values of 
$\beta$, we start the simulation with the low-resolution grids and run for a sufficiently long time in the nonlinearly saturated state until 
$\langle \mathrm {d} W_\mathrm {tot}/\mathrm {d} t \rangle \simeq 0$ is satisfied before restarting with the higher-resolution grids.
Figure 3 shows the time evolution of the free energy (
$W_\mathrm {AW}$ and 
$W_\mathrm {compr}$), the power balance (
$I_\mathrm {AW}$, 
$I_\mathrm {compr}$, 
$D_\mathrm {AW}$, 
$D_\mathrm {compr}$), the compressive-to-Alfvénic energy-injection ratio 
$I_\mathrm {compr}/I_\mathrm {AW}$ and the dissipation ratio 
$D_\mathrm {compr}/D_\mathrm {AW}$. From panel (a), one finds that the linear-growth phase occurs at 
$10 \lesssim \varOmega t \lesssim 20$ and is followed by the nonlinearly saturated turbulent phase. While the Alfvénic energy consists predominantly of 
$\delta B_{\perp}$, the compressive energy has almost the same contribution from 
$u_{\|}$ and 
$\delta B_{\|}$. We have confirmed that this trend is the same for 
$\beta = 0.1$ and 
$10$. Panel (b) shows that the energy injection balances with the energy dissipation. Interestingly, the amount of Alfvénic injection 
$I_\mathrm {AW}$ balances with Alfvénic dissipation 
$D_\mathrm {AW}$, and likewise, the compressive injection 
$I_\mathrm {compr}$ and dissipation 
$D_\mathrm {compr}$ are in balance. So, in the saturated state, there is, on average, barely any net nonlinear energy exchange between the two components of the turbulence – even at larger scales, where they are not formally decoupled. We have confirmed that this is also the case for 
$\beta = 0.1$ and 
$10$. As we will see later, this may be due to the fact that the critical scale at which the Alfvénic and compressive fluctuations decouple is located close to the injection scale (see figures 6 and 7). Panel (c) shows the evolution of 
$I_\mathrm {compr}/I_\mathrm {AW}$ and 
$D_\mathrm {compr}/D_\mathrm {AW}$. One finds that both ratios are 
${\simeq } 2\unicode{x2013}2.5$ in the nonlinear state. These values are almost the same for the runs with low-resolution grids (solid lines), medium-resolution grids (dashed lines) and high-resolution grids (dash-dotted lines).
Figure 3. Time evolution of the 
$\beta = 1$ run: (a) each component of the free energy (2.4a) and (2.4b) normalized by the total energy averaged over the nonlinearly saturated state, i.e. over the time interval 
$285\le \varOmega t \le 330$; (b) injection and dissipation rates of Alfvénic and compressive fluctuations normalized by the total injection power averaged over the nonlinearly saturated state; (c) the compressive-to-Alfvénic ratio of injection power 
$I_\mathrm {compr}/I_\mathrm {AW}$ and dissipation rate 
$D_\mathrm {compr}/D_\mathrm {AW}$. The solid, dashed, and dash-dotted lines correspond to the runs with the low-, medium- and high-resolution grids, respectively. The shaded region indicates the interval used for the time averaging.
Figure 4 shows snapshots of the turbulent fields. Structures are elongated in the 
$x$ direction, corresponding to the remnants of ‘channel flows’ driven by MRI, also seen in other shearing-box simulations of full MHD (e.g. Hawley & Balbus Reference Hawley and Balbus1992; Hirai et al. Reference Hirai, Katoh, Terada and Kawai2018). Note that our 
$\hat {{\boldsymbol {y}}}$ direction is almost vertical within the accretion disc, 
$\hat {{\boldsymbol {y}}} \simeq -{\boldsymbol {\varOmega }}/\varOmega$ (see figure 1). For the Alfvénic fields, one finds that 
${\boldsymbol {u}}_{\perp}$ has smaller-scale filamentary structures than 
$\delta {\boldsymbol {B}}_{\perp}$. In contrast, for the compressive fields, the level of filamentation is the same between 
$u_{\|}$ and 
$\delta B_{\|}$. We have found this tendency also for the 
$\beta = 0.1$ and 10 cases.
Figure 4. Snapshots of (a–d) 
$|{\boldsymbol {u}}_{\perp}|$, 
$|\delta {\boldsymbol {B}}_{\perp}|$, 
$|u_{\|}|$ and 
$|\delta B_{\|}|$, each normalized by its own root-mean-square value. These snapshots are taken at 
$\varOmega t = 395$, 
$z = 0$ and for 
$\beta = 1$.
The difference of filamentation levels is more transparent in figure 5, which shows the energy spectra of all fields vs 
$k_{\perp}$, compensated by 
$k_{\perp}^{3/2}$. Here, the energy spectrum of each integrand in (2.4a) and (2.4b) is denoted by 
$E$ with the corresponding subscript. We find that, for the compressive fields, both 
$E_{u_{\|}}$ and 
$E_{\delta B_{\|}}$ have 
$\simeq -3/2$ slope, while the slopes of the Alfvénic fields are not identifiable with the current numerical resolutionFootnote 5 . Independently of 
$\beta$, 
$E_{u_{\perp}}$ is subdominant compared with 
$E_{\delta B_{\perp}}$ at the injection scales, whereas 
$E_{u_{\|}}$ and 
$E_{\delta B_{\|}}$ have almost the same amplitudes throughout the whole 
$k_{\perp}$ range. It is well known that full-MHD simulations of MRI turbulence yield magnetically dominated spectra at large scales (e.g. Lesur & Longaretti Reference Lesur and Longaretti2011; Kimura et al. Reference Kimura, Toma, Suzuki and Inutsuka2016; Walker, Lesur & Boldyrev Reference Walker, Lesur and Boldyrev2016; Sun & Bai Reference Sun and Bai2021), due to generation of azimuthal magnetic field through the shear-flow effect. However, this mechanism cannot explain 
$E_{\delta B_{\perp}} \gg E_{u_{\perp}}$ in RRMHD because the shear flow does not directly produce 
$\delta {\boldsymbol {B}}_{\perp}$, as one can see in (2.3a). Instead we can explain the dominance of 
$E_{\delta B_{\perp}}$ by the linear relation 
$\varPsi /\varPhi = k_{\|} v_A/\gamma$ given by (2.3a), where 
$\gamma$ is the growth rate of MRI. For the fastest-growing mode, 
$k_{\|} v_A/\varOmega \simeq 1$ and 
$\gamma /\varOmega < 1$ (see figure 2), meaning that the linear MRI in RRMHD excites 
$\delta B_{\perp}$ preferentially over 
$u_{\perp}$. One also finds from figure 5 that the disparity between 
$\delta B_{\perp}$ and 
$u_{\perp}$ gets smaller as 
$\beta$ increases. More specifically, at 
$k_{\perp} L_{\perp} = 1$, 
$E_{\delta B_{\perp}}/E_{u_{\perp}} \simeq 27$, 12 and 10 for 
$\beta = 0.1$, 1 and 10, respectively, being consistent with the fact that 
$\gamma _\mathrm {max}/\varOmega$ is an increasing function of 
$\beta$. Nonetheless, the absolute values of the ratio are somewhat different from the linear estimate 
$(\varPsi /\varPhi )^2 \simeq 14$, 3 and 2 for 
$\beta = 0.1$, 1 and 10, respectively, for the fastest-growing mode. This indicates that the nonlinear effect is important and, indeed, as we will see in figure 8, the partition of energy flux between Alfvénic and compressive fluctuations is different between the linear calculation and nonlinear simulations.
Figure 5. Magnetic and kinetic spectra compensated by 
$k_{\perp}^{3/2}$ and averaged over the time interval shown by the shaded area in figure 3 for high-resolution runs with (a) 
$\beta = 0.1$, (b) 
$\beta = 1$ and (c) 
$\beta = 10$. The dashed lines indicate the 
$-$3/2 and 
$-$5/3 slopes.
It is worthwhile to compare our spectra with the incompressible MRI simulation by Walker et al. (Reference Walker, Lesur and Boldyrev2016), which is the highest-resolution shearing-box turbulence to date. They found that the slope of the magnetic field spectrum was close to 
$-$3/2 when the azimuthal component 
$B_Y$ was subtracted. They also found nearly 
$-$3/2 spectral slope for the velocity field as well. These spectra bear a resemblance to our spectra (figure 5). Note, however, that 
$B_Y$ is not necessarily the true mean magnetic field 
${\boldsymbol {B}}_0$, and thus, their magnetic spectrum 
$B_X^2 + B_Z^2$ is presumably a mixture of parallel and perpendicular fluctuations.
In order to investigate the decoupling of Alfvénic and compressive fluctuations, we compare the spectra of energy injection via MRI (
$I_\mathrm {MRI}$), the energy exchange between the Alfvénic and compressive fluctuations (
$I_\mathrm {AW}$) and the nonlinear energy transfer, defined below. Since the coupling between the Alfvénic and compressive fluctuations exists only through the linear terms, the two types of fluctuations are decoupled when the nonlinear energy transfer overwhelms 
$I_\mathrm {AW}$. The nonlinear energy transfer from all modes with wavenumber magnitudes smaller than 
$k_{\perp}$ are defined by (Alexakis, Mininni & Pouquet Reference Alexakis, Mininni and Pouquet2005; Grete et al. Reference Grete, O'Shea, Beckwith, Schmidt and Christlieb2017; St-Onge et al. Reference St-Onge, Kunz, Squire and Schekochihin2020)
\begin{align} \mathcal{N}_\mathrm{AW}^{{<}k_{\perp}} &\equiv \sum_{|q_{\perp}| = k_{\perp} } [-{\boldsymbol{u}}_{\perp{\boldsymbol{q}}}^*\boldsymbol{\cdot}({\boldsymbol{u}}_{\perp}\boldsymbol{\cdot}\boldsymbol{\nabla}_{\perp}{\boldsymbol{u}}_{\perp}^{[{<}k_{\perp}]})_{\boldsymbol{q}} + {\boldsymbol{u}}_{\perp{\boldsymbol{q}}}^*\boldsymbol{\cdot}({\boldsymbol{b}}_{\perp}\boldsymbol{\cdot}\boldsymbol{\nabla}_{\perp}{\boldsymbol{b}}_{\perp}^{[{<}k_{\perp}]})_{\boldsymbol{q}}\end{align}
\begin{align} &\quad - {\boldsymbol{b}}_{\perp{\boldsymbol{q}}}^*\boldsymbol{\cdot}({\boldsymbol{u}}_{\perp}\boldsymbol{\cdot}\boldsymbol{\nabla}_{\perp}{\boldsymbol{b}}_{\perp}^{[{<}k_{\perp}]})_{\boldsymbol{q}} + {\boldsymbol{b}}_{\perp{\boldsymbol{q}}}^*\boldsymbol{\cdot}({\boldsymbol{b}}_{\perp}\boldsymbol{\cdot}\boldsymbol{\nabla}_{\perp}{\boldsymbol{u}}_{\perp}^{[{<}k_{\perp}]})_{\boldsymbol{q}}] \end{align}
\begin{align} \mathcal{N}_\mathrm{compr}^{{<}k_{\perp}} &\equiv \sum_{|q_{\perp}| = k_{\perp} }[{-}u_{\|{\boldsymbol{q}}}^*\boldsymbol{\cdot} ({\boldsymbol{u}}_{\perp}\boldsymbol{\cdot}\boldsymbol{\nabla}_{\perp} u_{\|}^{[{<}k_{\perp}]})_{\boldsymbol{q}} + u_{\|{\boldsymbol{q}}}^*\boldsymbol{\cdot} ({\boldsymbol{b}}_{\perp}\boldsymbol{\cdot}\boldsymbol{\nabla}_{\perp} b_{\|}^{[{<}k_{\perp}]})_{\boldsymbol{q}}\end{align}
\begin{align} &\quad - b_{\|{\boldsymbol{q}}}^*\boldsymbol{\cdot}({\boldsymbol{u}}_{\perp}\boldsymbol{\cdot}\boldsymbol{\nabla}_{\perp} b_{\|}^{[{<}k_{\perp}]})_{\boldsymbol{q}} + b_{\|{\boldsymbol{q}}}^*\boldsymbol{\cdot}({\boldsymbol{b}}_{\perp}\boldsymbol{\cdot}\boldsymbol{\nabla}_{\perp} u_{\|}^{[{<}k_{\perp}]})_{\boldsymbol{q}}], \end{align}
where 
$f^{[< k_{\perp}]}({\boldsymbol {r}}) \equiv \sum _{q_z}\sum _{q_{\perp} < k_{\perp}} f_{\boldsymbol {q}}\,\mathrm {e}^{\mathrm {i}{\boldsymbol {q}}\boldsymbol {\cdot }{\boldsymbol {r}}}$, 
${\boldsymbol {b}} \equiv v_A \delta {\boldsymbol {B}}/B_0$. We also define the spectra of MRI injection 
$I_\mathrm {MRI} = \sum _{\boldsymbol {k}}\mathcal {I}_\mathrm {MRI}({\boldsymbol {k}})$ and energy exchange 
$I_\mathrm {AW} = \sum _{\boldsymbol {k}}\mathcal {I}_\mathrm {AW}({\boldsymbol {k}})$. The top panels of figure 6 show the perpendicular spectra of injection, exchange and nonlinear energy transfer. Both the injection and exchange peak near the box scale and drop quickly at smaller scales, while the nonlinear energy transfer is relatively constant throughout the 
$k_{\perp}$-range. Consequently, the nonlinear energy transfer overwhelms the injection and coupling immediately below the box scale. The bottom panels of figure 6 show the 
$z$ spectra of the same quantities. The peak of the injection is located around 
$k_z v_A/\varOmega \simeq 1$ and, thus, the injection scale corresponds to the fastest-growing modes. In the same way as the perpendicular spectra, the injection and exchange drop quickly at scales smaller than that of the fastest-growing mode and are overwhelmed by the nonlinear energy transfer. Therefore, in the small scales of our simulations, the coupling between Alfvénic and compressive fluctuations is negligible.
Figure 6. The spectra of energy injection via MRI, energy exchange between Alfvénic and compressive fluctuations, dissipation of Alfvénic and compressive fluctuations and nonlinear transfer vs (a–c) 
$k_\perp$ and (d–f) 
$k_z$ for (a,d) 
$\beta = 0.1$, (b,e) 
$\beta = 1$ and (c, f) 
$\beta = 10$. The spectra are normalized by 
$\langle I_\mathrm {MRI} \rangle$, integrated over (a–c) 
$k_z$ and (d–f) 
$k_{\perp}$, and averaged over the time interval shown by the shaded area in figure 3.
While the spectral comparison shown in figure 6 is the most direct proof of the decoupling of Alfvénic and compressive fluctuations, we expect that the ratio between the eddy turnover rate and the angular velocity of the accretion disc can also be a proxy for the measurement of the decouplingFootnote 6 . In figure 7, we plot the eddy turnover rate 
$\omega _\mathrm {nl} \sim |{\boldsymbol {u}}_{\perp}\boldsymbol {\cdot }\boldsymbol {\nabla }_{\perp}| \sim k_{\perp} u_{\perp} \sim k_{\perp}^{3/2}E_{u_{\perp}}^{1/2}$ normalized by 
$\varOmega$. One finds that this value is an increasing function of 
$k_{\perp} L_{\perp}$ and exceeds unity at some scale. As mentioned above, when 
$\omega _\mathrm {nl}/\varOmega$ is much larger than unity, the effect of differential rotation is expected to be negligible, and the turbulence obeys the standard RMHD where the Alfvénic and compressive fluctuations are decoupled (Schekochihin et al. Reference Schekochihin, Cowley, Dorland, Hammett, Howes, Quataert and Tatsuno2009). The scale at which 
$\omega _\mathrm {nl}/\varOmega \simeq 1$ is not much smaller than the injection scale, which is consistent with the fact that the nonlinear energy transfer overwhelms the MRI injection immediately below the injection scale, as shown in figure 6. In general, 
$\omega _\mathrm {nl}/\varOmega$ is easier to use as an indicator of decoupling because computing the nonlinear energy transfer is numerically cumbersome.
Figure 7. Normalized eddy turnover frequency 
$\omega _\mathrm {nl}/\varOmega$ vs 
$k_{\perp} L_{\perp}$ averaged over the time interval shown by the shaded area in figure 3 for higher-resolution runs. The effect of differential rotation is negligible where the value is greater than unity.
As the decoupling of Alfvénic and compressive fluctuations in our simulations has been demonstrated in figures 6 and 7, we now calculate the partition of energy flux carried by these two types of fluctuations. Figure 8 shows the dependence of 
$\langle D_\mathrm {compr} \rangle /\langle D_\mathrm {AW} \rangle$ on 
$\beta$. We find that 
$\langle D_\mathrm {compr} \rangle /\langle D_\mathrm {AW} \rangle$ is between 2 and 2.5 for all values of 
$\beta$ that we studied, without an obvious trend. Since Alfvénic and compressive fluctuations are decoupled, the result in figure 8 would not be changed for a finer-resolution simulation. Indeed, we have found almost identical values of 
$\langle D_\mathrm {compr} \rangle /\langle D_\mathrm {AW} \rangle$ in our simulations conducted at all resolutions, from low to high. This is because, as seen in figure 6, even the low-resolution grid runs resolve the critical scale where the nonlinear energy transfer dominates the linear coupling.
Figure 8. Partition of energy flux between Alfvénic and compressive fluctuations vs 
$\beta$: (black) 
$I_\mathrm {compr}/I_\mathrm {AW}$ calculated by the eigenfunctions of linear dispersion relation (3.2) and 
$D_\mathrm {compr}/D_\mathrm {AW}$ calculated by the nonlinear simulation with (blue) low-resolution grids, (orange) medium-resolution grids, (green) high-resolution grids. Error bars for the nonlinear simulations are estimated by calculating the standard deviation over the averaging interval.
Note that the values of 
$\langle D_\mathrm {compr} \rangle /\langle D_\mathrm {AW} \rangle$ obtained from our nonlinear simulations are different from the values of 
$I_\mathrm {compr}/I_\mathrm {AW}$ (see (B3) for the definition) computed ‘quasilinearly’ for the fastest-growing linear MRI modes (black dashed line in figure 8), the latter value being close to unity. This indicates that, even though the decoupling of Alfvénic and compressive fluctuations starts relatively near the injection scale, the preferential excitation of compressive fluctuations in MRI turbulence is the consequence of nonlinear effects, i.e. of the way in which the nonlinearity removes the energy injected by the MRI from the injection scale and transfers it into the two turbulent cascade.
5 Application to ion-to-electron heating prescription in hot accretion flows
In this section, we discuss the application of our findings to hot accretion flows, such as M87 and Sgr A*, together with some important caveats. Numerical simulations of gyrokinetic turbulence have revealed that the partition between ion and electron heating is crucially sensitive to the compressive-to-Alfvénic injection power ratio 
$P_\mathrm {compr}/P_\mathrm {AW}$ at the ion Larmor scale (Kawazura et al. Reference Kawazura, Schekochihin, Barnes, TenBarge, Tong, Klein and Dorland2020)Footnote 7 . Since the compressive and Alfvénic energy fluxes computed in our simulations are supposed to cascade down to the ion Larmor scale independently, it is straightforward to infer that 
$P_\mathrm {compr}/P_\mathrm {AW}$ is equal to 
$D_\mathrm {compr}/D_\mathrm {AW} \simeq 2\unicode{x2013}2.5$, as we found in figure 8. Therefore, we can combine the results of this paper with our previous study of gyrokinetic turbulence to formulate the ion-to-electron heating prescription that incorporates both driving of turbulence via MRI at MHD scales and the dissipation at kinetic scales. Substituting 
$P_\mathrm {compr}/P_\mathrm {AW} = 2$ in (14) of Kawazura et al. (Reference Kawazura, Schekochihin, Barnes, TenBarge, Tong, Klein and Dorland2020), one obtains
\begin{equation} \frac{Q_\mathrm{i}}{Q_e}(\beta_\mathrm{i},\ T_\mathrm{i}/T_e) = \frac{35}{1 + (\beta_\mathrm{i}/15)^{{-}1.4} \exp({-0.1/(T_\mathrm{i}/T_e)})} + 2, \end{equation}
where 
$T_\mathrm {i}/T_e$ is the ion-to-electron background temperature ratio, and 
$\beta _\mathrm {i}$ is the ion beta.
This prescription is a step forward from the currently used heating prescription and may help improve the quality of hot accretion flow modelling (e.g. Chael et al. Reference Chael, Rowan, Narayan, Johnson and Sironi2018; Chael, Narayan & Johnson Reference Chael, Narayan and Johnson2019). However, one must bear in mind that a number of heating channels are missing in (5.1). First, we do not consider spiral density waves (Heinemann & Papaloizou Reference Heinemann and Papaloizou2009) which are outside the RMHD ordering as they have no vertical structure, i.e. 
$k_Z \simeq 0$. The excitation of these waves may change the partition between Alfvénic and compressive fluctuations. Note that these waves form weak shocks and dissipate into thermal energy, but the amount of heating due to this dissipation is very little.
Second, while we have only considered collisional MRI in this study, the mean free path of hot accretion flows is almost equal to, or longer than, the scale height of the disc, meaning that MRI is supposed to be collisionlessFootnote 8 . When MRI is collisionless, the viscous stress due to pressure anisotropy gives rise to a new heating channel. Approximately 50 % of total injected power may be directly converted into heat at large scales by this viscous stress, which would not cascade down to the ion Larmor scale (Sharma et al. Reference Sharma, Quataert, Hammett and Stone2007; Kempski et al. Reference Kempski, Quataert, Squire and Kunz2019)Footnote 9 .
Third, even if these additional heating channels at large scales are absent, there are other heating channels at the ion Larmor scale that are not captured by standard gyrokinetics (see § II A in Kawazura et al. (Reference Kawazura, Schekochihin, Barnes, TenBarge, Tong, Klein and Dorland2020), for a detailed discussion), e.g. cyclotron heating (Cranmer, Field & Kohl Reference Cranmer, Field and Kohl1999), stochastic heating (Chandran et al. Reference Chandran, Li, Rogers, Quataert and Germaschewski2010) and background pressure anisotropy (Kunz et al. Reference Kunz, Abel, Klein and Schekochihin2018).
Thus, our heating prescription (5.1) is only the simplest possible model that considers both MRI injection and kinetic dissipation. Including the missing heating channels is an important task for future work.
6 Conclusions
In this study, we have calculated the energy partition between Alfvénic and compressive fluctuations in turbulence driven by MRI with near-azimuthal mean magnetic fields. The fastest-growing MRI modes are correctly captured by RMHD with differential rotation (RRMHD) because they satisfy 
$|k_Z/k_Y| \sim |k_{\perp}/k_{\|}| \gg 1$ when the background field is nearly azimuthal (Balbus & Hawley Reference Balbus and Hawley1992a). In RRMHD, the Alfvénic and compressive fluctuations are coupled only through the linear terms that are proportional to the angular velocity of the accretion disc. We have carried out nonlinear simulations of RRMHD and showed that the nonlinear energy transfer overwhelms the linear coupling immediately below the injection scale. Thus, the two kinds of fluctuations are decoupled at the small scales in our simulations. This is because, below the injection scale, the eddy turnover time is much shorter than the disc rotation time, i.e. 
$\omega _\mathrm {nl}/\varOmega \gg 1$. Most importantly, the energy flux carried by the compressive fluctuations is more than double that carried by the Alfvénic fluctuations at the decoupled scales – a result reflecting the interaction between MRI injection and nonlinearity at the injection scale and distinct from a ‘quasilinear’ estimate (which suggests near equipartition).
While these findings suggest that RRMHD is a useful model for studying MRI turbulence, we would like to stress the following two limitations of the RMHD approach for MRI-driven turbulence in accretion flows. First, we assume a near-azimuthal constant mean magnetic field. This may be quite restrictive: e.g. global MHD simulations (e.g. Suzuki & Inutsuka Reference Suzuki and Inutsuka2014) sometimes exhibit non-azimuthal components of magnetic field. Secondly, we assume that 
$k_{\|}/k_{\perp} \ll 1$ is already satisfied at a larger scale than the critical scale where 
$\omega _\mathrm {nl}/\varOmega \sim 1$. If this were not to hold, the rotation effects in full MHD may become negligible at scales larger than those where the RMHD approximation is already satisfied, and our RRMHD model would not be a good model of MRI turbulence at the decoupling scale. In such a case, the turbulence in the RMHD range would not be driven by MRI, but by the cascade from the full-MHD scales. A simulation of full MHD with extreme resolutions is necessary to explore this possibility.
Acknowledgements
Y.K. thanks M. Kunz for fruitful discussions. Y.K., A.A.S., M.B. and S.A.B. were supported in part by the STFC grant ST/N000919/1; the work of A.A.S. and M.B. was also supported in part by the EPSRC grant EP/R034737/1. Y.K. was supported by JSPS KAKENHI grants JP19K23451 and JP20K14509. Numerical computations reported here were carried out on Cray XC50 at Centre for Computational Astrophysics in National Astronomical Observatory of Japan, on the computing resource at Kyushu University, on Oakforest-PACS and Oakbridge-CX at the University of Tokyo, and on Flow at Nagoya University.
Editor Steve Tobias thanks the referees for their advice in evaluating this article.
Declaration of interests
The author report no conflict of interest.
Appendix A. Derivation of RRMHD model
Here, we explicitly derive (2.3a)–(2.3d) from (2.1a)–(2.1d). The way we do it is mostly the same as the derivation of (17), (18), (25) and (26) in Schekochihin et al. (Reference Schekochihin, Cowley, Dorland, Hammett, Howes, Quataert and Tatsuno2009), but with account taken of the differential rotation of the disc. We start by considering the following ordering:
\begin{equation} \frac{{\boldsymbol{u}}}{v_A} \sim \frac{\delta{\boldsymbol{B}}}{B_0} \sim \frac{k_{\|}}{k_{\perp}} \sim \sin\theta \sim \epsilon,\quad \frac{{\partial} }{{\partial} t} \sim \varOmega \sim k_{\|} v_A \equiv \omega. \end{equation}Then, the order of each term in (2.1a) is estimated as follows:
\begin{align} & \underbrace{\frac{{\partial} }{{\partial} t}\frac{\delta \rho}{\rho_0}}_{\epsilon^1\omega} + \underbrace{u_{\|}\frac{{\partial} }{{\partial} z}\frac{\delta \rho}{\rho_0}}_{\epsilon^2\omega} + \underbrace{{\boldsymbol{u}}_{\perp}\boldsymbol{\cdot}\boldsymbol{\nabla}_{\perp}\frac{\delta \rho}{\rho_0}}_{\epsilon^1\omega} - \underbrace{q\varOmega x\sin\theta\left(\frac{{\partial} }{{\partial} y} + \frac{1}{\tan\theta}\frac{{\partial} }{{\partial} z}\right)}_{\epsilon^2\omega} \frac{\delta \rho}{\rho_0} \nonumber\\ &\quad ={-}\left(\underbrace{\frac{{\partial} u_{\|}}{{\partial} z}}_{\epsilon^1\omega} + \underbrace{\boldsymbol{\nabla}_{\perp}\boldsymbol{\cdot}{\boldsymbol{u}}_{\perp}}_{\epsilon^0\omega}\right) - \frac{\delta \rho}{\rho_0}\left(\underbrace{\frac{{\partial} u_{\|}}{{\partial} z}}_{\epsilon^2\omega} + \underbrace{\boldsymbol{\nabla}_{\perp}\boldsymbol{\cdot}{\boldsymbol{u}}_{\perp}}_{\epsilon^1\omega}\right). \end{align}
To order 
$\mathcal {O}(\epsilon ^0\omega )$, we obtain 
$\boldsymbol {\nabla }_{\perp}\boldsymbol {\cdot }{\boldsymbol {u}}_{\perp} = 0$. Likewise, to lowest order, 
$\boldsymbol {\nabla }\boldsymbol {\cdot } \delta {\boldsymbol {B}} = 0$ gives 
$\boldsymbol {\nabla }_{\perp}\boldsymbol {\cdot }\delta {\boldsymbol {B}}_{\perp} = 0$. Therefore, we may write 
${\boldsymbol {u}}_{\perp}$ and 
$\delta {\boldsymbol {B}}_{\perp}$ in terms of stream and flux functions
\begin{equation} {\boldsymbol{u}}_{\perp}{=} \hat{{\boldsymbol{z}}}\times\boldsymbol{\nabla}_{\perp}\varPhi,\quad \frac{\delta{\boldsymbol{B}}_{\perp}}{B_0} = \frac{\hat{{\boldsymbol{z}}}\times\boldsymbol{\nabla}_{\perp}\varPsi}{v_A}. \end{equation}
Then, the 
$\mathcal {O}(\epsilon ^1\omega )$ terms in (A2) yield
\begin{equation} \left( \frac{{\partial} }{{\partial} t} + {\boldsymbol{u}}_{\perp}\boldsymbol{\cdot}\boldsymbol{\nabla}_{\perp} \right)\frac{\delta \rho}{\rho_0} ={-}\frac{{\partial} u_{\|}}{{\partial} z}. \end{equation}Note that the shearing term, viz., the fourth term on the left-hand side of (A2), is ordered out. As we will show shortly, the shearing terms in other equations are also ordered out.
Under the same ordering, terms in (2.1b) are ordered as follows:
\begin{align} &
\underbrace{\frac{{\partial}
{\boldsymbol{u}}}{{\partial}
t}}_{\epsilon^1\omega v_A} +
\underbrace{u_{\|}\frac{{\partial}
{\boldsymbol{u}}}{{\partial}
z}}_{\epsilon^2\omega v_A} +
\underbrace{{\boldsymbol{u}}_{\perp}\boldsymbol{\cdot}\boldsymbol{\nabla}_{\perp}{\boldsymbol{u}}}_{\epsilon^1\omega
v_A} - \underbrace{q\varOmega
x\sin\theta\left(\frac{{\partial}
}{{\partial} y} +
\frac{1}{\tan\theta}\frac{{\partial}
}{{\partial} z}\right)}_{\epsilon^2\omega
v_A}{\boldsymbol{u}} \nonumber\\ &\quad
={-}\hat{{\boldsymbol{z}}}\frac{{\partial}
}{{\partial}
z}\left[\underbrace{\frac{c_S^2}{\varGamma} \frac{\delta
p}{p_0}}_{\epsilon^1\omega v_A} +
v_A^2\left(\underbrace{\frac{1}{2}\frac{|\delta{\boldsymbol{B}}|^2}{B_0^2}}_{\epsilon^2\omega
v_A} + \underbrace{\frac{\delta
B_{\|}}{B_0}}_{\epsilon^1\omega v_A}\right)\right]
\nonumber\\ &\qquad
-\boldsymbol{\nabla}_{\perp}\left[\underbrace{\frac{c_S^2}{\varGamma}
\frac{\delta p}{p_0}}_{\epsilon^0\omega v_A} +
v_A^2\left(\underbrace{\frac{1}{2}\frac{|\delta{\boldsymbol{B}}|^2}{B_0^2}}_{\epsilon^1\omega
v_A} + \underbrace{\frac{\delta
B_{\|}}{B_0}}_{\epsilon^0\omega v_A}\right)\right]
\nonumber\\ &\qquad
-\frac{\delta\rho}{\rho_0}\hat{{\boldsymbol{z}}}\frac{{\partial}
}{{\partial}
z}\left[\underbrace{\frac{c_S^2}{\varGamma} \frac{\delta
p}{p_0}}_{\epsilon^2\omega v_A} +
v_A^2\left(\underbrace{\frac{1}{2}\frac{|\delta{\boldsymbol{B}}|^2}{B_0^2}}_{\epsilon^3\omega
v_A} + \underbrace{\frac{\delta
B_{\|}}{B_0}}_{\epsilon^2\omega v_A}\right)\right]
\nonumber\\ &\qquad
-\frac{\delta\rho}{\rho_0}\boldsymbol{\nabla}_{\perp}\left[\underbrace{\frac{c_S^2}{\varGamma}
\frac{\delta p}{p_0}}_{\epsilon^1\omega v_A} +
v_A^2\left(\underbrace{\frac{1}{2}\frac{|\delta{\boldsymbol{B}}|^2}{B_0^2}}_{\epsilon^2\omega
v_A} + \underbrace{\frac{\delta
B_{\|}}{B_0}}_{\epsilon^1\omega v_A}\right)\right]
\nonumber\\ &\qquad
+v_A^2\left(\underbrace{\frac{{\partial}
}{{\partial}
z}\frac{\delta{\boldsymbol{B}}}{B_0}}_{\epsilon^1\omega
v_A} + \underbrace{\frac{\delta
B_{\|}}{B_0}\frac{{\partial}
}{{\partial}
z}\frac{\delta{\boldsymbol{B}}}{B_0}}_{\epsilon^2\omega
v_A} +
\underbrace{\frac{\delta{\boldsymbol{B}}_{\perp}}{B_0}\boldsymbol{\cdot}\boldsymbol{\nabla}_{\perp}\frac{\delta{\boldsymbol{B}}}{B_0}}_{\epsilon^1\omega
v_A}\right) \nonumber\\ &\qquad -
2\varOmega(-\underbrace{\cos\theta\hat{{\boldsymbol{y}}}}_{\epsilon^2\omega
v_A} +
\underbrace{\sin\theta\hat{{\boldsymbol{z}}}}_{\epsilon^1\omega
v_A})\times{\boldsymbol{u}} + q\varOmega u_x
(\underbrace{\sin\theta\hat{{\boldsymbol{y}}}}_{\epsilon^2\omega
v_A} +
\underbrace{\cos\theta\hat{{\boldsymbol{z}}}}_{\epsilon^1\omega
v_A}). \end{align}
From the 
$\mathcal {O}(\epsilon ^0\omega v_A)$ terms in (A5), one gets the pressure balance
\begin{equation} \frac{c_S^2}{\varGamma} \frac{\delta p}{p_0} + v_A^2\frac{\delta B_{\|}}{B_0} = 0, \end{equation}which, when combined with (2.1d), becomes
\begin{equation} \frac{\delta \rho}{\rho_0} + \frac{v_A^2}{c_S^2} \frac{\delta B_{\|}}{B_0} = 0. \end{equation}
From the 
$\mathcal {O}(\epsilon ^1\omega v_A)$ terms in (A5), we obtain
\begin{align} \frac{{\partial} {\boldsymbol{u}}}{{\partial} t} + {\boldsymbol{u}}_{\perp}\boldsymbol{\cdot}\boldsymbol{\nabla}_{\perp}{\boldsymbol{u}} &={-} \boldsymbol{\nabla}_{\perp}\left(\frac{v_A^2}{2}\frac{|\delta{\boldsymbol{B}}|^2}{B_0^2}\right) + v_A^2\left(\frac{{\partial} }{{\partial} z}\frac{\delta{\boldsymbol{B}}}{B_0} + \frac{\delta{\boldsymbol{B}}_{\perp}}{B_0}\boldsymbol{\cdot} \boldsymbol{\nabla}_{\perp}\frac{\delta{\boldsymbol{B}}}{B_0}\right) \nonumber\\ &\quad + 2\varOmega\hat{{\boldsymbol{y}}}\times{\boldsymbol{u}} + q\varOmega u_x\hat{{\boldsymbol{z}}}, \end{align}
where we have used 
$\cos \theta \simeq 1$ and neglected all terms containing 
$\sin \theta \ll 1$. The desired perpendicular and parallel momentum equations (2.3b) and (2.3c) are recovered as 
$\hat {{\boldsymbol {z}}}\boldsymbol {\cdot }[\boldsymbol {\nabla }_{\perp}\times$(A8)] and 
$\hat {{\boldsymbol {z}}}\boldsymbol {\cdot }$(A8), respectively.
Next, the ordering of terms in (2.1c) is as follows:
\begin{align} & \underbrace{\frac{{\partial} }{{\partial} t}\frac{\delta {\boldsymbol{B}}}{B_0}}_{\epsilon^1\omega} + \underbrace{u_{\|}\frac{{\partial} }{{\partial} z}\frac{\delta {\boldsymbol{B}}}{B_0}}_{\epsilon^2\omega} + \underbrace{{\boldsymbol{u}}_{\perp}\boldsymbol{\cdot}\boldsymbol{\nabla}_{\perp}\frac{\delta {\boldsymbol{B}}}{B_0}}_{\epsilon^1\omega} - \underbrace{q\varOmega x\sin\theta\left(\frac{{\partial} }{{\partial} y} + \frac{1}{\tan\theta}\frac{{\partial} }{{\partial} z}\right)}_{\epsilon^2\omega} \frac{\delta {\boldsymbol{B}}}{B_0} \nonumber\\ &\qquad + \left(\underbrace{\hat{{\boldsymbol{z}}}}_{\epsilon^1\omega}+ \underbrace{\frac{\delta{\boldsymbol{B}}}{B_0}}_{\epsilon^2\omega}\right)\frac{{\partial} u_{\|}}{{\partial} z} \nonumber\\ &\quad = \underbrace{\frac{{\partial} {\boldsymbol{u}}}{{\partial} z}}_{\epsilon^1\omega} + \underbrace{\frac{\delta B_{\|}}{B_0}\frac{{\partial} {\boldsymbol{u}}}{{\partial} z}}_{\epsilon^2\omega} + \underbrace{\frac{\delta{\boldsymbol{B}}_{\perp}}{B_0}\boldsymbol{\cdot}\boldsymbol{\nabla}_{\perp}{\boldsymbol{u}}}_{\epsilon^1\omega} - q\varOmega\frac{\delta B_x}{B_0}( \underbrace{\sin\theta\hat{{\boldsymbol{y}}}}_{\epsilon^2\omega} + \underbrace{\cos\theta\hat{{\boldsymbol{z}}}}_{\epsilon^1\omega}). \end{align}
Together with (A4) and (A7), the 
$\mathcal {O}(\epsilon ^1\omega )$ terms in this equation yield
\begin{equation} \left( \frac{{\partial} }{{\partial} t} + {\boldsymbol{u}}_{\perp}\boldsymbol{\cdot}\boldsymbol{\nabla}_{\perp} \right)\left( \frac{\delta{\boldsymbol{B}}}{B_0} + \hat{{\boldsymbol{z}}}\frac{v_A^2}{c_S^2}\frac{\delta B_{\|}}{B_0} \right) = \left( \frac{{\partial} }{{\partial} z} + \frac{\delta{\boldsymbol{B}}}{B_0}\boldsymbol{\cdot}\boldsymbol{\nabla}_{\perp} \right){\boldsymbol{u}} - q\varOmega\frac{\delta B_x}{B_0}\hat{{\boldsymbol{z}}}. \end{equation}
Finally, we obtain the perpendicular and parallel magnetic field equations (2.3a) and (2.3d) as 
$\hat {{\boldsymbol {z}}}\boldsymbol {\cdot }[\mathrm {curl}^{-1}$(A10)] and 
$\hat {{\boldsymbol {z}}}\boldsymbol {\cdot }$(A10), respectively.
Appendix B. Compressive-to-Alfvénic energy-injection ratio for a single linear MRI mode in RRMHD
Substituting the solution to the dispersion relation (3.2) back into the linearized RRMHD equations (2.3a)–(2.3d), one gets the linear relations
$$\begin{gather} \frac{\delta B_{\|}}{B_0} = \lambda\frac{k_y\varPsi}{v_A}, \end{gather}$$
$$\begin{gather}\varPhi ={-}\frac{\omega}{k_{\|} v_A}\varPsi, \end{gather}$$
$$\begin{gather}u_{\|} ={-}\frac{\varOmega}{k_{\|} v_A}\left[ \left( 1 + \frac{v_A^2}{c_S^2} \right)\lambda\frac{\omega}{\varOmega} + q \right]k_y\varPsi, \end{gather}$$where
\begin{align} \lambda &= \frac{5\mathrm{i}\beta}{2\sqrt{2}(k_y/k_{\perp})^2(6 + 5\beta)^{3/2}} \left[{-}6\bar{k}_{\|}^2 + 7(k_y/k_{\perp})^2(6 + 5\beta) \vphantom{\sqrt{36\bar{k}_{\|}^4 + (k_y/k_{\perp})^2(6 + 5\beta)^2 + 4(k_y/k_{\perp})^2\bar{k}_{\|}^2(6 + 5\beta)(3 + 20\beta)}}\right.\nonumber\\ &\quad \left.-\sqrt{36\bar{k}_{\|}^4 + (k_y/k_{\perp})^2(6 + 5\beta)^2 + 4(k_y/k_{\perp})^2\bar{k}_{\|}^2(6 + 5\beta)(3 + 20\beta)}\right] \nonumber\\ &\quad \times\left[\bar{k}_{\|}^2(6 + 10\beta) + (k_y/k_{\perp})^2(6 + 5\beta) \vphantom{\sqrt{36\bar{k}_{\|}^4 + (k_y/k_{\perp})^2(6 + 5\beta)^2 + 4(k_y/k_{\perp})^2\bar{k}_{\|}^2 (6 + 5\beta)(3 + 20\beta)}}\right.\nonumber\\ &\quad \left.-\sqrt{36\bar{k}_{\|}^4 + (k_y/k_{\perp})^2(6 + 5\beta)^2 + 4(k_y/k_{\perp})^2\bar{k}_{\|}^2 (6 + 5\beta)(3 + 20\beta)}\right]^{{-}1}, \end{align}
with 
$\bar {k}_{\|} = k_{\|}v_A/\varOmega$. For the fastest-growing mode, 
$\lambda$ reduces to 
$\sqrt {-5\beta /(5\beta + 6)}$. Substituting (B1a)–(B1c) into (2.5a) and (2.5b), one obtains
\begin{equation} \frac{I_\mathrm{compr}}{I_\mathrm{AW}} ={-}\frac{ q\lambda( k_{\|} v_A )^2/[ ( 1 + v_A^2/c_S^2 ) \lambda\varOmega/\omega + q] + (2 - q)\varOmega\omega^*}{2\varOmega\omega^*}, \end{equation}
where the superscript star denotes the complex conjugate. Note that, when the rotation is not sheared, i.e. 
$q = 0$, this becomes the conservation of energy 
$I_\mathrm {compr} + I_\mathrm {AW} = 0$, i.e. the Alfvénic and compressive fluctuations exchange their energy via unsheared rotation.
Open access
