STOCHASTIC PARTICLE ACCELERATION IN TURBULENCE GENERATED BY MAGNETOROTATIONAL INSTABILITY

First, we perform ideal MHD simulations to obtain turbulent fields generated by MRI. We ignore back reaction by CRs for simplicity. We use the isothermal equation of state and do not solve the energy equation. We use the shearing box approximation without vertical gravity for simplicity (Hawley et al. 1995). The shearing box is a local approximation of differentially rotating plasma. We use the co-rotational frame and ignore the curvature effect. Then, the set of MHD equations is written as

$$\begin{eqnarray}&&\displaystyle \frac{\partial \rho }{\partial T}+{\rm{\nabla }}\cdot (\rho {{\boldsymbol{V}}}_{{\rm{fl}}})=0,\end{eqnarray} \tag{ 1 }$$

$$\begin{eqnarray}&&\displaystyle \frac{\partial {\boldsymbol{B}}}{\partial T}=-{\rm{\nabla }}\times ({{\boldsymbol{V}}}_{{\rm{fl}}}\times {\boldsymbol{B}}),\end{eqnarray} \tag{ 2 }$$

$$\begin{eqnarray}&&\displaystyle \frac{\partial {{\boldsymbol{V}}}_{{\rm{fl}}}}{\partial T}+{{\boldsymbol{V}}}_{{\rm{fl}}}\cdot {\rm{\nabla }}{{\boldsymbol{V}}}_{{\rm{fl}}}=-\displaystyle \frac{1}{\rho }\left(P+\displaystyle \frac{{B}^{2}}{8\pi }\right)\\ &&\quad +\displaystyle \frac{({\boldsymbol{B}}\cdot {\rm{\nabla }}){\boldsymbol{B}}}{4\pi \rho }-2{{\boldsymbol{\Omega }}}_{{\rm{K}}}\times {{\boldsymbol{V}}}_{{\rm{fl}}}+3{{\rm{\Omega }}}_{{\rm{K}}}^{2}x{{\boldsymbol{e}}}_{x},\end{eqnarray} \tag{ 3 }$$

where ${{\boldsymbol{e}}}_{i}$ is the unit vector of the x_i direction, ${{\boldsymbol{V}}}_{{\rm{fl}}}$ is the fluid velocity, ${{\boldsymbol{\Omega }}}_{{\rm{K}}}={{\rm{\Omega }}}_{{\rm{K}}}{{\boldsymbol{e}}}_{z}$ is the angular velocity at the box center (${{\rm{\Omega }}}_{{\rm{K}}}$ is the Keplerian angular velocity), T is the time for the MHD simulation, and the other letters have their usual meanings. In the shearing box, the directions $x,\;y,$ and z correspond to $r,\;\phi$, and z in the cylindrical coordinate, respectively. In this paper, we ignore the effect of resistivity, which may affect particle acceleration through the generation of electric fields (see Section 4.3).

The boundary conditions of the shearing box for the magnetic fields are described as

$$\begin{eqnarray}&&{\boldsymbol{B}}(x,y,z)={\boldsymbol{B}}(x+{L}_{x},\;y-1.5{{\rm{\Omega }}}_{{\rm{K}}}{L}_{x}T,\;z),\end{eqnarray} \tag{ 4 }$$

$$\begin{eqnarray}&&{\boldsymbol{B}}(x,y,z)={\boldsymbol{B}}(x,\;y+{L}_{y},\;z),\end{eqnarray} \tag{ 5 }$$

$$\begin{eqnarray}&&{\boldsymbol{B}}(x,y,z)={\boldsymbol{B}}(x,\;y,\;z+{L}_{z}).\end{eqnarray} \tag{ 6 }$$

The same conditions are used for the velocity fields and density. The background shear velocity in a shearing box is ${{\boldsymbol{V}}}_{{\rm{shear}}}=-1.5{{\rm{\Omega }}}_{{\rm{K}}}x{{\boldsymbol{e}}}_{y}$, where the factor of 1.5 comes from the index of the radial dependence of angular velocity, ${{\rm{\Omega }}}_{{\rm{K}}}\propto {r}^{-1.5}$. This generates relative velocities in the y direction between the adjacent boxes in the x direction, so that the y coordinate in Equation (4) depends on the MHD time T. As the initial condition for MHD simulations, we set the plasma to be of uniform density ${\rho }_{0}$. The initial plasma has a vertical magnetic field, ${B}_{0}=(0,\;0,\;{B}_{z})$ with ${\beta }_{{\rm{pl}}}={10}^{4}$, where ${\beta }_{{\rm{pl}}}=8\pi {\rho }_{0}{c}_{{\rm{s}}}^{2}/{B}_{0}^{2}$ (${c}_{{\rm{s}}}$ is the sound speed). We perform MHD simulations with the above setup using a modified version of an MHD code used for investigating MRI turbulence in protoplanetary disks (Suzuki & Inutsuka 2009; Suzuki et al. 2010) which adopts a second-order Godunov-CMoCCT scheme.

The most unstable wavelength of MRI, ${\lambda }_{{\rm{cr}}}$, is represented as (e.g., Balbus & Hawley 1998; Sano et al. 2004)

$$\begin{eqnarray}&&{\lambda }_{{\rm{cr}}}=\displaystyle \frac{8\pi }{\sqrt{15}}\displaystyle \frac{{v}_{{\rm{A,0}}}}{{{\rm{\Omega }}}_{{\rm{K}}}}\simeq 0.0649H{\left(\displaystyle \frac{{\beta }_{{\rm{pl}}}}{{10}^{4}}\right)}^{-1/2},\end{eqnarray} \tag{ 7 }$$

where ${v}_{{\rm{A,0}}}={B}_{0}/\sqrt{4\pi {\rho }_{0}}$ is the Alfvén velocity and $H=\sqrt{2}{c}_{{\rm{s}}}/{{\rm{\Omega }}}_{{\rm{K}}}$ is the scale height. For the particle simulation, the box size and mesh number are fixed as $({L}_{x}/{N}_{x},\;{L}_{y}/{N}_{y},\;{L}_{z}/{N}_{z})$ = (2H/128, 4H/256, 2H/128), for which ${\lambda }_{{\rm{cr}}}$ is resolved by several grids. We briefly discuss the effect of the box size in the next subsection.

The nonlinear growth of MRI generates strongly turbulent fields. After several orbital periods, a quasi-steady state is realized. The strength of the turbulent fields is indicated by the α parameter introduced by Shakura & Sunyaev (1973), which is expressed as

$$\begin{eqnarray}&&\alpha =\displaystyle \frac{{w}_{{xy}}}{{P}_{0}}=\displaystyle \frac{1}{{P}_{0}}{\left\langle \rho \delta {V}_{x}\delta {V}_{y}-\displaystyle \frac{{B}_{x}{B}_{y}}{4\pi }\right\rangle }_{{\rm{v}}}\end{eqnarray} \tag{ 8 }$$

where w_xy is the xy component of the stress tensor, ${\boldsymbol{\delta }}{\boldsymbol{V}}$ is the turbulent velocity, ${P}_{0}={\rho }_{0}{c}_{{\rm{s}}}^{2}$ is the initial pressure, and $\langle {\rangle }_{{\rm{v}}}$ denotes the volume average quantities. The turbulent velocity is defined as ${\boldsymbol{\delta }}{\boldsymbol{V}}={{\boldsymbol{V}}}_{{\rm{fl}}}-{{\boldsymbol{V}}}_{{\rm{shear}}}$. Figure 1 shows the evolution of α for two models of different box sizes, $({L}_{x},\;{L}_{y},\;{L}_{z})$ = (2H, 4H, 2H) and (0.8H, 1.6H, 0.8H). We can see that for the model with the larger box, α is almost constant ($\sim 0.02-0.03$) for $T\gtrsim 7\;{T}_{{\rm{rot}}}$, where ${T}_{{\rm{rot}}}=2\pi /{{\rm{\Omega }}}_{{\rm{K}}}$ is the rotation time. On the other hand, α varies between 0.01 and 0.06 for the smaller box. This indicates that the smaller the box size is, the stronger are the fluctuations in the quasi-steady state (Bodo et al. 2008; Jiang et al. 2013). Hereafter, we will use the data set of the larger box model in order to avoid the effect of the temporal fluctuation of turbulent fields. The initial value of ${\beta }_{{\rm{pl}}}$ also affects the saturation level of MRI, but we only treat the data set for ${\beta }_{{\rm{pl}}}={10}^{4}$ as the first step of our study.

Zoom In Zoom Out Reset image size

We tabulate the components of the averaged magnetic field $\langle {B}_{i}^{2}{\rangle }_{{\rm{v}}}$, turbulent velocity $\langle \delta {V}_{i}^{2}{\rangle }_{{\rm{v}}}$, and electric field $\langle {E}_{i}^{2}{\rangle }_{{\rm{v}}}$ in Table 1, for which we use MHD data of $T=20\;{T}_{{\rm{rot}}}$. The electric field is computed from the ideal MHD condition:

$$\begin{eqnarray}&&{\boldsymbol{E}}=-\displaystyle \frac{{{\boldsymbol{V}}}_{{\rm{fl}}}\times {\boldsymbol{B}}}{c}.\end{eqnarray} \tag{ 9 }$$

These fields satisfy $\langle {B}_{y}^{2}{\rangle }_{{\rm{v}}}$ > $\langle {B}_{x}^{2}{\rangle }_{{\rm{v}}}$ ≳ $\langle {B}_{z}^{2}{\rangle }_{{\rm{v}}}$, $\langle \delta {V}_{y}{\rangle }_{{\rm{v}}}\gg \langle \delta {V}_{x}{\rangle }_{{\rm{v}}}$ ≳ $\langle \delta {V}_{z}{\rangle }_{{\rm{v}}}$, and $\langle {E}_{z}^{2}{\rangle }_{{\rm{v}}}\gtrsim \langle {E}_{x}^{2}{\rangle }_{{\rm{v}}}\gg \langle {E}_{y}^{2}{\rangle }_{{\rm{v}}}$. The magnetic field is stretched by the shear motion, so that B_y is stronger than the other directions. The shear velocity is faster than turbulent velocity, $\langle {V}_{{\rm{shear}}}^{2}{\rangle }_{{\rm{v}}}\gt \langle \delta {V}_{i}^{2}{\rangle }_{{\rm{v}}}$. The electric field is then mainly generated by the shear motion, so that $\langle {E}_{y}^{2}{\rangle }_{{\rm{v}}}$ is weaker than $\langle {E}_{x}^{2}{\rangle }_{{\rm{v}}}$ and $\langle {E}_{z}^{2}{\rangle }_{{\rm{v}}}$.

Table 1.  Mean Values of Background Fields

$\langle {B}_{x}^{2}{\rangle }_{{\rm{v}}}/(8\pi {P}_{0})$	$6.58\times {10}^{-3}$
$\langle {B}_{y}^{2}{\rangle }_{{\rm{v}}}/(8\pi {P}_{0})$	$3.62\times {10}^{-2}$
$\langle {B}_{z}^{2}{\rangle }_{{\rm{v}}}/(8\pi {P}_{0})$	$3.23\times {10}^{-3}$
$\langle \delta {V}_{x}^{2}{\rangle }_{{\rm{v}}}{/c}_{{\rm{s}}}^{2}$	$1.79\times {10}^{-2}$
$\langle \delta {V}_{y}^{2}{\rangle }_{{\rm{v}}}{/c}_{{\rm{s}}}^{2}$	$2.66\times {10}^{-2}$
$\langle \delta {V}_{z}^{2}{\rangle }_{{\rm{v}}}{/c}_{{\rm{s}}}^{2}$	$1.10\times {10}^{-2}$
$\langle {E}_{x}^{2}{\rangle }_{{\rm{v}}}/(8\pi {P}_{0})$	$2.63\times {10}^{-5}$
$\langle {E}_{y}^{2}{\rangle }_{{\rm{v}}}/(8\pi {P}_{0})$	$5.07\times {10}^{-7}$
$\langle {E}_{z}^{2}{\rangle }_{{\rm{v}}}/(8\pi {P}_{0})$	$5.45\times {10}^{-5}$

Download table as:  ASCIITypeset image

We calculate the power spectra of turbulent magnetic, velocity, and electric fields. First, we take the Fourier transformation of the variables

$$\begin{eqnarray}&&{{\boldsymbol{B}}}_{{\boldsymbol{k}}}=\int {dxdydz}\mathrm{exp}(i{\boldsymbol{k}}\cdot {\boldsymbol{x}}){\boldsymbol{B}}(x,y,z).\end{eqnarray} \tag{ 10 }$$

Then, we define the one-dimensional power spectrum as

$$\begin{eqnarray}&&{B}_{{k}_{x}}^{2}=\int {{dk}}_{y}{{dk}}_{z}| {{\boldsymbol{B}}}_{{\boldsymbol{k}}}{| }^{2},\end{eqnarray} \tag{ 11 }$$

where we write the definition of the power spectrum of magnetic field for the x direction. The power spectra of the other variables and the other directions are also defined in the same manner. The power spectra of ${\boldsymbol{B}}$, ${\boldsymbol{\delta }}{\boldsymbol{V}}$, and ${\boldsymbol{E}}$ for the x and y directions of the data of $T=20\;{T}_{{\rm{rot}}}$ are shown in Figure 2. The Nyquist wavenumber is ${k}_{{\rm{N}}}\simeq 201/H$ for all of the directions, and then the small-scale turbulence ${k}_{i}\gtrsim {k}_{{\rm{N}}}/4\simeq 50/H$, shown by the vertical dotted line, suffers from numerical dissipation. On the other hand, the large-scale turbulence for ${k}_{i}\lesssim 10/H$ is probably affected by the injection process of turbulence. This scale becomes larger as MRI grows in the nonlinear regime, although it is given by Equation (7) in the linear regime, whose scale is close to the vertical dashed line in Figure 2 (e.g., Sano & Inutsuka 2001). Here, we consider $10\lesssim {k}_{i}H\lesssim 50$ as the inertial range. The shapes of the power spectra in the inertial range are similar for all of the variables. The power spectra have anisotropy. The spectra for the x and z directions are almost the same, which is almost consistent with the Kolmogorov turbulence, ${B}_{{k}_{x}}^{2}\sim {B}_{{k}_{z}}^{2}\propto {k}_{x}^{-5/3}$. On the other hand, ${B}_{{k}_{y}}^{2}\propto {k}_{y}^{-2}$, which is steeper than those for the other directions. Since the magnetic field is almost parallel to the y direction, these properties are consistent with the theory of incompressible MHD turbulence (Goldreich & Sridhar 1995), which is also confirmed by the previous simulations (e.g., Cho & Lazarian 2003). These properties of the turbulent fields are unchanged even for snapshot data of different MHD times. Note that in order to see the waves moving with the background shear, one should change the coordinates and wavenumbers according to Hawley et al. (1995). Although this slightly modifies the inertial range of the power spectra, we do not have to deal with it because the motions of the test particles are calculated in the coordinates of the shearing box (see Section 3.1).

Zoom In Zoom Out Reset image size

We calculate the orbits of CRs as test particles in the turbulence generated in the MHD simulations. The equation of motion for a test particle is

$$\begin{eqnarray}&&\displaystyle \frac{d{\boldsymbol{p}}}{{dt}}=e\left({{\boldsymbol{E}}}_{{\rm{mod}}}+\displaystyle \frac{{\boldsymbol{v}}\times {{\boldsymbol{B}}}_{{\rm{mod}}}}{c}\right),\end{eqnarray} \tag{ 12 }$$

where ${\boldsymbol{p}}=\gamma m{\boldsymbol{v}}$ and ${\boldsymbol{v}}$ are the momentum and velocity of the test particle, respectively (γ is the Lorentz factor of the particle). The tidal and Coriolis forces are included in ${E}_{{\rm{mod}}}$ and ${B}_{{\rm{mod}}}$, respectively (Hoshino 2013),

$$\begin{eqnarray}&&{{\boldsymbol{E}}}_{\mathrm{mod}}={\boldsymbol{E}}+\displaystyle \frac{3\gamma m}{e}{{\rm{\Omega }}}_{{\rm{K}}}^{2}x{{\boldsymbol{e}}}_{x},\end{eqnarray} \tag{ 13 }$$

$$\begin{eqnarray}&&{{\boldsymbol{B}}}_{\mathrm{mod}}={\boldsymbol{B}}+\displaystyle \frac{2\gamma m}{e}{{\boldsymbol{\Omega }}}_{{\rm{K}}}.\end{eqnarray} \tag{ 14 }$$

Note that the tidal and Coriolis forces are much weaker than the electromagnetic force, and they have little effect on the orbits of CRs in our simulations. We assume CRs as protons and neglect radiative and Coulomb losses. We use the Boris method (e.g., Birdsall & Langdon 1991) to integrate this equation, which is often used in PIC simulations.

Since the MHD data are defined only on the grid points, it is necessary to interpolate the data to the particle position. The MHD data consist of ${\boldsymbol{B}}$ and ${{\boldsymbol{V}}}_{{\rm{fl}}}$. First, we compute ${\boldsymbol{B}}$ and ${{\boldsymbol{V}}}_{{\rm{fl}}}$ at the particle position through linear interpolation. Then, we calculate ${\boldsymbol{E}}$ according to Equation (9). This procedure guarantees ${\boldsymbol{E}}\cdot {\boldsymbol{B}}=0$, which is necessary to avoid artificial acceleration due to interpolation (cf., Roh et al. 2016). In reality, ${\boldsymbol{B}}$ and ${\boldsymbol{E}}$ evolve on an MHD timescale. However, as the first step, we use the snapshot MHD data for simplicity.

As the initial conditions, the particle positions and velocities are given such that a uniform and isotropic distribution is realized in the calculation box (hereafter, we call this box the initial box). We set monoenergetic CRs whose gyro radii ${r}_{{\rm{gyro,0}}}$ satisfy

$$\begin{eqnarray}&&{r}_{{\rm{gyro,0}}}\equiv \displaystyle \frac{{E}_{p,0}}{{e\langle B\rangle }_{{\rm{v}}}}={\epsilon }_{{\rm{gyro}}}{\rm{\Delta }}x,\end{eqnarray} \tag{ 15 }$$

where ${E}_{p,0}$ is the initial energy of CRs, ${\epsilon }_{{\rm{gyro}}}$ is a parameter of the initial energy of CRs, and ${\rm{\Delta }}x={L}_{x}/{N}_{x}$ is the mesh size. The initial Lorentz factor is then represented as

$$\begin{eqnarray}&&{\gamma }_{0}=\displaystyle \frac{e\langle B{\rangle }_{{\rm{v}}}{\epsilon }_{{\rm{gyro}}}}{{{mc}}^{2}}\displaystyle \frac{{L}_{x}}{{N}_{x}}.\end{eqnarray} \tag{ 16 }$$

The value of ${\gamma }_{0}$ depends on the models of the accretion flows. We tabulate ${\gamma }_{0}$ in Table 2, supposing that the accretion flow has ${c}_{{\rm{s}}}$ ≃ 2.2 × 10⁹ cm s⁻¹, ${{\rm{\Omega }}}_{{\rm{K}}}$ ≃ 4.4 × 10⁻⁶ s⁻¹, and ${n}_{0}={\rho }_{0}/m\simeq 9.0\times {10}^{7}\;{\mathrm{cm}}^{-3}$. These values correspond to the self-similar solution for RIAFs (Yuan & Narayan 2014) with $m={10}^{8}$ (black hole mass in unit of M_⊙), ${\dot{m}}_{{\rm{BH}}}={10}^{-3}$ (mass accretion rate in unit of the Eddington ratio), r = 30 (disk radius in unit of the Schwarzschild radius), s = 0 (parameter of mass loss rate), and $\alpha =0.03$.⁴ In reality, the energies of these CRs are so high that they probably escape from the flow in the vertical direction. It is desirable to calculate CRs of much lower energies, such as ${\gamma }_{0}\sim 10-100$. However, the gyro-scale of such low-energy CRs is too small compared to the grid scale of the MHD simulations. These particles do not significantly change their energies in interactions with the large-scale turbulence that is likely to determine the maximum energy of the CRs in accretion flows. To investigate the interaction between CRs and large-scale turbulence, we calculate the orbital motions of the high-energy CRs by prohibiting them from escaping in the vertical direction. Since these CRs are ultrarelativistic, we can write ${E}_{p,0}={p}_{0}c$, where p₀ is the initial momentum of CRs.

Table 2.  Model Parameters and Physical Quantities

Model	T	${\epsilon }_{{\rm{gyro}}}$	${C}_{{\rm{esc}}}$	${\gamma }_{0}$	${t}_{{\rm{end}}}$ a	$\delta t$ a	Dxb	${D}_{y}$ b	${D}_{z}$ b	q	D0c	A
A1	$20{T}_{{\rm{rot}}}$	4	2	$3.4\times {10}^{8}$	416	50	2.3	17	1.6	2.38	$1.59\times {10}^{-4}$	0.30
A2	$20{T}_{{\rm{rot}}}$	1	2	$8.5\times {10}^{7}$	5628	200	2.4	24	1.6	1.91	$3.17\times {10}^{-5}$	0.25
A3	$20{T}_{{\rm{rot}}}$	8	2	$6.8\times {10}^{8}$	94	40	2.0	16	2.0	2.79	$5.38\times {10}^{-4}$	0.25
B1	$15{T}_{{\rm{rot}}}$	4	2	$3.5\times {10}^{8}$	447	50	2.2	19	1.5	2.38	$1.50\times {10}^{-4}$	0.31
C1	$25{T}_{{\rm{rot}}}$	4	2	$3.7\times {10}^{8}$	438	50	2.2	17	1.6	2.46	$1.45\times {10}^{-4}$	0.28
X1	$20{T}_{{\rm{rot}}}$	4	$\infty$	$3.4\times {10}^{8}$	$\infty$	200	3.1	19	2.0	0.969	$5.83\times {10}^{-5}$	0.54
X2	$20{T}_{{\rm{rot}}}$	1	$\infty$	$8.5\times {10}^{7}$	$\infty$	200	2.5	20	1.5	1.31	$2.88\times {10}^{-5}$	0.46
X3	$20{T}_{{\rm{rot}}}$	8	$\infty$	$6.8\times {10}^{8}$	$\infty$	200	3.5	21	3.0	0.654	$6.44\times {10}^{-5}$	0.58

Notes.

^aIn unit of ${t}_{{\rm{gyro,0}}}$. ^bIn unit of ${D}_{{\rm{Bohm}}}$. ^cIn unit of ${\bar{D}}_{p}$.

Download table as:  ASCIITypeset image

Due to Equations (4)–(6), we can continue to calculate the orbits of CRs outside of the initial box. Because of the relative velocities between boxes, we should convert the energies and momenta of CRs crossing the box boundary in the x direction by Lorentz transformation as

$$\begin{eqnarray}&&{E}_{p}^{\prime }={\rm{\Gamma }}({E}_{p}-{\beta }_{{\rm{box}}}{{cp}}_{y}),\end{eqnarray} \tag{ 17 }$$

$$\begin{eqnarray}&&{p}_{y}^{\prime }={\rm{\Gamma }}\left({p}_{y}-{\beta }_{{\rm{box}}}\displaystyle \frac{{E}_{p}}{c}\right),\end{eqnarray} \tag{ 18 }$$

$$\begin{eqnarray}&&\ {p}_{x}^{\prime }={p}_{x},{p}_{z}^{\prime }={p}_{z},\end{eqnarray} \tag{ 19 }$$

where ${\beta }_{{\rm{box}}}=\mp 1.5{{\rm{\Omega }}}_{{\rm{K}}}{L}_{x}/c$ is the relative velocity between boxes (we use the − sign for the CRs crossing $x=+{L}_{x}/2$ boundary, and vice versa) and ${\rm{\Gamma }}={(1-{\beta }_{{\rm{box}}}^{2})}^{-1/2}$. These Γ and ${\beta }_{{\rm{Box}}}$ depend on ${c}_{{\rm{s}}}$ and ${L}_{x}/H$. We get ${\beta }_{{\rm{box}}}\simeq 0.306$ and ${\rm{\Gamma }}-1$ ≃ 5.01 × 10⁻² for ${c}_{{\rm{s}}}$ ≃ 2.2 × 10⁹ cm s⁻¹ and ${L}_{x}/H=2$. Since ${\boldsymbol{v}}={c}^{2}{\boldsymbol{p}}/{E}_{p}$, v_x and v_z are also changed by the transformation. We solve Equation (12) in the rest frame of a box where each CR is located.

This Lorentz transformation causes an artificial acceleration for sufficiently energetic CRs. This acceleration happens when the energetic CRs continuously cross several boxes without changing their directions of motion (see Section 4.1 for details). In order to avoid this, we remove those CRs that have crossed a boundary at $x={x}_{{\rm{esc}}}$ as escape particles. We define ${x}_{{\rm{esc}}}=({C}_{{\rm{esc}}}+0.5){L}_{x}$, where ${C}_{{\rm{esc}}}$ is the escape parameter, and set ${C}_{{\rm{esc}}}=2$ such that the relative Lorentz factor between the escape boundaries at $x=-{x}_{{\rm{esc}}}$ and ${x}_{{\rm{esc}}}$ is less than two. Although the Lorentz contraction may also affect the length scale that CRs travel, we ignore its effect for simplicity.

We note that the boundary condition for CRs described above is different from that for the MHD simulations. The latter use the Galilean transformation, since the fluid velocity is non-relativistic, ${\rm{\Gamma }}-1\ll 1$. We cannot use the Galilean transformation as the boundary condition for CRs because the velocity of CRs after the Galilean transformation could exceed the speed of light, ${v}_{y}^{\prime }={v}_{y}+1.5{{\rm{\Omega }}}_{{\rm{K}}}{L}_{x}\sim c+1.5{{\rm{\Omega }}}_{{\rm{K}}}{L}_{x}\gt c$.

We set the time step as ${\rm{\Delta }}t={\rm{min}}({\rm{\Delta }}{t}_{{\rm{gyro}}},\;{\rm{\Delta }}{t}_{{\rm{cell}}})$, where ${\rm{\Delta }}{t}_{{\rm{gyro}}}={C}_{{\rm{safe}}}{E}_{p,0}/({{ceB}}_{{\rm{max}}})$ and ${\rm{\Delta }}{t}_{{\rm{cell}}}={C}_{{\rm{safe}}}{\rm{\Delta }}x/c$. Note that the velocities of CRs are always almost the speed of light because we focus on the ultrarelativistic particles. We use the maximum value of the magnetic field in the box ${B}_{{\rm{max}}}$ to estimate ${\rm{\Delta }}{t}_{{\rm{gyro}}}$, and set ${C}_{{\rm{safe}}}=0.01$.

We show the results of the simulations and discuss the evolution of the distribution function. The parameter sets are tabulated in Table 2. The letters A, B, and C represent the different times T of snapshot data, and we use the numbers 1, 2, and 3 to distinguish between the initial energy ${\epsilon }_{{\rm{gyro}}}$. Group X will be discussed in Section 4.2. We use ${N}_{p}={2}^{15}=32768$ CR particles and calculate their orbits until half of CRs escape from the system, the times of which are tabulated in Table 2 in units of the initial gyro-period, ${t}_{{\rm{gyro,0}}}=2\pi {r}_{{\rm{gyro,0}}}/c$. For all the models, the CRs randomly gain or lose their energies through interactions with turbulent fields, and diffuse in both the configuration and momentum spaces.

3.2.1. Lab Frame, Box Frame, and Shear Frame

There are three frames for evaluating the positions and momenta of CRs. One is the rest frame of the initial box where the CRs are initially located (lab frame), another is the rest frame of a box where each CR is located at the evaluation time (box frame), and the other is the rest frame of the MHD fluid element in the mean flow (i.e., the unperturbed flow) in each box (shear frame).

Figure 3 shows the evolution of the energy of a CR in the shear frame (thick solid line) and the box frame (dashed line). When we measure the energy in the box frame, the CR energy jumps due to the Lorentz transformation. In Figure 3, we can see two jumps of energy at $t\sim 240\;{t}_{{\rm{gyro,0}}}$ and $t\sim 243\;{t}_{{\rm{gyro,0}}}$, which coincide with the CR's crossing the box boundary (shown by the dotted lines) in the x direction. The CR position is represented by the thin solid line. On the other hand, the energy measured in the shear frame does not have such jumps but smoothly evolves with time. Since the box boundaries are not special surfaces in nature, we use the shear frame for discussing the evolution of CR energy unless otherwise noted.

Zoom In Zoom Out Reset image size

As mentioned above, we find that CRs randomly gain or lose a small amount of energy through their interaction with turbulent fields. A small fraction of CRs continuously gain (lose) energies, so that they can reach several times higher (lower) energies than their initial energies. Figure 4 shows the long-term evolution of the energies of such CRs in the shear frame. The most energetic CR at $t=400\;{t}_{{\rm{gyro,0}}}$ has about six times higher energy than the average value. This gradual change of particle energy implies that there is no "hot spot," where CRs gain energy efficiently, in the MRI turbulence. Note that the average energy, shown as the dotted line in Figure 4, is gradually increasing. This is consistent with the quasi-linear theory of stochastic acceleration (e.g., Stawarz & Petrosian 2008).

Zoom In Zoom Out Reset image size

3.2.2. Evolution in Configuration Space

The diffusion coefficient in configuration space is estimated to be

$$\begin{eqnarray}&&{D}_{{x}_{i}}=\displaystyle \frac{\langle \delta {x}_{i}^{2}{\rangle }_{p}}{2\delta t},\end{eqnarray} \tag{ 20 }$$

where $\delta {x}_{i}={x}_{i}(t+\delta t)-{x}_{i}(t)$ is the displacement of a particle (${x}_{i}=x,\;y$, or $\;z$), $\delta t$ is the time span for estimating the diffusion coefficient, and $\langle {\rangle }_{p}$ denotes the average over the CRs. The estimate of the diffusion coefficient is affected by the gyro-motion for too short $\delta t$, while it is affected by excluding the escaped CRs for too long $\delta t$. We calculate ${D}_{{x}_{i}}$ for several values of $\delta t$ and choose suitable $\delta t$ for estimating the diffusion coefficient, which is tabulated in Table 2.

We show the temporal evolution of ${D}_{x},{D}_{y},$ and D_z for model A1 in the upper panel of Figure 5. We normalize ${D}_{{x}_{i}}$ by the Bohm coefficient,

$$\begin{eqnarray}&&{D}_{{\rm{Bohm}}}=\displaystyle \frac{1}{3}{r}_{{\rm{gyro}}}c=\displaystyle \frac{\langle p{\rangle }_{p}{c}^{2}}{3{e\langle B\rangle }_{{\rm{v}}}}.\end{eqnarray} \tag{ 21 }$$

The diffusion coefficients in configuration space have anisotropy. Since the magnetic field is stretched in the y direction, CRs easily move in the y direction. Thus, ${D}_{y}\gt {D}_{x}\sim {D}_{z}$ is satisfied. The diffusion coefficients are almost constant in time, and so super diffusion is not observed in our simulations (see Xu & Yan 2013; Lazarian & Yan 2014; Roh et al. 2016).

Zoom In Zoom Out Reset image size

We do not discuss the diffusion coefficient in the direction parallel or perpendicular to the averaged magnetic field. Although we can define the volume-averaged magnetic field $\langle {\boldsymbol{B}}{\rangle }_{{\rm{v}}}$, it does not represent the local (in the scale of gyro radius) direction of the magnetic field due to the strong turbulence.

Dividing CRs into a few momentum bins, we estimate the momentum dependence of ${D}_{{x}_{i}}$. We use 20 momentum bins of equal interval in the logarithmic space for $0.1\leqslant p/{p}_{0}\leqslant 20$. To reduce the statistical fluctuation, we use only those momentum bins which have more than 100 CRs. In the lower panel of Figure 5, the momentum dependence of ${D}_{{x}_{i}}$ is shown at $t=400{t}_{{\rm{gyro,0}}}$. We can see that ${D}_{{x}_{i}}\propto p$, which is the same dependence of the Bohm coefficient. These are common features for all of the models. The diffusion coefficients in configuration space for the other models are tabulated in Table 2. All of the models have similar values of the diffusion coefficients, ${D}_{y}\sim 20{D}_{{\rm{Bohm}}}$ and ${D}_{x}\sim {D}_{z}\sim 2{D}_{{\rm{Bohm}}}$.

3.2.3. Evolution in Momentum Space

We discuss the evolution of the momentum distribution function $f({\boldsymbol{p}})$ of CRs. Figure 6 shows the momentum distribution of each direction ${dN}/{{dp}}_{i}$, normalized as $\int ({dN}/{{dp}}_{i}){{dp}}_{i}=1$ (${p}_{i}={p}_{x},\;{p}_{y},$ or p_z). CRs diffuse in all of the directions in momentum space. From the upper panel, it can be seen that the momentum distribution is isotropic in the shear frame. The turbulent fields isotropize the motion of CRs in the shear frame. The dispersion of the momentum distribution monotonically increases with time. The momentum distribution in the box frame is almost the same as that in the shear frame because the relative shear velocity inside of a box is sufficiently smaller than the particle velocities. On the other hand, in the lab frame, the distribution is anisotropic due to the velocity difference between the boxes. The dispersion in the y direction is larger than in the other directions. In this section, we will use the shear frame when discussing the evolution of the distribution function. Note that the isotropy of the momentum distribution does not conflict with the anisotropic diffusion in configuration space. This situation is possible when the timescale of pitch angle scattering is much shorter than the diffusion timescale in configuration space.

Zoom In Zoom Out Reset image size

Since the momentum distribution is isotropic, we can write $f({\boldsymbol{p}})=f(p)$, where $p=| {\boldsymbol{p}}|$. Below, we show that the evolution of f(p) obtained by our simulations is well described by the diffusion equation in momentum space:

$$\begin{eqnarray}&&\displaystyle \frac{\partial f(p)}{\partial t}=\displaystyle \frac{1}{{p}^{2}}\displaystyle \frac{\partial }{\partial p}\left({p}^{2}{D}_{p}\displaystyle \frac{\partial f(p)}{\partial p}\right)-\displaystyle \frac{f(p)}{{t}_{{\rm{esc}}}}+{\dot{f}}_{{\rm{inj}}},\end{eqnarray} \tag{ 22 }$$

where D_p is the diffusion coefficient in momentum space, ${t}_{{\rm{esc}}}={({C}_{{\rm{esc}}}+0.5)}^{2}{L}_{x}^{2}/(2{D}_{x})$ is the escape time, and ${\dot{f}}_{{\rm{inj}}}={N}_{p}\delta (t)\delta (p-{p}_{0})$ is the injection term. We write $p={E}_{p}/c=\gamma {mc}$ because CRs are ultrarelativistic. Equation (22) is expected to describe the evolution of the distribution function when the fractional energy change per scattering is small.

To calculate the evolution of f(p) using Equation (22), we estimate the diffusion coefficient in momentum space D_p from our simulation results. The diffusion coefficient is likely to be represented as

$$\begin{eqnarray}&&{D}_{p}=A\displaystyle \frac{\langle \delta {p}^{2}{\rangle }_{p}}{\delta t},\end{eqnarray} \tag{ 23 }$$

where $\delta p=p(t+\delta t)-p(t)$ and A is the numerical factor. We estimate D_p using the momentum bins defined in the previous subsection. Although $A=1/6$ when we consider the isotropic random walk in three-dimensional space, the values of A are not obvious due to the dependence of D_p on p. In fact, the values of A are different among the models (see Table 2). We plot $\langle \delta {p}^{2}{\rangle }_{p}/\delta t$ for model A1 in Figure 7. We can fit $\langle \delta {p}^{2}{\rangle }_{p}/\delta t$ by a power-law function, $\langle \delta {p}^{2}{\rangle }_{p}/\delta t\propto {(p/{p}_{0})}^{q}$, for all of the models. The resultant q is tabulated in Table 2.

Zoom In Zoom Out Reset image size

We compute A as follows. If the temporal evolution of f(p) is described by Equation (22), then the temporal evolution of the averaged momentum $\langle \dot{p}{\rangle }_{p}$ is written as (see Becker et al. 2006)

$$\begin{eqnarray}&&\langle \dot{p}{\rangle }_{p}=\displaystyle \frac{d\langle p{\rangle }_{p}}{{dt}}=\displaystyle \frac{4\pi }{{N}_{p}}\int {p}^{3}\displaystyle \frac{\partial f}{\partial t}{dp}=\displaystyle \frac{4\pi }{{N}_{p}}\int p\displaystyle \frac{\partial }{\partial p}\left({p}^{2}{D}_{p}\displaystyle \frac{\partial f}{\partial p}\right),\end{eqnarray} \tag{ 24 }$$

where we use Equation (22) and ignore the terms ${\dot{f}}_{{\rm{inj}}}$ and $f(p)/{t}_{{\rm{esc}}}$. Integrating this equation by part twice, we obtain

$$\begin{eqnarray}&&\langle \dot{p}{\rangle }_{p}=\displaystyle \frac{4\pi }{{N}_{p}}\int f\displaystyle \frac{\partial }{\partial p}({p}^{2}{D}_{p}).\end{eqnarray} \tag{ 25 }$$

Writing ${D}_{p}={D}_{0}{(p/{p}_{0})}^{q}$, ${dN}/{dp}=4\pi {p}^{2}f$, and $p/{p}_{0}=\xi$, we obtain

$$\begin{eqnarray}&&\displaystyle \frac{{D}_{0}}{{p}_{0}^{2}}=\displaystyle \frac{{N}_{p}\langle \dot{p}{\rangle }_{p}}{2+q}\left(\int d\xi \displaystyle \frac{{dN}}{d\xi }{\xi }^{q-1}\right).\end{eqnarray} \tag{ 26 }$$

We compute D₀ taking the time average of this equation over the whole calculation time. We tabulate the values of D₀ and A in Table 2, where D₀ is normalized by $\bar{{D}_{p}}={p}_{0}^{2}/{t}_{{\rm{gyro,0}}}$. Since $\bar{{D}_{p}}/{D}_{0}={t}_{{\rm{accel}}}/{t}_{{\rm{gyro,0}}}$, where we define ${t}_{{\rm{accel}}}={p}_{0}^{2}/{D}_{0}$, we find that the acceleration time is about ${10}^{3}\mbox{--}{10}^{4}$ times longer than the gyro-period in our simulations.

We solve Equation (22) with D_p obtained using the above procedure, and compare the distribution functions calculated by Equation (22) to those obtained by the particle simulations. In Figure 8, we plot f(p) of $t=100\;{t}_{{\rm{gyro,0}}}$ and $t=400\;{t}_{{\rm{gyro,0}}}$ calculated by both the simulation and the diffusion equation for model A1. We find that the distribution functions of the diffusion equation are always in agreement with those of the particle simulation. We confirm that the diffusion equation reproduces the evolution of the distribution function in all of the models. We can conclude that the stochastic acceleration inside the accretion flows can be described by Equation (22).

Zoom In Zoom Out Reset image size

We also calculate the orbits of CRs with snapshots in different times, $T=15\;{T}_{{\rm{rot}}}$ (model B1) and $T=25\;{T}_{{\rm{rot}}}$ (model C1). The differences of the values in Table 2 among models A1, B1, and C1 are less than 15%, which indicates that our results are almost unchanged by the temporal fluctuation of the turbulent fields in our current setup of the particle simulations. However, we should note that ${t}_{{\rm{gyro,0}}}\simeq 6.36\times {10}^{-3}\;{T}_{{\rm{rot}}}$ in model A1, and then the calculation time of CR orbits is longer than the time step of the MHD simulation. Thus, in reality, CRs should move in dynamically evolving turbulence. Although the current simulations with the snapshots of different times provide almost the same results, it is unclear how much the dynamically evolving fields affect the orbits of CRs. This effect should be investigated as the next step in a separate paper.

Finally, we briefly mention the scales of turbulence interacting with CRs. We show the initial gyro-scales of CRs as the vertical dashed and dotted lines for models A1 and A3, respectively, in Figure 2. Except for model A3, the initial gyro-scales are in the dissipation scale of the MHD simulation. Although the turbulence in the dissipation scale with steep spectra is artificially realized, the stochastic acceleration of CRs through interaction with such turbulence is physical. Accelerated CRs can interact with the turbulence in the inertial range due to their larger gyro-scales. However, highly accelerated CRs suffer from artificial acceleration due to the Lorentz transformation (see Section 4.1). Thus, it is difficult to keep a large fraction of CRs in the inertial range. We need MHD data with much higher resolution to observe the CRs that continue to interact with the realistic turbulent fields.

As briefly discussed in Section 3.1, sufficiently energetic CRs in the shearing box simulation become runaway particles. When CRs of ${p}_{y}\gt 0$ (${p}_{y}\lt 0$) cross the $x=+{L}_{x}/2$ ($x=-{L}_{x}/2$) boundary, their energies in the box frame increase. If their gyro radii are large enough to cross several boxes without changing the direction of motion, the CRs continuously gain energies.

To observe the behavior of the runaway particles, we perform a simulation without an escape boundary. Figure 9 shows the temporal evolution of the energy for a runaway particle in the box frame. It is seen that the energy of the runaway particle increases stepwise every time it crosses the box boundaries in the x direction. We find that the runaway particle has an almost constant box crossing time, ${\rm{\Delta }}t={L}_{x}/{v}_{x}\simeq {L}_{x}/({\beta }_{x}c)$. The energy gain per crossing boundary is ${\rm{\Delta }}{E}_{p}={\rm{\Gamma }}({E}_{p}+{\beta }_{{\rm{Box}}}{{cp}}_{y})-{E}_{p}=\{{\rm{\Gamma }}(1+{\beta }_{{\rm{box}}}{\beta }_{y})-1\}{E}_{p}$. Thus,

$$\begin{eqnarray}&&\displaystyle \frac{{{dE}}_{p}}{{dt}}\approx \displaystyle \frac{{\rm{\Delta }}{E}_{p}}{{\rm{\Delta }}t}\simeq \{{\rm{\Gamma }}(1+{\beta }_{{\rm{box}}}{\beta }_{y})-1\}{\beta }_{x}\displaystyle \frac{{E}_{p}c}{{L}_{x}}.\end{eqnarray} \tag{ 27 }$$

Solving this equation, we obtain ${E}_{p}\propto \mathrm{exp}(t/{t}_{{\rm{grow}}})$, where ${t}_{{\rm{grow}}}\simeq \{{\rm{\Gamma }}(1+{\beta }_{{\rm{box}}}{\beta }_{y})-1\}{\beta }_{x}{L}_{x}/c\simeq 0.30{T}_{{\rm{rot}}}$. Here, we use ${\beta }_{{\rm{box}}}\simeq 0.306$, ${\beta }_{x}\simeq 0.31c$, and ${\beta }_{y}\simeq 0.92c$, which are obtained using the simulation result. This solution is consistent with the temporal evolution of the energy of the runaway particle, as shown in Figure 9. One should care about such runaway particles when using the shearing box approximation without the escape boundaries. A more realistic global simulation without a shearing box approximation does not suffer from this kind of unphysical effect. Particle simulations with global turbulent fields are desirable to obtain a solid conclusion, although it is computationally too expensive to perform with the required resolution.

Zoom In Zoom Out Reset image size

In this subsection, we discuss the physical interpretation of D_p obtained in our simulations. As seen below, we find that the shear acceleration works efficiently, although the origin of D_p is not fully understood.

The shear acceleration is expected to be efficient in accretion flows (e.g., Katz 1991; Subramanian et al. 1999). The shear motion seems to be important in our simulation because the electric field is mainly produced by shear motion. To understand its relevance clearly, we perform particle simulations assuming ${V}_{{\rm{shear}}}=0$ and ${\beta }_{{\rm{box}}}=0$. In this treatment, we do not deal with the Lorentz transformation for the orbital evolutions of CRs (Equations (17) and (18)), and the electric field induced by shear motion vanishes. Since the runaway particle does not appear due to ${\beta }_{{\rm{box}}}=0$, we set ${C}_{{\rm{esc}}}=\infty$. As a result, we find that the evolution of the distribution function is described by Equation (22), and D_p and ${D}_{{x}_{i}}$ are obtained in the same manner as in Section 3.2. The resultant values are tabulated in Table 2 as model group X. We can see that both q and D₀ for group X are lower than those for group A. The differences of q and D₀ between groups X and A are larger since the initial gyro radius is larger. This implies that shear acceleration works efficiently for particles with larger gyro radii and that it is a dominant process of particle acceleration in model A3.

For the shearing box, D_p by the shear acceleration is analytically estimated as (see Earl et al. 1988; Rieger & Duffy 2004; Ohira 2013)

$$\begin{eqnarray}&&{D}_{p,\mathrm{shear}}\approx \displaystyle \frac{1}{15}{\left(\displaystyle \frac{\partial {V}_{y}}{\partial x}\right)}^{2}\tau {p}^{2}=\displaystyle \frac{3}{20}{p}^{2}\tau {{\rm{\Omega }}}_{{\rm{K}}}^{2},\end{eqnarray} \tag{ 28 }$$

where τ is the collision time. Since we found that the diffusion coefficient in configuration space is proportional to the Bohm coefficient, we may assume

$$\begin{eqnarray}&&\tau \approx \displaystyle \frac{2{D}_{x}}{{v}_{x}^{2}}=\displaystyle \frac{\eta }{\pi }\left(\displaystyle \frac{p}{{p}_{0}}\right){t}_{{\rm{gyro,0}}}\end{eqnarray} \tag{ , 29 }$$

where $\eta ={D}_{x}/{D}_{{\rm{Bohm}}}$ and ${v}_{x}^{2}={c}^{2}/3$. Then, we can write

$$\begin{eqnarray}&&{D}_{p,\mathrm{shear}}\approx \displaystyle \frac{3\pi \eta }{5}{\left(\displaystyle \frac{{t}_{{\rm{gyro,0}}}}{{T}_{{\rm{rot}}}}\right)}^{2}{\left(\displaystyle \frac{p}{{p}_{0}}\right)}^{3}\bar{{D}_{p}}.\end{eqnarray} \tag{ 30 }$$

We find ${D}_{p,\mathrm{shear}}\propto {p}^{3}$ and ${D}_{p,\mathrm{shear}}/\bar{{D}_{p}}\simeq 6.1\times {10}^{-4}$ for $p={p}_{0}$, using the parameter set of model A3 ($\eta \simeq 2.0$ and ${t}_{{\rm{gyro,0}}}/{T}_{{\rm{rot}}}\simeq 1.3\times {10}^{-2}$). The momentum dependence and value of this estimate are almost consistent with D_p of the simulation result, which also shows that the shear acceleration is the dominant process for model A3.

Then, a remaining problem is the physical origin of D_p in model X, in which CRs are accelerated by the electric field induced only by the turbulent fluid motion with the velocity ${\boldsymbol{\delta }}{\boldsymbol{V}}$. According to the quasi-linear theory of gyro-resonance by isotropic Alfvén modes, the power index of D_p is the same as the index of the power spectrum of the turbulent magnetic field (e.g., Stawarz & Petrosian 2008), in which $q=5/3$ is expected for the Kolmogorov turbulence. For our turbulence, the power spectra are steeper for smaller scales (see Figure 2), and so this theory would predict higher values of q for lower p₀ models. Model group A shows the opposite trend, while model group X is qualitatively consistent with this feature. However, quantitatively, $q\sim 4$ for model X1 (large k due to low ${E}_{p,0}$) and $q\sim 5/3$ for model X3 (small k due to high ${E}_{p,0}$) are expected in this theory, which are far from the particle simulation results.

Our turbulence has anisotropic power spectra (see Figure 2). Yan & Lazarian (2002) investigate the scattering frequency of such anisotropic MHD turbulence and show that the scattering frequency by the Alfvén modes is sufficiently suppressed compared to the case of isotropic turbulences. Instead, gyro-resonance by fast modes plays a crucial role in scattering CRs under MHD turbulence. The fast modes are almost isotropic and the index of the power spectrum is 1.5 (Cho & Lazarian 2003), which leads to q = 1.5 (see Cho & Lazarian 2006). Model X2 has a value close to q = 1.5, but it is unlikely that the fast modes have substantial power in the MRI turbulence because velocity fluctuations are sub-sonic in our turbulent field.

Lynn et al. (2014) claim that for sub-sonic MHD turbulence, the fast modes have only a small fraction of the turbulent power and slow modes dominate over the other modes in the acceleration of relativistic CRs. They consider the mirror force acting on CRs, and analytically derive ${D}_{p}\propto {p}^{2}$ for stochastic acceleration in the weakly compressible turbulence. In this case, q should always be equal to 2, which is often called the hard-sphere model, but this theory does not seem to be applicable to our model since the resultant q of our simulations is not close to 2 for both groups X and A.

As the first step of our study, we use the shearing box approximation for RIAFs, which was also used in the previous studies (e.g., Sharma et al. 2006; Hoshino 2013). However, the shearing box approximation has a few points which are at odds with RIAFs. First, the local approximation is not good for RIAFs because the aspect ratio, H/r, is close to unity due to the the high temperature. The centrifugal balance is also violated due to the strong pressure gradient force in RIAFs. The angular velocity is usually sub-Keplerian, which seems to weaken the shear acceleration. In addition, the radial advection of heat is important for the global structure of RIAFs. The high temperature makes viscous stress strong under the α prescription, which makes the radial velocity fast compared to the standard disk (e.g., Narayan & Yi 1994; Kimura et al. 2014). The global simulations can consistently take into account these effects.

Magnetic reconnection is supposed to play an important role in the saturation of MRI turbulence (e.g., Sano & Inutsuka 2001). Recent PIC simulations also suggest that dissipation by the reconnection of stretched magnetic fields is important for both the saturation of MRI turbulence and the generation of CRs (Hoshino 2015). However, since we solve a set of ideal MHD equations in which magnetic reconnection occurs due to numerical resistivity, we ignore the electric field generated by finite resistivity in the simulation of CRs. In reality, the electric field generated at the reconnection site is expected to enhance particle acceleration in MRI turbulence. To include this effect, the MHD simulation that includes the resistivity model based on collisionless plasma physics is required, which is beyond the scope of the present paper.

The effect of back reaction by CRs on the nonlinear evolution for MRI is not well understood. Kuwabara & Ko (2015) analyze the linear growth of MRI including the effect of CR diffusion. They show that CR diffusion barely affects the linear growth of MRI in shearing box calculations, even if the CR pressure is comparable to the thermal pressure. However, the effect of CRs on the saturation level of MRI turbulence has not yet been studied. In the case close to ideal MHD, the saturation level is supposed to be determined by the parasitic instability of the KH modes that induce the dissipation by magnetic reconnection (Goodman & Xu 1994; Sano & Inutsuka 2001; Sano et al. 2004; Pessah 2010). Since the growth rate of KH instability is modified by CRs if they have comparable energy density to the thermal particles (Suzuki et al. 2014), CRs probably affect the saturation level of MRI turbulence. To investigate this, three-dimensional MHD simulations including back reaction by CRs are necessary (see Bai et al. 2015).

Although we assume that the thermal component is described by the MHD formalism, it is unclear whether or not collisionless plasma is described by MHD. As an alternative method of describing the collisionless plasma in accretion flows, Quataert et al. (2002) used the anisotropic MHD formalism where the plasma has anisotropic pressure. In this formalism, as the magnetic field becomes stronger, the pressure anisotropy increases, which suppresses the growth of MRI (Sharma et al. 2006). In reality, the thermal particles are expected to be isotropized by some plasma instabilities, such as the mirror mode, so that the pressure anisotropy decreases. This allows MRI to grow in the nonlinear regime, which has been confirmed by recent PIC simulations with the shearing box approximation (Riquelme et al. 2012; Hoshino 2013, 2015). However, as noted in Section 1, there is the scale gap between the gyro-scales of electrons and the scale heights of accretion flows. The hybrid simulations in which we treat protons as particles and electrons as a fluid may be useful for investigating larger-scale turbulences with kinetic effects (see Shirakawa & Hoshino 2014), but it is currently difficult to bridge the scale gap completely.