Acceleration and Escape Processes of High-energy Particles in Turbulence inside Hot Accretion Flows

TL;DR

This study uses magnetohydrodynamic simulations to explore how turbulence in hot accretion flows accelerates high-energy particles, revealing diffusion in energy space and superdiffusion in space, with implications for cosmic ray acceleration in active galactic nuclei.

Contribution

It provides a detailed analysis of particle acceleration mechanisms in hot accretion flows using simulations, highlighting the roles of turbulence and magnetic field configurations.

Findings

Particle energy diffusion follows a $D_     $

Superdiffusive behavior in configuration space indicates faster-than-normal radial displacement of particles.

High-energy particles can be accelerated up to 0.1-10 PeV in active galactic nuclei with hot accretion flows.

Abstract

We investigate acceleration and propagation processes of high-energy particles inside hot accretion flows. The magnetorotational instability (MRI) creates turbulence inside accretion flows, which triggers magnetic reconnection and may produce non-thermal particles. They can be further accelerated stochastically by the turbulence. To probe the properties of such relativistic particles, we perform magnetohydrodynamic simulations to obtain the turbulent fields generated by the MRI, and calculate orbits of the high-energy particles using snapshot data of the MRI turbulence. We find that the particle acceleration is described by a diffusion phenomenon in energy space with a diffusion coefficient of the hard-sphere type: $D_\epsilon\propto \epsilon^2$, where $\epsilon$ is the particle energy. Eddies in the largest scale of the turbulence play a dominant role in the acceleration process. On the other hand, the stochastic behaviour in configuration space is not usual diffusion but superdiffusion: the radial displacement increases with time faster than that in the normal diffusion. Also, the magnetic field configuration in the hot accretion flow creates outward bulk motion of high-energy particles. This bulk motion is more effective than the diffusive motion for higher energy particles. Our results imply that typical active galactic nuclei that host hot accretion flows can accelerate CRs up to $\epsilon\sim 0.1-10$ PeV.

Tables

Table 1. Table 1: Parameters and physical quantities for MHD simulations. Here, ⟨ X ⟩ 𝒱 subscript delimited-⟨⟩ 𝑋 𝒱 \langle X\rangle_{\mathcal{V}} indicates the volume-averaged values of 0.1 ≤ R / R c ≤ 0.6 0.1 𝑅 subscript 𝑅 𝑐 0.6 0.1\leq R/R_{c}\leq 0.6 .

run	A	B	C	D
$(N_r,N_{\theta },N_{\varphi })$ for	(640, 320, 768)	(384, 192, 512)	(256, 128, 384)	(960, 192, 512)
$(N_r,N_{\theta },N_{\varphi })$ for	(320, 320, 768)	(192, 192, 512)	(127, 128, 384)	(166, 192, 512)
$\Delta R_{i+1}/\Delta R_i$	1.006	1.01	1.015	1
${\alpha }_{\mathrm{SS},M}$	0.022	0.025	0.023	0.021
${\alpha }_{\mathrm{mag}}$	0.45	0.46	0.38	0.47
${⟨B_R^2⟩}_{𝒱}/{⟨B_{\varphi }^2⟩}_{𝒱}$	0.16	0.16	0.10	0.16
${⟨B^2⟩}_{𝒱}/{⟨8\pi P⟩}_{𝒱}$	0.050	0.053	0.061	0.044
${⟨B_R^2⟩}_{𝒱}/{⟨8\pi P⟩}_{𝒱}$	$6.7\times {10}^{-3}$	$7.0\times {10}^{-3}$	$5.5\times {10}^{-3}$	$5.7\times {10}^{-3}$
${⟨B_{\theta }^2⟩}_{𝒱}/{⟨8\pi P⟩}_{𝒱}$	$2.8\times {10}^{-3}$	$2.9\times {10}^{-3}$	$2.0\times {10}^{-3}$	$2.2\times {10}^{-3}$
${⟨B_{\varphi }^2⟩}_{𝒱}/{⟨8\pi P⟩}_{𝒱}$	0.040	0.043	0.053	0.036
${⟨c{E}_R^2⟩}_{𝒱}/{⟨8\pi P⟩}_{𝒱}$	$6.9\times {10}^{-3}$	$6.7\times {10}^{-3}$	$4.7\times {10}^{-3}$	$5.0\times {10}^{-3}$
${⟨c{E}_{\theta }^2⟩}_{𝒱}/{⟨8\pi P⟩}_{𝒱}$	0.015	0.015	0.011	0.012
${⟨c{E}_{\varphi }^2⟩}_{𝒱}/{⟨8\pi P⟩}_{𝒱}$	$5.6\times {10}^{-5}$	$7.4\times {10}^{-5}$	$6.8\times {10}^{-5}$	$5.5\times {10}^{-5}$

Equations

$$\displaystyle\frac{\partial\rho}{\partial T}+\mbox{\boldmath${\nabla}$}\cdot\left(\rho\mbox{\boldmath${V}$}\right)$$

$$\displaystyle\frac{\partial(\rho\mbox{\boldmath${V}$})}{\partial T}+\mbox{\boldmath${\nabla}$}\cdot\left(\rho\mbox{\boldmath${V}$}\mbox{\boldmath${V}$}-{\mbox{\boldmath${B}$}\mbox{\boldmath${B}$}\over 4\pi}+P^{*}\mathbb{I}\right)$$

$$\displaystyle\frac{\partial E_{\rm tot}}{\partial T}+\mbox{\boldmath${\nabla}$}\cdot\left[\left(E_{\rm tot}+P^{*}\right)\mbox{\boldmath${V}$}-{\mbox{\boldmath${B}$}\cdot\mbox{\boldmath${V}$}\over 4\pi}\mbox{\boldmath${B}$}\right]$$

$$\displaystyle\frac{\partial\mbox{\boldmath${B}$}}{\partial T}-\mbox{\boldmath${\nabla}$}\times\left(\mbox{\boldmath${V}$}\times\mbox{\boldmath${B}$}\right)$$

$$E_{\mathrm{tot}}=E_{\mathrm{th}}+\frac{1}{2}\rho V^2+\frac{B^2}{8\pi },$$

$$\frac{P}{\rho }=\frac{GM}{(n+1)R_c}\left[\frac{R_c}{R}-\frac{1}{2}{\left(\frac{R_c}{R\sin \theta }\right)}^2-\frac{1}{2d}\right],$$

$$⟨X{⟩}_{\mathcal{V}}=\frac{\int \rho X{R}^2\sin \theta dRd\theta d\varphi }{\int \rho R^2\sin \theta dRd\theta d\varphi }.$$

$${\alpha }_{\mathrm{SS},\mathrm{M}}=\frac{⟨B_r{B}_{\varphi }{⟩}_{\mathcal{V}}}{⟨4\pi P{⟩}_{\mathcal{V}}},$$

$${\alpha }_{\mathrm{mag}}=\frac{⟨2B_r{B}_{\varphi }{⟩}_{\mathcal{V}}}{⟨B^2{⟩}_{\mathcal{V}}}$$

$$\mbox{\boldmath${E}$}=-{\mbox{\boldmath${V}$}\times\mbox{\boldmath${B}$}\over c}.$$

$$⟨X{⟩}_{\mathcal{S}}(R)=\frac{\int \rho X\sin \theta d\theta d\varphi }{\int \rho \sin \theta d\theta d\varphi }.$$

$$⟨V_{\varphi }{⟩}_{\mathcal{L}}(R,\varphi )=\frac{\int \rho V_{\varphi }\varphi }{\int \rho d\varphi }.$$

$$X_m=\frac{1}{\sqrt{2\pi }}\int X\exp (-im\varphi )d\varphi ,$$

$${\mathcal{P}}_m=\frac{\int |X_m{|}^{2}RdRd\theta }{\int RdRd\theta },$$

$${d\mbox{\boldmath${p}$}\over dt}=e\left(\mbox{\boldmath${E}$}+{\mbox{\boldmath${v}$}\times\mbox{\boldmath${B}$}\over c}\right),$$

$${\mathcal{L}}_u=\chi R_s,$$

$${\mathcal{T}}_u=\sqrt{\frac{{\mathcal{L}}_u^3}{GM}},$$

$${\mathcal{M}}_u=\eta {\dot{M}}_{\mathrm{Edd}}{\mathcal{T}}_u,$$

$${\mathcal{L}}_u\simeq 1.5\times 10^{15}{M}_8{\chi }_{1.7}\mathrm{cm}$$

$${\mathcal{T}}_u\simeq 4.9\times 10^5{M}_8{\chi }_{1.7}^{3/2}\mathrm{s},$$

$${\mathcal{M}}_u\simeq 6.9\times 10^{30}{M}_8^2{\chi }_{1.7}^{3/2}{\eta }_0\mathrm{g},$$

$$\mbox{\boldmath${V}$}_{\rm bg}=V_{{\rm bul},\phi}\mbox{\boldmath${e}$}_{\phi},$$

$$\frac{\partial f}{\partial t}=\frac{1}{p^2}\frac{\partial }{\partial p}\left(p^2{D}_p\frac{\partial f}{\partial p}\right).$$

$$\frac{\partial N_{ϵ}}{\partial t}=\frac{\partial }{\partial ϵ}\left(D_{ϵ}\frac{\partial N_{ϵ}}{\partial ϵ}\right)-\frac{\partial }{\partial ϵ}\left(\frac{2D_{ϵ}}{ϵ}{N}_{ϵ}\right).$$

$${\mu }_{ϵ}=\frac{1}{N_{\mathrm{actv}}}\int N_{ϵ}ϵdϵ\approx {ϵ}_{\mathrm{ini}}+v_{{ϵ}_{\mathrm{ini}}}t$$

$${\sigma }_{ϵ}^2=\frac{1}{N_{\mathrm{actv}}}\int N_{ϵ}{ϵ}^{2}dϵ-{\mu }_{ϵ}^2\approx 2D_{{ϵ}_{\mathrm{ini}}}t$$

$$D_{ϵ,\mathrm{FTB}}$$

$$D_{ϵ,\mathrm{TTD}}$$

$$t_{\mathrm{acc}}\approx \frac{{ϵ}^2}{2D_{ϵ,\mathrm{FTB}}}\sim \frac{3}{2}{\left(\frac{c}{V_{R,\mathrm{tur}}}\right)}^2{t}_{\mathrm{int}}\sim 1.7\times 10^7{M}_8{\chi }_{1.7}^{2}s.$$

$${\mu }_{ϵ}\approx {ϵ}_{\mathrm{ini}}+2v_{{ϵ}_{\mathrm{ini}}}t.$$

$${\mu }_{\delta R}=\frac{1}{N_{\mathrm{actv}}}\sum _j\delta R_j,$$

$${\sigma }_{\delta R}^2=\frac{1}{N_{\mathrm{actv}}}\sum _j\delta R_j^2-{\mu }_{\delta R}^2,$$

$${\sigma }_{\delta R}^2\sim r_L^2{\left(\frac{\zeta t}{t_L}\right)}^{\xi },$$

$${\beta }_{\mathrm{CR},\mathrm{bul}}=\frac{c}{N_{\mathrm{actv}}{\mu }_{ϵ}}\int p_R\frac{d{N}_p}{d{p}_R}d{p}_R.$$

$${\beta }_{\mathrm{CR},\mathrm{bul}}\simeq 2.9\times 10^{-3}{ϵ}_{\mathrm{PeV}}^{0.62}{M}_8^{-0.31}{\chi }_{1.7}^{0.16}{\eta }_0^{-0.31},$$

Full text

††pubyear: 2018††pagerange: Acceleration and escape processes of high-energy particles in turbulence inside hot accretion flows–16

1 Introduction

Mass accretion onto a compact object powers broadband emissions from active galactic nuclei (AGN) and Galactic X-ray binaries. Hot accretion flows are formed when an accretion rate is sufficiently lower than the Eddington accretion rate, $\dot{M}_{\rm Edd}=L_{\rm Edd}/c^{2}$ (Ichimaru, 1977; Narayan & Yi, 1994), which are believed to be realized in our Galactic center (Sgr A*; Narayan et al., 1995; Manmoto et al., 1997), low-luminosity active galactic nuclei (LLAGNs; Nemmen et al., 2006, 2014), and X-ray binaries in the low-hard state (Esin et al., 1997; Yuan et al., 2005). Plasma in hot accretion flows can be collisionless in the sense that the thermalization timescale is longer than the dynamical timescale (Takahara & Kusunose, 1985; Mahadevan & Quataert, 1997; Quataert & Gruzinov, 1999), which motivates us to consider the existence of non-thermal particles and the resulting emissions (Özel et al., 2000; Toma & Takahara, 2012; Kimura et al., 2014, 2015; Chael et al., 2017). Although the spectrum of Sgr A* can be fit by the hot accretion flow models only with thermal electrons (Narayan et al., 1995; Manmoto et al., 1997), the models with non-thermal component match the observations better (Yuan et al., 2003; Ball et al., 2016). Non-thermal protons (hereafter, we call them cosmic rays; CRs) are more commonly expected to exist in hot accretion flows. They interact with thermal protons and photons, leading to gamma-ray and neutrino production. Some predictions have been made using one-dimensional modeling (Mahadevan et al., 1997; Oka & Manmoto, 2003; Niedźwiecki et al., 2013), and observations of nearby Seyferts and LLAGNs by the Fermi satellite give some hints and put limits on the CR production efficiency (Wojaczyński et al., 2015; Wojaczyński & Niedźwiecki, 2017).

Most of the previous studies used a single power-law distribution for the non-thermal components, $dN/d\epsilon\propto\epsilon^{-s}$ with $2\lesssim s\lesssim 4$, where $\epsilon$ is the particle energy. Such a power-law distribution is expected if CRs are produced by first-order Fermi mechanisms, such as the diffusive shock acceleration (Bell, 1978; Blandford & Ostriker, 1978). However, it is unclear whether such a single power-law distribution is achieved, because the accretion flows are unlikely to have a strong shock. Although shocked accretion flows may be formed in hot accretion flows (e.g., Le & Becker, 2005; Becker et al., 2008), we do not observe such structures in the multi-dimensional global hydrodynamic simulations (Yuan & Narayan, 2014). In the accretion flows without shocks, CRs are expected to be produced by magnetic reconnection (e.g., Hoshino, 2012) and/or stochastic acceleration by turbulence (e.g., Lynn et al., 2014). Inside accretion flows, magnetrotational instability (MRI) generates strong turbulence and induces magnetic reconnection (e.g., Balbus & Hawley, 1991, 1998; Sano & Inutsuka, 2001). Recent Particle-In-Cell (PIC) simulations show that when MRI takes place in collisionless plasma, magnetic reconnection produces non-thermal particles (Riquelme et al., 2012; Hoshino, 2013, 2015; Kunz et al., 2016). These non-thermal particles can further be accelerated stochastically through interactions with larger scale eddies. However, current PIC simulations cannot track such a late time phase because of the computational limitation, although recent developments of computational resources and techniques partially enable us to simulate particle acceleration in turbulence (Comisso & Sironi, 2018; Zhdankin et al., 2018). The stochastic particle acceleration by magnetohydrodynamic (MHD) turbulence is often modeled as a diffusion phenomenon in energy space (e.g., Blandford & Eichler, 1987), which has been applied to various astrophysical objects such as galaxy clusters (e.g., Blasi, 2000; Brunetti & Lazarian, 2007; Fujita et al., 2016), gamma-ray bursts (e.g., Asano & Terasawa, 2009; Murase et al., 2012a), radio-lobes of radio galaxies (e.g. Hardcastle et al., 2009; O’Sullivan et al., 2009), and blazars (e.g., Katarzyński et al., 2006; Asano et al., 2014). Engaging this stochastic acceleration model to the hot accretion flow at the Galactic center, we can explain flares of Sgr A* (Liu et al., 2004), TeV gamma rays from the Galactic Center (Liu et al., 2006; Fujita et al., 2015), and perhaps PeV CRs observed at the Earth (Fujita et al., 2017). In addition, Kimura et al. (2015) showed that using the acceleration model, hot accretion flows in LLAGNs can reproduce the high-energy neutrinos detected by IceCube. Note that the model leads to a very hard spectrum, $-1\leq s\leq 0$, compared to the shock acceleration.

In the stochastic acceleration model, the diffusion coefficient in energy space is approximated by a power-law function of energy, $D_{\epsilon}\approx D_{0}(\epsilon/\epsilon_{0})^{q}$. The values of $q$ and $D_{0}$ depend on the power spectrum of the MHD turbulence and interaction processes between CRs and MHD waves (e.g., Cho & Lazarian, 2006). For example, gyro resonant scattering by Alfven modes makes the value of $q$ equal to the slope of the power spectrum of the turbulence (e.g., Dermer et al., 1996; Becker et al., 2006; Stawarz & Petrosian, 2008). The turbulent strength, $(\delta B/B_{0})^{2}$, is related to $D_{0}$, and analytic theories used in the works above assume that the turbulent strength is smaller than unity. However, this condition is likely to be violated in weakly magnetized accretion flows according to MHD simulations (e.g., Stone & Pringle, 2001; McKinney, 2006; Suzuki et al., 2014). Applicability of the analytic models to the strong turbulence has been investigated using test particle simulations, but it is still controversial. The turbulence is usually provided by a superposition of plane waves in the Fourier space (e.g., O’Sullivan et al., 2009; Fatuzzo & Melia, 2014; Teraki et al., 2015), or driven by some algorithms (e.g., Dmitruk et al., 2003; Teaca et al., 2014; Lynn et al., 2014). These studies are useful to investigate features of the stochastic acceleration owing to their controllablity of the turbulence. However, each astronomical object has a different driving mechanism of turbulence, which may lead to a different behaviour of the CR particles (see Roh et al. 2016 for supernova remnants and Porth et al. 2016 for pulsar wind nebulae).

Kimura et al. (2016) performed test-particle simulations in the MRI turbulence using the shearing box approximation (Hawley et al., 1995). However, the shearing box approximation has a few inconsistencies with the hot accretion flows, such as geometrical thickness and non-negligible advection cooling (Narayan & Yi, 1994). More importantly, escape of CRs cannot be implemented in a realistic manner. In this paper, we present results of global simulations, which enables us to investigate behaviours of the high-energy CRs more consistently. We perform MHD simulations to model hot and turbulent accretion flows, and solve orbits of test particles using the snapshot data of the MHD simulations. This paper is organized as follows. First, we describe the global MHD simulations dedicated to the hot accretion flows in Section 2. Then, we show the results of the test-particle simulations in Section 3. We discuss implications and future directions in Section 4, and summarize our results in Section 5.

2 Properties of the MRI turbulence

2.1 Setup for MHD simulations

We use the Athena++ code111https://princetonuniversity.github.io/athena/ to solve the set of the ideal MHD equations (Stone et al., 2008, ; Stone et al. in prep.):

[TABLE]

where $T$ is the time for the MHD calculations, $\rho$ is the density, ${V}$ is the velocity of the MHD fluid, ${B}$ is the magnetic field, $P^{*}=P+B^{2}/(8\pi)$ is the total pressure, $P$ is the gas pressure, $\mathbb{I}$ is the unit tensor, and $\Phi$ is the gravitational potential. The total energy of the fluid is written as

[TABLE]

and we use the equation of state for ideal gas, $P=(\gamma_{s}-1)E_{\rm th}$ ($\gamma_{s}=5/3$ is the specific heat ratio and $E_{\rm th}$ is the thermal energy). We solve the MHD equations in the spherical polar coordinate, $(R,~{}\theta,~{}\phi)$, using the second-order van-Leer integrator, the second-order piecewise linear reconstruction, the Harten-Lax-van Leer Discontinuities (HLLD) approximate Riemann solver (Miyoshi & Kusano, 2005), and the constrained transport scheme. We use the Newtonian gravitational potential, $\Phi=-GM/R$, where $G$ is the gravitational constant and $M$ is the mass of the central black hole (BH). With this potential, we do not have to specify the values of the BH mass and the distance from the BH (see the last paragraph of this subsection about the unit for MHD calculations). Hence, we can use the same MHD data set to multiple particle simulations with various parameters by scaling the physical quantities (see Section 3.1). On the other hand, with the pseudo-Newtonian potential commonly used in the literature (Paczyńsky & Wiita, 1980), we need to give specific values of both the mass of the BH and the radial distance from the BH to simulate the MHD turbulence. This requires more MHD data sets than those with the Newtonian potential, making it impossible to investigate wide parameter space due to a limited computational resource.

The initial condition for the MHD simulations is an equilibrium torus with poloidal magnetic field loops embedded in a non-rotating gas of uniform density without magnetic fields as in Stone & Pringle (2001) (see also Papaloizou & Pringle, 1984; Stone et al., 1999). Within the torus, the pressure is related to the density as $P=K\rho^{\gamma_{s}}$. The density distribution is represented as

[TABLE]

where $n=(\gamma_{s}-1)^{-1}$ is the polytropic index, $R_{c}$ is the radius of the density maximum of the torus, and $d$ is the distortion parameter. We set $\rho=\rho_{\rm env}$ and $P=\rho_{\rm env}/R$ outside the torus. We give the magnetic field through the vector potential: $A_{\phi}=\rho^{2}/\beta_{0}$ and $A_{r}=A_{\theta}=0$ in the torus and \mbox{\boldmath{A}}=0 otherwise. This vector potential produces magnetic fields parallel to the density contours.

Our computational grids extend from $R_{\rm min}=0.1R_{c}$ to $R_{\rm max}=4.0R_{c}$, from $\theta_{\rm min}=\pi/6$ to $\theta_{\rm max}=5\pi/6$, and all the azimuthal angle (i.e. $0\leq\phi\leq 2\pi$). Alfven velocity at the polar region is expected to become very high as the calculation proceeds (e.g., McKinney, 2006), so the polar regions $\theta<\pi/6$ and $\theta>5\pi/6$ are truncated to reduce the computational time. We use the outflow boundary conditions for both inner and outer boundaries of the $R$ and $\theta$ directions. We performed some simulations with other values of $\theta_{\rm max}$ and $\theta_{\rm min}$, and confirmed that they do not affect the features of turbulence discussed in Section 2.2. For our fiducial run, the grid points for the $\theta$ and $\phi$ directions are distributed equally in linear space, while the grids in the radial direction are equally spaced logarithmically.

We perform the simulations with the unit of $R_{c}=1$, $GM=1$, and $\rho_{c}=1$, where $\rho_{c}$ is the maximum density at the torus. We set $\beta_{0}=10^{3}$ and $d=1.5$. For this value of $d$, the innermost radius of the initial torus is located at $R_{\rm trc}\simeq 0.63R_{c}$. We show the results of four runs; they are different only in resolutions and grid spacing as tabulated in Table 1. For runs A, B, and C, the radial grids are equally spaced logarithmically, while for run D, they are equally distributed in the linear space.

2.2 Results of MHD simulations

We are interested in features of the turbulence inside the disc rather than the lower density corona located above the disc. Thus, we calculate mass-weighted average values of physical quantities. We use the snapshot data at $T\Omega_{c}=20\pi$, where $\Omega_{c}$ is the Kepler angular velocity at $R_{c}$, unless otherwise noted. We also analyzed the snapshot data at $T\Omega_{c}=18\pi$ and confirmed that the basic features are the same.

First, we compute volume-averaged quantities of the MRI turbulence:

[TABLE]

The radial integration is performed from $R=R_{\rm min}$ to $R=0.6R_{c}$, where the outer radius is chosen so as to be smaller than the innermost radius of the initial torus. The integration for the other two directions are performed over all the computational region. The grid numbers within the integration region are tabulated in Table 1.

Figure 1 shows the time evolution of Maxwell stresses normalized by the gas pressure, $\alpha_{\rm SS,M}$, and by the magnetic pressure, $\alpha_{\rm mag}$. They are defined as

[TABLE]

respectively. The subscript SS represents Shakura & Sunyaev (1973) in which the $\alpha$ parameter is introduced for the first time222The definition of the $\alpha$ parameter in Shakura & Sunyaev (1973) includes the Reynolds stress, but we do not include it to $\alpha_{\rm SS,M}$ because it is sub-dominant (e.g., Suzuki & Inutsuka, 2014).. We can see that these values converge at $\alpha_{\rm SS,M}\simeq 0.02-0.03$ and $\alpha_{\rm mag}\simeq 0.4-0.5$ for $T\Omega_{c}\gtrsim 50$. This displays that the quasi-steady state is realized. We compute the mass accretion rates and find that they are roughly constant in both time and radius, which also indicates the realization of the quasi-steady state. $\alpha_{\rm mag}$ is used as an indicator of the numerical convergence of the MRI turbulence (Hawley et al., 2011; Hawley et al., 2013). We tabulate the values of $\alpha_{\rm mag}$ at the end of the calculations. Sufficiently high-resolution simulations give $\alpha_{\rm mag}\simeq 0.45$, which is seen in runs A, B, and D. The ratio of the radial magnetic energy to the azimuthal one, $\langle B_{R}^{2}\rangle_{\mathcal{V}}/\langle B_{\phi}^{2}\rangle_{\mathcal{V}}$ is also useful in diagnosis (Hawley et al., 2011; Hawley et al., 2013). The values of the ratio converge to 0.15–0.16 for runs A, B and D as tabulated in Table 1. Hence, the grid numbers for our simulations are high enough to follow the features of the MRI turbulence except for run C.

The inverse of plasma beta, $B^{2}/(8\pi P)$, and each component of magnetic and electric fields are also tabulated in Table 1. Although the MRI amplifies the magnetic field, the gas pressure is still stronger than the magnetic pressure. In the accretion flow, $B_{\phi}$ dominates over the other two components because the shear motion stretches the magnetic field. The electric field is computed by the ideal MHD condition,

[TABLE]

Since $B_{\phi}$ and $V_{\phi}$ are dominant, $E_{r}$ and $E_{\theta}$ are stronger than $E_{\phi}$. Note that the electric field in the accretion flow is much weaker than the magnetic field, since $V/c<1$.

The radial profiles of the velocities are estimated to be

[TABLE]

The fluid velocity strongly affects the orbits of the test particles, since the drift velocities of the test particles are the same as the fluid velocity perpendicular to the magnetic field. Each component of the fluid velocity is divided into two parts, $V_{i}=V_{{\rm bul},i}+V_{{\rm tur},i}$. For $V_{R}$ and $V_{\phi}$, assuming that $\langle V_{{\rm tur},i}\rangle_{\mathcal{S}}=0$, we obtain $V_{{\rm bul},i}=\langle V_{i}\rangle_{\mathcal{S}}$. For $V_{\theta}$, we define $V_{{\rm bul},\theta}=\langle V_{\theta}{\rm sgn}(\cos\theta)\rangle_{\mathcal{S}}$. Then, the turbulent component for each direction is estimated to be $V_{{\rm tur},i}=\sqrt{\langle V_{i}^{2}\rangle_{\mathcal{S}}-V_{{\rm bul},i}^{2}}$. We plot $V_{{\rm bul},i}$ and $V_{{\rm tur},i}$ in the upper panel of Figure 2. The turbulent components dominate over the bulk components for $V_{R}$ and $V_{\theta}$. On the other hand, $V_{{\rm bul},\phi}$ is much higher than its turbulent component. The turbulent components are comparable each other, but $V_{{\rm tur},R}$ is slightly higher than the other components.

The lower panel of Figure 2 shows comparison of the turbulent velocity ($V_{{\rm tur},R}$), rotation velocity ($V_{{\rm bul},\phi}$), sound velocity ($C_{s}=\sqrt{\gamma_{s}P/\rho}$), and Alfven velocity ($V_{A}=B/\sqrt{4\pi\rho}$). We see that the rotation velocity is the fastest of the four and its value is similar to the Keplerian velocity ($V_{K}=\sqrt{GM/R}$) shown as the thin line. The sound velocity is the second fastest, and Alfven velocity follows it. This means that the accretion flow consists of high-$\beta$ plasma and that the turbulence is sub-sonic and sub-Alfvenic. These results are consistent with the previous simulations (e.g., Stone & Pringle, 2001; Machida & Matsumoto, 2003).

Next, we briefly discuss two-dimensional structures. The left and middle panels of Figure 3 show the colormaps of the density and the magnetic energy on the $\phi=0$ plane, respectively. We can see from the left panel that the accretion flow is geometrically thick with the aspect ratio $H/R\sim 0.5$, and the density at midplane is much higher than that around the polar boundaries. The middle panel shows that the magnetic field is strongly turbulent not only in the disc region ($|\cos\theta|\lesssim 0.45$) but also in the corona region ($|\cos\theta|\gtrsim 0.45$). The right panel represents contours of the azimuthally averaged $V_{\phi}$:

[TABLE]

We can see that the iso-$V_{\phi}$ surface depends on the spherical radius, $R$, rather than the cylindrical radius, $R\sin\theta$. This allows us to use $V_{\rm{bul},\phi}$ as the background velocity for analyses of the test-particle simulations in Section 3.2.

Figure 4 plots the colormaps of the density (upper) and the magnetic energy (lower) on the equatorial plane. The magnetic fields are frozen in the differentially rotating fluid elements that fall to the BH. This creates the spiral structure as seen in the figure. We can also see the fluctuation of the density is much smaller than that of the magnetic field energy density. This implies that the fast modes are a sub-dominant component in the MRI turbulence.

To clarify the importance of the modes of the MHD waves (fast, slow, and Alfven), we evaluate the Pearson correlation coefficients between the fluctuations of the density, $\delta\rho(R,~{}\theta,~{}\phi)=\rho-\langle\rho\rangle_{\mathcal{L}}$, and the magnetic energy, $\delta B^{2}(R,~{}\theta,~{}\phi)=B^{2}-\langle B^{2}\rangle_{\mathcal{L}}$. According to the linear MHD wave theory, the fast mode has a positive correlation, the slow mode has a negative correlation, and the Alfven mode has no correlation. We evaluate the correlation coefficients as a function of $R$ and $\theta$, and average over them with weights associated with the area in the meridional plane. The resulting coefficients indicate that the density and magnetic energy are weakly anti-correlated: the value of the coefficient is $-0.22$ in the disc region ($|\cos\theta\lesssim 0.45|$) for run A. The lower resolution runs have higher coefficients, i.e., the anti-correlations are weaker, but no run has a positive correlation. Therefore, the fast modes do not play an important role in this system. This result is natural in the sub-Alfvenic and sub-sonic turbulence.

Finally, we discuss the azimuthal power spectra of the turbulence (cf. Sorathia et al. 2012; Suzuki & Inutsuka 2014; see Parkin & Bicknell 2013 for three-dimensional power spectra). We take the Fourier transformation in the azimuthal direction,

[TABLE]

where $m=k_{\phi}R$ ($k_{\phi}$ is the wavenumber in the $\phi$ direction). Then, we take the average of the power spectrum over the disc region:

[TABLE]

where the integration region is set to be $0.1R_{c}\leq R\leq 0.6R_{c}$ and $|\cos\theta|\leq 0.45$. We plot the power spectra, $m\mathcal{P}_{m}$, for the magnetic field in Figure 5. We can see that all the data sets have similar values for a larger scale of $m\lesssim 10$. The spectra for $B_{r}$ and $B_{\theta}$ are $m\mathcal{P}_{m}\propto m$, while those for $B_{\phi}$ are roughly $m\mathcal{P}_{m}\propto m^{1/2}$. For a smaller scale of $m\gtrsim N_{\phi}/10$, the spectra decrease with $m$ very rapidly for all the data sets because of the numerical dissipation. The power spectra peak at intermediate scale of $m\sim 10-20$, depending on the resolution and component. These features are consistent with the previous calculations (Sorathia et al., 2012; Suzuki & Inutsuka, 2014).

The fastest growing mode of the MRI is approximated to be $L_{\rm MRI}\sim 2\pi V_{A}/\Omega$, where $\Omega$ is the angular velocity. Saturation of MRI turbulence is expected to be controlled either by the large-scale magnetic reconnection (Sano & Inutsuka, 2001; Sano et al., 2004) or by the growth of the parasitic instabilities of Kelvin-Helmholtz modes (Goodman & Xu, 1994; Pessah, 2010). These phenomena occur inside the disc, where the largest scale is the scale height, $H\approx C_{s}/\Omega$. Hence, the characteristic scale of the saturated MRI turbulence should be the smaller one of the two, $L_{\rm tur}\approx{\rm min}(L_{\rm MRI},~{}H)$. From Figure 2, we roughly see $V_{A}\sim V_{\phi,\rm bulk}/7$ and $C_{s}\sim V_{\phi,\rm bulk}/2$, leading to $L_{\rm MRI}\approx 2\pi R/7>H\approx R/2$. Hence, $L_{\rm tur}=H\approx R/2$. This scale corresponds to $m\sim 13$, which is consistent with the peaks of the power spectra.

For the intermediate scale, we narrowly see that the spectra gradually decrease with $m$. Theoretically, fully developed Alfven turbulence results in $P_{k}\propto k_{\perp}^{-5/3}$ and $P_{k}\propto k_{\parallel}^{-2}$, where $P_{k}$ is the power spectrum, $k_{\perp}$ and $k_{\parallel}$ are the perpendicular and parallel wave numbers to the background magnetic field, respectively (Goldreich & Sridhar, 1995). Such an anisotropic cascade takes place with respect to the local magnetic field. In strong turbulence where the large-scale magnetic field is significantly tilted, the direction of the local magnetic field is not aligned. Then, the global Fourier analysis would smear out the local anisotropy, resulting in $P_{k}\propto k^{-5/3}$ in all the directions (Cho & Vishniac, 2000). However, we cannot clearly see the power-law shape in the power spectra of our simulations, due to the insufficient dynamic range. Simulations with a higher resolution and a higher-order reconstruction scheme are necessary to determine the power-law index in the inertial range. To observe the anisotropic feature, even more dedicated analyses reconstructing coordinates based on the local magnetic field will also be required. Note that the shape of these power spectra is independent of the integration range of $R$, because the turbulence is generated by the same mechanism at all the radii.

3 Behaviours of High-Energy Particles

3.1 Setup for particle simulations

We calculate orbits of relativistic particles to investigate behaviour of high-energy particles in the accretion flows. We ignore CR injection mechanisms because they are related to small-scale plasma processes. They should be investigated by other methods, such as PIC simulations (Hoshino, 2015; Kunz et al., 2016), which is beyond the scope of this paper.

We solve the relativistic equation of motion for each CR particle:

[TABLE]

where $t$ is the time for particle calculation, $c$ is the speed of light, and \mbox{\boldmath{p}}=\gamma m_{p}\mbox{\boldmath{v}}, ${v}$, $e$, $m_{p}$, and $\gamma=\sqrt{1-(v/c)^{2}}$ are the momentum, velocity, charge, mass, and Lorentz factor of the CR particle, respectively. Here, we neglect the gravity acting on the CR particle, since it is typically weaker than the electromagnetic force by more than ten orders of magnitude. This equation is integrated using the Boris method (e.g., Birdsall & Langdon, 1991), which is often used in PIC simulations. In the particle simulations, we use $m_{p}$ and $e$ for protons, but we can scale our simulation results to the heavy nuclei using the rigidity $\mathcal{R}=\epsilon/Z$.

The snapshot data of the MHD simulations shown in Section 2.2 are used to obtain ${E}$ and ${B}$. Since the MHD data contain the values of ${V}$ and ${B}$ only at the discrete grid points, we first interpolate ${B}$ and ${V}$ at the position of the particle using quadratic functions333Although we use the quadratic functions for the interpolation, the results are very similar if we use the linear interpolation.. Then, we compute ${E}$ through Equation (10) using the interpolated ${B}$ and ${V}$. This procedure guarantees \mbox{\boldmath{E}}\cdot\mbox{\boldmath{B}}=0, so artificial acceleration due to the interpolation is avoided.

We initially distribute particles on a ring of $R=R_{\rm ini}$ and $\theta=\pi/2$. The energy distribution of the initial particles is monoenergetic and isotropic in the fluid frame (see Section 3.2 for the definition of the fluid frame). The initial radius is fixed at $R_{\rm ini}=0.3R_{c}$ for simplicity. We performed the simulations with $R_{\rm ini}=0.2R_{c}$, and checked that the results are almost unchanged. The initial energy of the particle, $\epsilon_{\rm ini}$, is given so that the Larmor radius of the particle is equal to $\lambda_{\rm ini}$ times the grid scale: $r_{L}=\epsilon_{\rm ini}/(ecB_{\rm ave})=\lambda_{\rm ini}\Delta x_{\rm ini}$, where $\Delta x_{\rm ini}={\rm min}(\Delta R_{\rm ini},~{}R_{\rm ini}\Delta\theta,~{}R_{\rm ini}\Delta\phi)$ is the grid scale at the initial ring. The timestep of the particle calculation is determined by $\Delta t={\rm min}(\Delta t_{L},\Delta t_{x})$, where $\Delta t_{L}=C_{\rm safe}t_{L,\rm min}=2\pi C_{\rm safe}\epsilon_{\rm ini}/(ecB_{\rm max})$ and $\Delta t_{x}=C_{\rm safe}\Delta x_{\rm min}/c$. Here, $B_{\rm max}$ is the maximum value of the magnetic field, $\Delta x_{\rm min}$ is the minimum length between the grids in the computational region, and $C_{\rm safe}$ represents the safety factor that determines the timestep. We set $C_{\rm safe}=0.01$. We performed some simulations with $C_{\rm safe}=0.001$, and confirmed that the results are unchanged by the values of $C_{\rm safe}$. As a fiducial value, we set $\lambda_{\rm ini}=4$. With a smaller value of $\lambda_{\rm ini}$, we cannot trace the resonant scattering process, while the particles escape from the computational region too quickly with a higher value of $\lambda_{\rm ini}$.

The computational region for the particle simulations is the same with the MHD simulations except for the outer boundary in the $R$ direction. Since the dynamical structures of the outer parts of the MHD simulations are affected by the initial conditions, we set the outer boundary of the particle simulations to $R_{\rm esc}=0.6R_{c}$. The particles that go beyond the computational region are removed from the simulation, and we stop the calculation when half of the particles escape from the computational region.

We solve the equations of motion for $N_{p}=2^{14}=16384$ particles using the MHD data sets shown in the previous section. To solve the equation of motion, we need to convert the units used in the MHD calculations to those of our interest. The units of the mass, length, and time for the MHD calculations are written as $\mathcal{L}_{u}=R_{c}$, $\mathcal{M}_{u}=\rho_{c}R_{c}^{3}$, and $\mathcal{T}_{u}=\sqrt{R_{c}^{3}/(GM)}$, respectively. For our particle simulations, we rescale these units as

[TABLE]

where $R_{s}=2GM/c^{2}$ is the Schwarzschild radius, $\dot{M}_{\rm Edd}=L_{\rm Edd}/c^{2}$ is the Eddington mass accretion rate ($L_{\rm Edd}$ is the Eddington luminosity), and $\chi$ and $\eta$ are the scaling factors of the length and the mass, respectively. The relation between $\eta$ and density is $\rho_{c}=\eta\dot{M}_{\rm Edd}\mathcal{T}_{u}/\mathcal{L}_{u}^{3}$, so a higher $\eta$ leads to a higher density.

We choose the reference parameter set for the particle simulations so as to be consistent with our assumptions: hot accretion flows in LLAGNs with Newtonian gravity. In our MHD simulations, mass accretion rate is written as $\dot{M}\sim\dot{m}_{\rm sim}\mathcal{M}_{u}\mathcal{T}_{u}^{-1}$, where $\dot{m}_{\rm sim}\sim 10^{-3}-10^{-2}$ is the resulting mass accretion rate in the MHD simulations. Then, the rescaled mass accretion rate is represented as $\dot{M}=\eta\dot{m}_{\rm sim}\dot{M}_{\rm Edd}$. For $\eta\lesssim 10$, this mass accretion rate is in the hot accretion flow regime, i.e., $\dot{M}\lesssim 0.1\dot{M}_{\rm Edd}$ (Narayan & Yi, 1995; Xie & Yuan, 2012). The scale factor for the length, $\chi$, should be large enough to be consistent with the Newtonian gravity. For $\chi\geq 20$, the initial radius, $R_{\rm ini}=0.3R_{c}$, is larger than $6R_{s}=2R_{\rm ISCO}$, where $R_{\rm ISCO}=3R_{s}$ is the innermost circular stable orbit (ISCO) for the Schwartzchild BH. Based on the considerations above, we set the reference parameters to $\chi=50$, $\eta=1$, and $M=10^{8}M_{\odot}$, which corresponds to typical low-luminosity AGNs, such as Seyferts or low-ionization nuclear emission-line regions (LINERs). This parameter set leads to

[TABLE]

where we use the notation $Q_{x}=10^{x}$ (unit for $M$ is M⊙). The speed of light is $c\simeq 10\chi_{1.7}^{1/2}\mathcal{L}_{u}\mathcal{T}_{u}^{-1}$ in this unit system. We use the MHD data set of run A with $T\Omega_{c}=20\pi$ unless otherwise noted. The Larmor radius and timescale are $r_{L}\simeq 1.0\times 10^{13}M_{8}\chi_{1.7}\lambda_{\rm ini,0.6}$ cm and $t_{L}\simeq 2.1\times 10^{3}M_{8}\chi_{1.7}\lambda_{\rm ini,0.6}$ s, respectively.

3.2 Results of particle simulations

3.2.1 Orbits and momentum distribution

The upper and lower panels of Figure 6 show orbits of the test particles projected in the $R-\theta$ and $R-\phi$ planes, respectively. The particles spread in all the directions, but majority of the particles move outward in the $R$ direction rather than fall onto the BH or escape to the vertical direction. The particles are likely to migrate to the direction at which the magnetic field is weak, partially due to the magnetic mirror force. The strong magnetic fields at the inner region prevent the particles from falling to the BH. Also, the magnetic fields at the high latitudes are not week, compare to those at the outer region (see Figure 3). At the end of the simulation, 68% of all the escaping particles go out through the radial boundary with $|\cos\theta|\leq 0.5$ for the case with $\lambda_{\rm ini}=4$. This fraction is higher for lower $\lambda_{\rm ini}$, and vice versa for higher $\lambda_{\rm ini}$. Higher energy particles can cross the magnetic field more easily, which may enhance the vertical diffusion.

From the lower panel, we find that the particles mostly travel along the $\phi$ direction, because the magnetic field is also directed to the $\phi$ direction. Interestingly, the outward-going particles tend to move the opposite direction to the background fluid. This arises from the magnetic field configuration. The accretion flow creates the spiral-shape magnetic fields as seen in Figure 4. When the CR particles stream along the field line outward, they counter-rotate with respect to the accretion flow.

The upper panel of Figure 7 shows the momentum distribution in each direction, $dN/dp_{i}$ ($i=r,~{}\theta,~{}\phi$) measured in the lab frame at time $t=10t_{L}$, where $t_{L}=2\pi\epsilon_{\rm ini}/(ecB_{\rm ini})$ and $B_{\rm ini}=\sqrt{\langle B^{2}\rangle_{\mathcal{S}}}$ at $R=R_{\rm ini}$. The momentum distribution is anisotropic: there is a bulk rotational motion. This is because the background fluid motion creates the electric field that induces \mbox{\boldmath{E}}\times\mbox{\boldmath{B}} drift.

We compute the momentum distribution in the fluid frame by performing Lorentz transformation of the particle momenta. Since $V_{\phi}$ dominates over the other component, we approximate the background velocity to be

[TABLE]

where \mbox{\boldmath{e}}_{\phi} is the unit vector of the $\phi$ direction and $V_{{\rm bul},\phi}$ is independent of $\theta$. The bottom panel shows the momentum distribution in the fluid frame, where we can see no bulk rotational motion. In the following sections, we use the energy distribution in the fluid frame. Note that the particle distribution is slightly anisotropic: the particles tend to have positive $p_{R}$ and negative $p_{\phi}$. This is because the particles tend to move radially outward along the spiral magnetic field, as discussed above. This anisotropy becomes stronger in later time and for higher energy particles (see Section 3.2.3). Since this anisotropy appears in the particle simulations with all the MHD data sets, the grid spacing and resolutions are not the cause of the anisotropy.

3.2.2 Diffusion in energy space

We examine evolution of the energy distribution function in the fluid frame. The time evolution of the energy distribution for $\lambda_{\rm ini}=4$ is shown in Figure 8. We can see that the width of the energy distribution increases with time. This motivates us to consider the diffusion equation in the energy space.

In general, the transport equation, including the diffusion and advection terms in both configuration and momentum spaces, describes the evolution of the distribution function for the particles with isotropic distribution in the fluid rest frame (e.g. Skilling, 1975; Strong et al., 2007). When the terms for configuration space and the advection term in momentum space are negligible, the transport equation may be simplified to the diffusion equation only in momentum space (e.g., Stawarz & Petrosian, 2008):

[TABLE]

Since the anisotropy in our system is not very strong, we apply this equation to our system. We focus on the ultrarelativistic regime, so the particle energy is approximated to be $\epsilon\approx pc$. Using the differential number density, $N_{\epsilon}=N_{p}/c=4\pi p^{2}f/c$, we can write the evolution of $N_{\epsilon}$ by the advection-diffusion equation in energy space:

[TABLE]

In our simulation, most of the particles are confined in a narrow energy range of $\epsilon\sim\epsilon_{\rm ini}$. Thus, we approximate that $D_{\epsilon}\approx D_{\epsilon_{\rm ini}}$ and $2D_{\epsilon}/\epsilon\approx v_{\epsilon_{\rm ini}}$ are constant. Then, after some algebra, we obtain the time evolution of the mean and variance of the energy distribution:

[TABLE]

[TABLE]

where $N_{\rm actv}$ is the number of the particles confined in the computational region (see Appendix A for derivation).

Figure 9 shows the time evolution of $\sigma^{2}_{\epsilon}$ for $\lambda_{\rm ini}=1$, 2, 4, and 8 from $t=0.1t_{L}$ to the end of the simulation. Initially, $\sigma^{2}_{\epsilon}$ rapidly increases with time for all the models. Due to the turbulent velocity component whose amplitude is about 10 % of the background velocity, the particles are not exactly in the rest frame of the fluid elements, which causes drift motions. This leads to the initial jump of $\sigma^{2}_{\epsilon}$. For $\lambda_{\rm ini}=1$, 2, and 4, $\sigma^{2}_{\epsilon}$ becomes almost constant at $\sigma^{2}_{\epsilon}\sim 10^{-4}\epsilon_{\rm ini}^{2}$ after $t\gtrsim 0.3t_{L}$, because the particles start to move with the local drift velocity. In the late time, the particles start interacting with the turbulence, and $\sigma^{2}_{\epsilon}\propto t$ is realized. For $\lambda_{\rm ini}=8$, $\sigma^{2}_{\epsilon}$ continues to increase and approaches $\sigma^{2}_{\epsilon}\propto t$, because the gyration period is comparable to the interaction timescale with the turbulence. From the results, we find that the interaction timescale can be expressed as $t_{\rm int}\sim 1\times 10^{3}-2\times 10^{3}M_{8}\chi_{1.7}$ s, including the parameter dependence. This is about a factor of 4 shorter than the crossing time of the turbulent length, $L_{\rm tur}/c\sim 7\times 10^{3}M_{8}\chi_{1.7}$ s, i.e., we can write $t_{\rm int}\sim L_{\rm tur}/(4c)$.

This result, $\sigma^{2}_{\epsilon}\propto t$ for $t>t_{\rm int}$, indicates that the particle acceleration in the MRI turbulence occurs through the diffusion in energy space. In the sub-sonic turbulence including the MRI turbulence, the slow modes are expected to play an important role in particle scattering (Lynn et al., 2014). We can consider two mechanisms that change the particle energy in such turbulence: the Fermi-type B mechanism (FTB; see e.g. Lynn et al., 2012) and the transit-time damping (TTD). In FTB, the particles stream along a curved magnetic field that has a velocity. Then, the particles gain or lose energy at the fluid frame after the magnetic field sufficiently change the direction (see Figure 1 of Lynn et al. (2012)). The mean velocity of the magnetic field is expected to be $V_{R,\rm tur}$ in our MHD simulation. Then, the energy change per “collision” is approximated to be $\Delta\epsilon\sim\epsilon V_{R,\rm tur}/c$. Using the interaction time with the turbulence, $t_{\rm int}\approx L_{\rm tur}/(4c)$, the diffusion coefficient in energy space can be estimated to be (e.g., Blandford & Eichler, 1987)

[TABLE]

In the last equation, we write down the parameter dependence using $\epsilon\approx\epsilon_{\rm ini}\propto\lambda_{\rm ini}M^{1/2}\chi^{-1/4}\eta^{1/2}$. We can use the relation because the energy of the particles does not change very much in each run of the particle simulations.

TTD requires the resonant condition: $v_{\rm pha}\simeq v_{\parallel}$, where $v_{\rm pha}$ is the phase velocity of the slow mode and $v_{\parallel}$ is the particle velocity parallel to the magnetic field. For the relativistic particles in weak sub-sonic turbulence, the condition for TTD cannot be satisfied, because $v_{\parallel}\sim c$ is always much faster than $v_{\rm pha}\sim V_{A}$. However, in strong turbulence, the relativistic particles can interact with the slow mode because non-linear effects broaden the energy range of the resonant particles (Yan & Lazarian, 2008; Lynn et al., 2014). If TTD is effective, the mean energy change per collision is typically $\epsilon V_{A}/c$. Then, the diffusion coefficient in energy space can be estimated to be

[TABLE]

The parameter dependence of $D_{\epsilon,\rm TTD}$ is the same as that of $D_{\epsilon,\rm FTB}$, while the normalization of $D_{\epsilon,\rm TTD}$ is higher than that of $D_{\epsilon,\rm FTB}$.

We calculate $D_{\epsilon}$ with various values of $\lambda_{\rm ini}=$ (0.5, 1, 2, 4, 8, 16), $M=$ ($10^{7}~{}M_{\odot},~{}10^{8}~{}M_{\odot}$), $\chi=$ (30, 50, 100, 200), and $\eta=$ (1, 10), and show the resulting $D_{\epsilon}$ in Figure 10. We combine the simulation results with various $\epsilon_{\rm ini}$ to discuss the energy dependence of $D_{\epsilon}$. For the calculations with $\epsilon_{\rm ini}\gtrsim 10^{3}$ PeV, the particles escape from the computational region before the condition $\sigma^{2}_{\epsilon}\propto t$ is realized, so we only plot the results with $\epsilon_{\rm ini}<10^{3}$ PeV. The parameter dependence of $D_{\epsilon}$ is consistent with both of the simple estimates above: $D_{\epsilon}\propto\epsilon^{2}$ in the upper panel and $D_{\epsilon}\propto\chi^{-5/2}\eta$ for $\lambda_{\rm ini}=4$ in the lower panel. The normalization of the simple estimates are consistent with the simulation results within a factor of 3, while $D_{\epsilon,\rm FTB}$ matches better than $D_{\epsilon,\rm TTD}$. For the rest of the paper, we use $D_{\epsilon,\rm FTB}$ as a diffusion coefficient in energy space. The acceleration time is estimated to be

[TABLE]

This acceleration time is independent of energy. Note that the fast modes have little influence on particle scattering in our simulation because they do not have enough power as discussed in Section 2.2.

The values of $D_{\epsilon}$ can be estimated in two ways using either $\mu_{\epsilon}$ or $\sigma_{\epsilon}$, and these can be different when the particle distribution is anisotropic. We evaluate the time evolution of $\mu_{\epsilon}$ via Equation (25), and confirm that the two methods are consistent with each other within a factor of 3 for $\lambda_{\rm ini}\lesssim 8$. So far, we have assumed that $D_{\epsilon}$ is constant, but our results indicate that $D_{\epsilon}\propto\epsilon^{2}$ is more realistic. For the case with $D_{\epsilon}\propto\epsilon^{2}$, we can derive the time evolution of $\mu_{\epsilon}$ and $\sigma_{\epsilon}^{2}$ without assuming that $D_{\epsilon}$ is constant. Then, as shown in Appendix, the evolution of $\sigma_{\epsilon}^{2}$ is unchanged, while the increasing rate of $\mu_{\epsilon}$ is twice higher than that given by Equation (25):

[TABLE]

We also estimate the values of $D_{\epsilon}$ based on the above Equation, and the results agree with those obtained by $\sigma_{\epsilon}^{2}$ within a factor of 2 for $\lambda_{\rm ini}\leq 8$ as shown in Appendix, implying the improvement compared to those based on Equation (25). For the models with $\lambda_{\rm ini}\gtrsim 8$, the anisotropy is large enough to affect the momentum diffusion, and the agreement becomes worse. However, this does not affect the discussion on the maximum energy in Section 4, because in reality, high-energy particles escape from the system before they attain the energy corresponding to $\lambda_{\rm ini}\gtrsim 8$.

3.2.3 Behaviour in configuration space

First, we discuss the displacement in $R$ direction, which is directly related to the escape process. We estimate time evolutions of the mean and the variance of the radial displacement:

[TABLE]

where $\delta R_{j}=R_{j}-R_{\rm ini}$ is the radial displacement of each particle and $N_{\rm actv}$ is the number of the confined particles. Summation is performed over the particles confined in the computational region. $\mu_{\delta R}$ represents the bulk motion of CR particles, while $\sigma^{2}_{\delta R}$ expresses the diffusive motion.

We show $\mu^{2}_{\delta R}$ and $\sigma^{2}_{\delta R}$ for the cases with various energies ($\lambda_{\rm ini}=2$, 4, and 8) in Figure 11. In the early phase, the diffusive motion is more efficient than the bulk CR motion. From the figure, we see that the stochastic behaviour of CRs in configuration space cannot be described as a usual diffusion. If the particles obey the usual diffusion, $\sigma^{2}_{\delta R}\propto t^{2}$ at the beginning, and $\sigma^{2}_{\delta R}\propto t$ after scattering timescale (e.g., Casse et al., 2002; Cohet & Marcowith, 2016). In our simulation, $\sigma^{2}_{\delta R}$ initially increases with $t^{2}$. After about a half of the Larmor timescale, we see a transition to $\sigma^{2}_{\delta R}\propto t^{1.4}$, and finally, $\sigma^{2}_{\delta R}$ becomes flat at the time when the escape fraction becomes non-negligible. We analyze the data with various computational regions, and find that the behaviour is essentially the same; $\sigma^{2}_{\delta R}$ rapidly increases initially, and it is flattened when the escape becomes effective, as seen in Figure 12. This trend is similar for the cases with different parameter sets and other sets of the snapshot data of the MHD simulations. Thus, the radial variance is approximately written as

[TABLE]

where $\xi\simeq 1.2-1.4$ and $\zeta\simeq 2-3$. Since $1<\xi<2$, the particle behaviour in configuration space in the MRI turbulence is superdiffusive. Note that $\sigma^{2}_{\delta R}$ starts to increase around $t\sim t_{L}/\zeta$, while $\sigma_{\epsilon}$ discussed in the previous subsection begins to grow at $t\sim t_{\rm int}$.

In later time, $\mu_{\delta R}^{2}$ rapidly increases with time, and becomes the dominant process for a higher energy run. This bulk motion originates from the anisotropic behaviour discussed in Section 3.2.1. Using the distribution function, we can estimate the radial component of the bulk velocity of the CR particles:

[TABLE]

For all the runs, $\beta_{\rm CR,bul}$ is not constant at the beginning. $\beta_{\rm CR,bul}$ starts to increase rapidly around $t\sim t_{\rm int}$, and for $t\gtrsim 5t_{\rm int}$, $\beta_{\rm CR,bul}$ becomes almost constant with time. We estimate this constant values of $\beta_{\rm CR,bul}$ for various cases, and plot its energy dependence in the upper panel of Figure 13. The higher energy runs have higher values of $\beta_{\rm CR,bul}$. The bottom panel shows the dependence on the other parameters for $\lambda_{\rm ini}=4$, where we see that $\beta_{\rm CR,bul}$ is independent of $M$, $\chi$, and $\eta$. According to our simulation results, the bulk radial velocity is represented as

[TABLE]

where $\epsilon_{\rm PeV}=\epsilon/\rm PeV$. This power-law representation fits the simulation data well as seen in Figure 13. This bulk CR motion is much faster than the background MHD fluid, $V_{R}/c\sim 2\times 10^{-3}\chi_{1.7}^{-1/2}$. Also, the bulk CR motion is outward, while the background fluid moves inward. Note that this energy dependence of the bulk CR motion appears due to our mono-energetic simulations. In reality, the CRs have an energy distribution, and convolution of CRs for all the energies will provide an energy-independent CR bulk velocity.

We estimate the escape time using $\sigma^{2}_{\delta R}$ and $\beta_{\rm CR,bul}$. If the bulk CR motion is effective, the escape time can be estimated to be

[TABLE]

On the other hand, for the diffusion-dominant cases, the particle starts to escape when $\sigma^{2}_{\delta R}\gtrsim(R_{\rm esc}-R_{\rm ini})^{2}$. Then, we can estimate the escape time to be

[TABLE]

We expect $t_{\rm esc,d}\propto\epsilon^{-0.43}-\epsilon^{-0.82}$ with $\xi\simeq 1.2-1.4$. For the demonstration purpose, we use $\xi=1.3$ and $\zeta=3$ for the rest of the paper, which leads to $t_{\rm esc,d}\propto\epsilon^{-0.54}$ for fixed values of $M$, $\chi$, and $\eta$ and $t_{\rm esc,d}\propto M\chi$ for a fixed value of $\lambda_{\rm ini}$. Note that the dependences on $M$, $\chi$, and $\eta$ for a fixed value of $\epsilon$ differ from those for a fixed $\lambda_{\rm ini}$.

We calculate a typical escape time, $t_{\rm esc}$, at which half of the particles escape from the system, and show it as a function of the particle energy in the upper panel of Figure 14444The definition of $t_{\rm esc}$ is arbitrary. If we define $t_{\rm esc}$ as the time when the number of the particles inside the computational region becomes $e^{-1}$ of the initial number, $t_{\rm esc}$ can be longer than those in Figure 14. The difference between these two definition is within a factor of 1.4. The higher-energy particles escape from the system faster than the lower ones. The dotted and dashed lines show the escape time given by Equation (36) and (37), respectively, which roughly fit the simulation results for $\epsilon\lesssim 10^{3}$ PeV. Interestingly, these two estimates give similar values. The slope becomes steeper for higher-energy particles at $\epsilon\gtrsim 10^{3}$ PeV, because these particles escape before interacting with the turbulent field. In the lower panel of Figure 14, we calculate $t_{\rm esc}$ with various parameter sets, and find that the parameter dependence of $t_{\rm esc}$ is $t_{\rm esc}\propto M\chi$ for a fixed value of $\lambda_{\rm ini}$, which is also consistent with both of the estimates. Note that this escape time is longer than the light crossing time, $0.3R_{c}/c\sim 1.5\times 10^{4}M_{8}\chi_{1.7}$ s as long as $\epsilon\lesssim 10^{3}$ PeV. Thus, the escape process is not ballistic. Note also that our results indicate that the advection by MHD fluid motion, either turbulent motion or outflow, is not a dominant process of the particle transport. The advection time by MHD fluid is represented as $t_{\rm MHD}\approx R/V_{\rm MHD}$, where $V_{\rm MHD}$ is the turbulent velocity or the outflow velocity. The parameter dependence of the advection time by the MHD fluid is $t_{\rm MHD}\propto M\chi^{3/2}$, which is inconsistent with our results.

Next, we discuss the anisotropic behaviour in configuration space. The variances of the polar and azimuthal displacements are estimated to be

[TABLE]

where $\delta\theta_{j}=\theta_{j}-\pi/2$ and $\delta\phi_{j}=\phi_{j}-\phi_{j0}$ are the displacements in $\theta$ and $\phi$ directions, respectively ($\phi_{j0}$ represents the initial azimuthal position), $\mu_{\delta\theta}$ and $\mu_{\delta\phi}$ are the means of the displacements. Figure 15 plots $\sigma^{2}_{\delta R}$, $R_{\rm ini}^{2}\sigma^{2}_{\delta\theta}$, and $R_{\rm ini}^{2}\sigma^{2}_{\delta\phi}$ for $\lambda_{\rm ini}=4$. Initially, all lines increase with $t^{2}$, and later, they changes the increasing rates to $t^{1.3}-t^{1.4}$. We can regard the break time as the mean free time for each direction. The mean free time for the azimuthal direction is around $10^{4}$ s, which is longer than those for the radial and polar directions ($\sim t_{L}/\zeta\simeq 7\times 10^{2}$ s). This means that the particles stream to $\phi$ direction without strongly being scattered for a longer time. Since the magnetic field is directed to the azimuthal direction, this result also indicates that particles tend to stream parallel to the magnetic field. The mean free time in $\phi$ direction is the same order as the time when $\beta_{\rm CR,bul}$ becomes constant, $\sim 5t_{\rm int}$. Since the particles escape to radial direction, $\sigma^{2}_{\delta\theta}$ and $\sigma^{2}_{\delta\phi}$ continue to increase even after the non-negligible fraction of the particles escape from the system. We also estimate $\mu_{\delta\theta}$ and $\mu_{\delta\phi}$, and find that the bulk motions in $\theta$ and $\phi$ directions are sub-dominant for the parameter space we investigated here.

To discuss whether scattering between waves and particles occurs, we often use the first adiabatic invariance, $\mu_{\rm ad}=p_{\perp,*}^{2}/(2mB_{*})$, where $p_{\perp,*}$ is the momentum perpendicular to the magnetic field and the subscript $*$ indicates the values at the fluid frame. In our simulation, since the magnetic field is turbulent in the scale of the Larmor radius, the adiabatic invariance is not conserved. Hence, we cannot use $\mu_{\rm ad}$ for measuring the mean free time.

4 Discussion

4.1 Maximum energy of CRs

The escape time is given by Equations (36) or (37), and the acceleration time is written in Equation (29). We here discuss implications of these results, and estimate the maximum achievable energy as an example. Since the acceleration time is longer than the escape time shown in Figures 14, the expected maximum energy is lower than those assumed in our simulations. For the lower energy particles, $t_{\rm esc,d}$ is shorter than $t_{\rm esc,b}$. Equating the acceleration time and the diffusive escape time, we obtain the maximum energy:

[TABLE]

where $r_{L,i}$ is the Larmor radius for the particles of $\epsilon=\epsilon_{i}$ and we use $\xi=1.3$ and $\zeta=3$ to obtain the value. A higher BH mass makes the system size larger, which helps to accelerate CRs to higher energy. The magnetic field is stronger for a higher $\eta$ or smaller $\chi$, leading to the higher $\epsilon_{\rm max}$. Note that even if we use $t_{\rm esc,b}$ instead of $t_{\rm esc,d}$, the estimate does not drastically change.

For $\chi=10$, which corresponds to $R_{\rm ini}=3R_{s}$, the maximum energy is $\epsilon_{\rm max}\simeq 2$ PeV. With this energy, we may expect production of high-energy neutrinos of $\sim 0.1$ PeV from accretion discs through inelastic hadronic collision and photomeson production (Kimura et al., 2015). Indeed, the astrophysical neutrinos of $0.01-10$ PeV are detected by IceCube (Aartsen et al., 2015a, b), and LLAGNs are a good candidate of the neutrino source (Kimura et al., 2015; Khiali & de Gouveia Dal Pino, 2016). On the other hand, the maximum energy for Galactic X-ray binaries of $10~{}M_{\odot}$ is at most a few TeVs according to Equation (40). This means that hot accretion flows in Galactic X-ray binaries cannot produce neutrinos detected by IceCube through the turbulent acceleration mechanism.

Our results imply that the hot accretion flow cannot accelerate ultrahigh-energy CRs (UHECRs). A higher $\eta$, i.e., higher mass accretion rate, results in forming a standard thin disc (Shakura & Sunyaev, 1973; Ohsuga & Mineshige, 2011), where the Coulomb loss prevents the particles from being accelerated (Kimura et al., 2014, 2015). This condition gives $\eta\lesssim 10$ (see Section 3.1). The initial radius should be larger than a few $R_{s}$ for the CRs to escape from the system. Otherwise, the CRs should fall to the BH because the accretion flow is supersonic (e.g., Abramowicz et al., 1996; Nakamura et al., 1997; Kimura et al., 2014). The maximum mass of a SMBH is expected to be $\sim 10^{10}~{}M_{\odot}$ (Netzer, 2003; Jun et al., 2015; Ichikawa & Inayoshi, 2017). Even with the extreme parameter set ($M=10^{10}~{}M_{\odot}$, $\chi=10$, $\eta=10$), the maximum CR energy is estimated to be $\sim 70$ PeV, which is too low to be the source of UHECRs (see e.g., Kotera & Olinto, 2011, for a review). They should be produced by other sites or sources, such as radio galaxies (e.g., Takahara, 1990; Murase et al., 2012b; Kimura et al., 2018) or gamma-ray bursts (e.g., Waxman, 1995; Murase et al., 2008; Zhang et al., 2018).

4.2 Comparison to the shearing box simulations

Kimura et al. (2016) performed MHD simulations with the shearing box approximation, and investigated the behaviour of CRs using the snapshot data of the shearing box MHD simulations. In the shearing box simulations, the particle acceleration is described by a diffusion phenomenon in energy space. The troidal magnetic field dominates over the poloidal field, so the particles tend to move to the azimuthal direction.

These features are common with the global simulations. However, we find a few discrepancies between the global and the shearing box simulations. The first point is the importance of the shear acceleration (see e.g., Berezhko & Krymskii, 1981; Katz, 1991; Subramanian et al., 1999). In the shearing box simulations, the shear acceleration is effective for higher energy particles. However, the shear acceleration is not observed in the global simulations. This is because the calculation region is finite for the global simulations. The shearing box approximation creates an unrealistic shear velocity due to the periodic boundary condition. With the global simulation data, the high-energy particles escape from the system before the acceleration by the shear. Therefore, the shear acceleration is not effective in the realistic accretion flows.

Another is the bulk outward motion discussed in Section 3.2.1 and 3.2.3. In the global simulations, the magnetic fields have spiral structure owing to the differentially rotating accretion flow. Since the magnetic field is stronger in the inner region, the magnetic mirror force prevents the particles from moving inward. Therefore, the particles tend to move outward when they stream along the spiral magnetic field. By contrast, the shearing boxes do not distinguish radially inward or outward. Although the local simulations can create spiral magnetic fields, there is no magnetic mirror force because the magnetic field strength should be the same in both directions. Thus, they cannot produce the bulk outward motion. Note that the quantitative physical understanding of the energy dependence of $\beta_{\rm CR,bul}$ remains as future work.

The other is the appearance of the superdiffusion. In the shearing box simulations, the spacial diffusion is described by the anisotropic Bohm diffusion, where the diffusion coefficient for the azimuthal direction is higher than those for the other two directions. On the other hand, the global simulations result in a superdiffusive transport. While the superdiffusive transports are observed in some of the previous simulations (Zimbardo et al., 2006; Xu & Yan, 2013; Roh et al., 2016) and have been discussed in the literature of interplanetary turbulence (Perri & Zimbardo, 2007, 2009; Sugiyama & Shiota, 2011), their situations are different from our setup of the turbulent fields. Understanding the cause of the superdiffusion is left as future work.

4.3 Future directions

Our results indicate $D_{\epsilon}\propto\epsilon^{2}$ with $\lambda_{\rm ini}\simeq 0.5-8$. This is different from the index of the turbulent power spectrum in the relevant scales ($m\gtrsim 100$). This means that the resonant scattering is ineffective in our simulation. However, the numerical dissipation suppresses the turbulent power in these small scales. If the turbulence has sufficient power in the smaller scales, the resonant scattering might be important for the lower energy particles. Higher resolution calculations with a higher-order reconstruction scheme is necessary to investigate effects of the waves at the smaller scales.

In the vicinity of the BH, our Newtonian treatment is no longer valid, and effects of general relativity must be important. Many general relativistic MHD (GRMHD) simulations are peformed (e.g., McKinney, 2006; Takahashi et al., 2016). However, GRMHD calculations require more computational time than Newtonian ones, which makes it difficult to perform high resolution simulations. Also, in the GRMHD simulations, the velocity of the background fluid is close to the speed of light. This situation disallows our use of the snapshot data, so we need to perform particle simulations with a dynamically evolving turbulence.

Although we setup the initial condition for MHD simulations with the zero net vertical magnetic flux, there can be a strong large-scale magnetic field (Narayan et al., 2003; Tchekhovskoy et al., 2011; McKinney et al., 2012). This type of magnetic field configuration is expected to be related to the relativistic jet production. The non-thermal particles can play an important role on the mass loading processes to the jet (Toma & Takahara, 2012; Kimura et al., 2014), which would be investigated by the test-particle simulations with GRMHD simulation data (Bacchini et al., 2018).

Hot accretion flows are expected to be composed of collisionless plasma (Takahara & Kusunose, 1985; Quataert & Gruzinov, 1999; Kimura et al., 2014), where the non-ideal MHD processes can be important. Resistivity creates an electric field parallel to a magnetic field in a dissipation scale, which might affect the energy of CRs. The features of the MRI turbulence can be affected by the anisotropic pressure (Quataert et al., 2002; Sharma et al., 2006; Kunz et al., 2014; Hirabayashi & Hoshino, 2017) and/or anisotropic heat transfer (Ressler et al., 2015; Foucart et al., 2017). Magnetic reconnections are expected to have some influence on the MRI turbulence, since they are important dissipation processes (e.g., Sano & Inutsuka, 2001). The reconnections also provide seed CRs (Kowal et al., 2011, 2012; Hoshino, 2012; Ball et al., 2017), which are necessary for the turbulent stochastic acceleration to work as CR sources.

5 Summary

We have investigated behaviours of high-energy relativistic particles in MRI turbulence in accretion flows. We have generated the turbulent fields via MHD simulations to model accretion flows around a black hole, calculated orbits of test-particles by using the snapshot MHD data, and investigated the particle acceleration and escape processes.

We have run four MHD simulations with the same initial condition of an equilibrium torus with different resolutions and grid spacing. For all the simulations, the MRI grows in a few rotation time of the initial torus, and a quasi-steady state is achieved after several rotation time. Due to the shear motion, $B_{\phi}$ is stronger than the other components. The turbulent velocity is sub-sonic and sub-Alfvenic, while it is faster than the radial advection velocity. The mass-density and magnetic-energy fluctuations are weakly anti-correlated, which means that the slow and Alfven modes have important influences on the turbulence. The power spectra of the magnetic fields are anisotropic: $mP_{m}\propto m$ for $B_{r}$ and $B_{\theta}$, while $mP_{m}\propto m^{1/2}$ for $B_{\phi}$, although they have a similar peak around $m\sim 10-20$. These features of the turbulence are common for all the runs, although the lowest-resolution run has weaker turbulent fields because of the lack of resolution.

Using the snapshot data of the MHD simulations, we have calculated a number of orbits of test-particles. We have found that the particles tend to escape through the outer radial boundary, rather than through the vertical boundary or falling to the BH. The particles have anisotropic momentum distribution due to the magnetic field configuration generated by the shearing accretion flow. The evolution of the distribution function in energy space is described by a diffusion phenomenon. The diffusion coefficient is described by Equation (3.2.2), which implies that large-scale waves have a dominant role for changing the particle energy for the energy range we investigated. The evolution in configuration space is not simply written by the diffusion. The particles spread faster than the usual diffusion (superdiffusion). The variance of the radial displacement is approximated by Equation (33). Also, the anisotropic momentum distribution creates the bulk outward motion of the CR particles. This bulk motion is faster for higher energy particles as shown in Figure 13. For higher energy particles, the bulk CR motion is more effective than the diffusive motion, and vice versa for lower energy particles. Physical interpretations of the emergence of the superdiffusion and the energy dependence of the bulk CR motion are left as future work. To obtain a solid conclusion, MHD simulations with a higher-order scheme and an order of magnitude higher resolution is essential.

Acknowledgements

S.S.K. thanks Matthew W. Kunz for useful comments. We are grateful to the anonymous referee for a careful reading of the manuscript and thoughtful comments. This work is supported by JSPS Oversea Research Fellowship, the IGC post-doctoral fellowship program (S.S.K.), Alfred P. Sloan Foundation, NSF Grant No. PHY-1620777 (K.M.), and JSPS KAKENHI Grant Numbers 16H05998 and 16K13786 (K.T.). Numerical computations were carried out on Cray XC30 and Cray XC50 at Center for Computational Astrophysics, National Astronomical Observatory of Japan.

Appendix A Derivation of Relation between $\sigma_{\epsilon}^{2}$ and $D_{\epsilon}$

Under the approximation of $D_{\epsilon}\approx D_{\epsilon_{\rm ini}}$ and $2D_{\epsilon}/\epsilon\approx v_{\epsilon_{\rm ini}}$, Equation (24) is expressed as

[TABLE]

The mean of the momentum is written as $\mu_{\epsilon}=\int N_{\epsilon}\epsilon d\epsilon/N_{\rm actv}$. Its time derivative is

[TABLE]

where we use a partial integration and $N_{\epsilon}\to 0$ for $\epsilon\to\infty$ and $\epsilon\to 0$. Integrating both sides with $t$, we obtain

[TABLE]

A similar calculation gives us the variance $\sigma_{\epsilon}^{2}$. The variance of the momentum is written as $\sigma_{\epsilon}^{2}=\int N_{\epsilon}\epsilon^{2}d\epsilon/N_{\rm actv}-\mu_{\epsilon}^{2}$. Its time derivative is

[TABLE]

Hence, we obtain

[TABLE]

For the special case of $D_{\epsilon}=D_{\epsilon_{\rm ini}}(\epsilon/\epsilon_{\rm ini})^{2}$, we can derive the time evolution of $\mu_{\epsilon}$ and $\sigma_{\epsilon}^{2}$ by similar algebra without assuming that $D_{\epsilon}$ is constant. Using Equation (24), the time derivative of $\mu_{\epsilon}$ is written by

[TABLE]

where $v_{\epsilon_{\rm ini}}=2D_{\epsilon_{\rm ini}}/\epsilon_{\rm ini}$. Then, we obtain

[TABLE]

Hence, the increasing rate of $\mu_{\epsilon}$ is twice higher than that in Equation (43). The time derivative of the mean of $\epsilon^{2}$, $\langle\epsilon^{2}\rangle$, is

[TABLE]

Integrating both sides, we get

[TABLE]

Using Equations (47) and (49), we obtain

[TABLE]

This $\sigma_{\epsilon}^{2}$ is the same as that in Equation (45).

Figure 16 shows the ratio of $D_{\epsilon}$ obtained by $\mu_{\epsilon}$, using either Equation (43) or Equation (47), to those given by $\sigma_{\epsilon}^{2}$. We can see that the values are consistent within a factor of 3 for $\lambda_{\rm ini}\leq 8$ (corresponding to $\epsilon\leq 3.6\times 10^{2}$ PeV for run A) for the cases that $D_{\epsilon}\approx D_{\epsilon_{\rm ini}}$ is constant (blue cross symbols). The agreement is improved for the cases with $D_{\epsilon}\propto\epsilon^{2}$ for $\lambda_{\rm ini}\leq 8$, where the values of $D_{\epsilon}$ agree each other within a factor of 2 (red plus symbols).