GRMHD simulation of black hole accretion disk

Characterization of the production of intense Alfvén pulses : GRMHD simulation of black hole accretion disks


The episodic dynamics of the magnetic eruption of a spinning black hole (BH) accretion disks and its associated intense shapeup of their jets is studied via three-dimensional general-relativistic magnetohydrodynamics (GRMHD). The embedded magnetic fields in the disk get amplified by the magnetorotational instability (MRI) so large as to cause an eruption of magnetic field (recconection) and large chunks of matter episodically accrete toward the roots of the jets upon such an event. We also find that the eruption events produce intensive Alfvén pulses, which propagate through the jets. After the eruption, the disk backs to the weakly magnetic states. Such disk activities cause short time variabilities in mass accretion rate at the event horizon as well as electric-magnetic luminosity inside the jet. Since the dimensionless strength parameter of these Alfvén wave pulses is extremely high for the solar masses central black hole and % Eddington accretion rate accretion flow, the bow wake acceleration model proposed by Ebisuzaki & Tajima (2014) certainly works to accelerate the ultra-high energy cosmic rays and electrons which finally emit gamma-rays. Since our GRMHD model has universality in its spatial and temporary scales, it is applicable to a wide range of astrophysical objects ranging from those of AGN (which is the primary target of this research), to micro-quasars, and down to young stellar objects. Properties such as time variabilities of blazar gamma-ray flares and spectrum observed by Fermi Gamma-ray Observatory are well explained by linear acceleration of electrons by the bow wake.

Subject headings:
magnetohydrodynamics (MHD) — accretion disks — quasars: supermassive black holes — galaxies: jets

1. Introduction

Active galactic Nuclei (AGNs) are high energy astronomical objects, so that they emit non-thermal radiation in any frequency ranges of radiation, in other words, radio, infrared, visible lights, ultraviolet, X-rays, and gamma-rays (Begelman et al., 1984; Hughes, 1991; Burgarella et al., 1993; Tsiganos, 1996; Ferrari, 1998). These central engines are believed to be accreting supermassive black holes, with relativistic jet whose bulk Lorentz factor is (Biretta et al., 1999; Asada et al., 2014; Boccardi et al., 2016). Their jets show strong time variabilities in the timescales from days to years (Fossati et al., 1998; Abdo et al., 2009a, 2010a, 2010b, 2011; Ackermann et al., 2010; Chen et al., 2013). In the extreme cases, blazars show bursts of hours (Ackermann et al., 2016; Britto et al., 2016). These radiations are believed to be emitted by a bunch of electrons with strongly relativistic motions.

The system of AGN jets is also believed as a cosmic ray accelerator (Dermer et al., 2009). Although it has a potential to accelerate highest energy to the energy of eV, it is still not well understood what is physical mechanism for the particle acceleration and where the acceleration site is. Many studies have been done for diffusive shock acceleration model based on conventional Fermi acceleration mechanism (Fermi, 1954). In the Fermi acceleration, charged particles interact with magnetized clouds, and are vented random directions by their magnetic field. Since the head-on collisions, which the particles gain the energy, are more frequent than rear-end collisions, which they loose their energies, the particles statistically gain the energy step by step and eventually obtain very high energy to be cosmic ray particles. However, this Fermi acceleration mechanism is difficult to explain highest energy particles eV, because of 1) the large number of scatterings necessary to reach highest energies, 2) energy losses through the synchrotron emission at the bending associated with scatterings, and 3) difficulty in the escape of particles which are initially magnetically confined in the acceleration domain (e.g., Kotera & Olinto (2011)).

On the other hand, Tajima & Dawson (1979) proposed that particles can be accelerated by the wakefield induced by an intense laser pulse, see review by Tajima et al. (2017). In particular, ponderomotive force which is proportional to gradient of works to accelerate the charged particles, effectively, where is electric field of the electric-magnetic wave. The acceleration towards relativistic regime by the ponderomotive force is confirmed by recent experiments by ultra intense lasers for electrons (Leemans et al., 2006; Nakamura et al., 2007). Positron acceleration up to GeV is also reported by Corde et al. (2015) driven by wakefields by electrons.

Takahashi et al. (2000); Chen et al. (2002); Chang et al. (2009) applied this mechanism to magnetowave induced plasma wakefield acceleration for ultra high energy cosmic rays. Recently, Ebisuzaki & Tajima (2014a, b) applied this wakefield acceleration theory to the relativistic jets launched from an accreting black hole. In such astrophysical context going far beyond the laboratory scales of wakefields (see, for example, Tajima et al. (2017)), the relativistic factors that characterize the dynamics becomes far greater than unity. This has not been achieved in the laboratory yet, though simulations began to peek into it. In this regime (Ebisuzaki & Tajima, 2014a, b) the ponderomotive acceleration has advantages over the Fermi mechanism. In close scrutiny, the wakefields are composed with two parts the frontal bow part and the following stern (wake) part. Here we may call simply the bow wake and stern wake. In the work of high simulation, it was shown in Lau et al. (2015) that the bow wake (it is driven directly by the ponderomotive force) is dominant over the stern wake. The advantages of this bow wake acceleration over the Fermi mechanism are:

  1. The ponderomotive field provides an extremely high accelerating field (including the wakefield).

  2. It does not require particle bending, which would cause strong synchrotron radiation losses in extreme energies.

  3. The accelerating fields and particles move in the collinear direction at the same velocity, the speed of light, so that the acceleration has a built-in coherence called “relativistic coherence” (Tajima, 2010); in contrast, the Fermi acceleration mechanism, based on multiple scatterings, is intrinsically incoherent and stochastic.

  4. No escape problem (Kotera & Olinto, 2011) exists. Particles can escape from the acceleration region since the accelerating fields naturally decay out.

They found that protons can be accelerated even above ZeV eV in the bow wake of a burst of Alfvén waves emitted by an accretion disk around a black hole with the mass of M. Ebisuzaki & Tajima (2014a) used three major assumptions based on the standard -disk model (Shakura & Sunyaev, 1973).

  • Assumption A: the magnetic field energy included in an Alfvén wave burst is assumed as:


    where is the magnetic field stored in the inner most regions of the accretion disk, is the Schwarzschild radius of the black hole, is the thickness of the disk, is gravitational radius, is the accretion rate normalized by the Eddington luminosity, and is the mass of the black hole in the unit of solar mass.

  • Assumption B: they assumed that the angular frequency of the Alfvén wave corresponds to that excited by magnetorotational instability (MRI (Velikhov, 1959; Chandrasekhar, 1960; Balbus & Hawley, 1991; Matsumoto & Tajima, 1995)), which takes place in a magnetized accretion disk, in other words:


    where is the wavelength of the Alfvén wave, and is speed of Alfvén wave.

  • Assumption C: the recurrence rate of the Alfvén burst is evaluated as:


    where is the episode-dependent parameter of the order of unity.

They found that the non-dimensional strength parameter is as high as for the case of and , where is electric charge, is the intensity of the electric field, is mass of electron, and is speed of light. The ponderomotive force of this extremely relativistic waves co-linearly accelerate to the jet particle up to the maximum energy:


where is the charge of the particle and is the bulk Lorentz factor of the jet. Recent one-dimensional particle in cell (PIC) simulation shows maximum energy gain via a ponderomotive force in the bow wake and the maximam energy is almost proportional to (Lau et al., 2015). Based on the above estimation, Ebisuzaki & Tajima (2014a) concluded that the accreting supermassive black hole is the ZeV ( eV) linear accelerator.

Major motivation of the present paper is to verify the assumptions of Ebisuzaki & Tajima (2014a) (Equations (2) and (3)) by the 3D GRMHD simulations of accretion disk around a supermassive black hole. This paper is organized as follows. We describe our physical models and numerical details in § 2. The results are shown in § 3. The application to the ultra high energy cosmic ray acceleration, blazars, and gravitational waves is discussed in § 4 and § 5.

2. General Relativistic Magnetohydrodynamic Simulation Method

2.1. Basic Equations

We numerically solve general relativistic magnetohydrodynamic equations, assuming a fixed metric around a black hole. The unit in which and are unity is adopted, where is gravitational constant, is mass of the central black hole. The scales of length and time are and , respectively. The mass and energy is scale free. The achieved mass accretion rate at the event horizon is used to scale of the mass and energy, for example. The metric around a rotating black hole whose dimensionless spin parameter is can be described by Boyer-Lindquist coordinate or Kerr-Schild coordinate as follows. The line element in BL coordinate is


where , , , , , , , and . We follow standard notation used in Misner et al. (1973), i.e, metric tensor for Minkowski space is diag.

The line element in Kerr-Schild coordinate is


where , , , , , , and .

We also use so-called modified Kerr-Schild coordinate so that the numerical grids are fine near the event horizon. The transformation between Kerr-Schild coordinate and modified Kerr-Schild coordinate is described as , , , and , where is a parameter which controls how the grids are concentrated around the equator but we set in this study, resulting in uniform grids in polar angle. The grid numbers are , , and which are uniformly spaced. The resolution is comparable to those used in recent 3D GRMHD simulations McKinney et al. (2012); Penna et al. (2010). Computational domain covers from inside the event horizon to , in polar angle, and in azimuthal angle.

Contravariant vectors in Boyer-Lindquist coordinate and Kerr-Schild coordinate are related with , , , and . Contravariant vectors in Kerr-Schild coordinate and modified Kerr-Schild coordinate are related with , , , and .

Mass and energy-momentum conservation laws are,


where is energy momentum tensor, is rest mass density, is fluid 4-velocity, is determinant of metric tensor, i.e., , and , and is Christoffel symbol which is defined as . Energy momentum tensor which includes matter and electromagnetic parts is defined as follows


where is specific enthalpy, is thermal energy density, is thermal pressure, and is the Faraday tensor and a factor of is absorbed into the definition of . The dual of the Faraday tensor is


where and is the completely antisymmetric Levi-Civita symbol (). The magnetic field observed by normal observer is


where is the normal observer’s four-velocity and is the lapse. Note the time component of is zero, since . Here we introduce another magnetic field which is used in Gammie et al. (2003); Noble et al. (2006); Nagataki (2009) as


The time component of is also zero. We also introduce four magnetic field which is measured by an observer at rest in the fluid,


and are related with


is satisfied. By using this magnetic four vector electric-magnetic component of energy momentum tensor and the dual of Faraday tensor can be written as


where the magnetic pressure is . The electric-magnetic part of the energy momentum tensor can be written by and as


The Maxwell equations are written as


By using Eqs. (19) these equations give


These are no-monopole constraint equation and time evolution of spacial magnetic field equations, i.e, the induction equations, respectively.

In order to close the equations, ideal gas equation of state is adopted, where is the specific heat ratio which is assumed to be constant (). We ignore self-gravity of the gas around the black hole and any radiative processes, assuming radiatively inefficient accretion flow (RIAF) in the disk (Narayan & Yi, 1994).

We numerically solve these equations by GRMHD code developed by one of authors (Nagataki, 2009, 2011). Magnetohydrodynamic equations are solved by using shock capturing method (HLL method), applying 2nd order interpolation to reconstruct of physical quantities at the cell surfaces and 2nd order time integration by using TVD (total variation diminishing) Runge-Kutta method, see also Gammie et al. (2003); Noble et al. (2006). The boundary conditions are zero gradient for and periodic one for .

2.2. Initial Condition

We adopt the Fishbone-Moncrief solution as an initial condition for hydrodynamic quantities as adopted for recent GRMHD simulations McKinney & Gammie (2004); McKinney (2006); McKinney et al. (2012); Shiokawa et al. (2012). This solution describes the hydro-static torus solution around a rotating black hole. The gravitational force by the central black hole, the centrifugal force and the pressure gradient force balance each other. There are some free parameters to give a solution. We assume the disk inner edge at the equator is at and a constant specific angular momentum (). The 4-velocity is firstly given at Boyer–Lindquist coordinate then transformed to Kerr-Schild one and to modified Kerr-Schild one. The dimensionless spin parameter is assumed to be . The radii of event horizon and the innermost stable circular orbits (ISCO) at the equator are at and , respectively. The initial disk profile is on the plane including polar axis shown in Fig. 1 which shows mass density contour. The disk is geometrically thick and is different from standard accretion disk (Shakura & Sunyaev, 1973) in which the geometrically thin disk is assumed.

Figure 1.— Log scaled rest mass density profile on plane at .

As we have discussed in Sec. 1 magnetic fields play an important role not only for the dynamics of the accretion flows but also for the dynamics of the outflows. We impose initially weak magnetic field inside the disk as a seed. This weak magnetic field violates the initial static situation and is expected to be amplified by winding and MRI.

Here we introduce the four-vector potential of the electromagnetic field. The Faraday tensor is defined by this vector potential as follows.


In this study we assume initially closed poloidal magnetic field, i.e., the toroidal component of the initial vector potential is,


where is maximum mass density in the initial torus. Other spacial components are zero, i.e., . Since the vector potential has only toroidal components, the poloidal magnetic field is imposed. To violate axis-symmetry maximally 5% amplitude random perturbation is imposed in the thermal pressure i.e., thermal pressure is , where is equilibrium thermal pressure derived by Fishbone-Moncrief solution and is random number in the range of . This perturbation violates axis-symmetry of the system and triggers generation of non axisymmetric mode.

3. Episodic Eruption of Disks and Jets

3.1. B-field amplification and mass accretion

The magnetic field in the disk is amplified by winding effect and MRI, as follows. The initially imposed poloidal magnetic field is stretched in the toroidal direction, generating toroidal components, since there is a differential rotation inside the accretion disk. Although initially imposed magnetic field is weak, i.e., the minimum plasma is 100 inside the disk, the strength of the magnetic field quickly increases by the winding effect (Duez et al., 2006). Fig. 2(a) shows the volume averaged strength of the magnetic field at the equator i.e., averaged over , , and from to . Since the magnetic field is stretched by the differential rotation in the disk, toroidal component dominates over the poloidal component, though both components show strong time variability. Fig. 2(b) shows the mass accretion rate at the event horizon , where the mass accretion rate at a radius is defined as


where is volume element, for example, , and the sign is chosen so that the case of mass inflow is positive mass accretion. The mass accretion rate also shows strong time variability, and synchronized with magnetic field strength near the ISCO.

Figure 2.— (a) Time evolution of strength of poloidal and toroidal magnetic fields at the equator averaged at and . (b) Time evolution of mass accretion rate () at the event horizon ().

Bottom two panels in Fig. 3 show the contours of at the equator at and . Between these two figures the state of the disk near the disk inner edge () transitioned from low state () to high state (). Such a transition between low and high states are pointed by Tajima & Gilden (1987); Shibata et al. (1990). Characteristics of these two states are also listed at the top of Fig. 3 (taken from Mineshige et al. (1995)). Middle two panels in Fig. 3 show the toroidal magnetic field lines of two disks from MHD simulations by Tajima & Gilden (1987). At the low state (left panels in Fig. 3) the toroidal magnetic field is stretched to the limit and magnetic field energy is stored. Since some field lines are almost anti-parallel and very close to each other, the reconnection will happen eventually. After magnetic field energy is dissipated via recconection, the system becomes high state (right panels in Fig. 3).

Bar structures near the disk inner edge can be seen (bottom panels in Fig. 3). The non-axis-symmetric mode is excited, as shown in global hydrodynamic and magnetohydrodynamic simulations of accretion disks, for example Tajima & Gilden (1987); Machida & Matsumoto (2003); Kiuchi et al. (2011); McKinney et al. (2012). Figure 4 shows the contours of (at the equator), mass density ( plane), and magnetic pressure ( plane) at two different times as shown in Fig. 3. Along the polar axis low density and highly magnetized region appears. It corresponds to the Poynting flux dominated jet. Thus baryonless and highly magnetized jet is formed along the polar axis. In this region the Alfvén speed is almost speed of light . A disk wind blows between the jet and the accretion disk. Filamentaly structures which are excited by MRI can be seen in the magnetic pressure contours. The thickness of filaments are , as shown in Fig. 4.

Parameter Low disk High disk
( supported) ( supported)
Configuration Optically thin disk Optically thick disk
(advection-dominated) (cooling-dominated)
consisting of blobs + corona
Dissipation of magnetic fields Reconnection Escape via buoyancy
and reconnection
Dissipation of energy Sporadic Continuous
Spectrum Hard power law Soft + hard tail
Fluctuations Large Small
Figure 3.— Table: Properties of low and high states of the disks taken from Mineshige et al. (1995). Middle panels : Toroidal magnetic field lines at low state (left) and high state taken from Tajima & Gilden (1987). Bottom panels : Inverse of plasma beta () contours at the equator shown by logarithmic scales at (low state, left) and at (high state, right.) The shadowed regions in the circle indicate inside the event horizon.
Figure 4.— Contours of inverse of the plasma (at the equator), mass density ( plane), and magnetic pressure ( plane) at two different times as shown in Fig. 3. The domain shows , . .

Both the averaged poloidal component of magnetic field near the ISCO at the equator and the mass accretion rate at the event horizon show synchronous time variability. This is because that the mass accretion rate at the event horizon is strongly affected by the activities of magnetic field amplification near the disk inner edge. The magnetic field amplification via MRI enhances the specific angular momentum transfer inside the accretion disk, resulting in the increase of the accretion rate. This means that the amplification of magnetic fields acts as a viscosity which is introduced as -viscosity in Shakura & Sunyaev (1973). Typical increase timescales are . While the magnetic field is amplified, the mass accretion rate at the event horizon rapidly increases. As we will show later, the outflow properties are also show intense time variability which is strongly related with particle acceleration via wakefield acceleration.

The mass accretion rate repeatedly shows sharp rises followed by gradual falling down. The rising timescale for the quick increase of mass accretion rate corresponds to the value of the growth timescale of MRI


where is order of unity, is growth rate of MRI, and for the fastest growing mode and is Newtonian Keplerian angular velocity. The timescale for the fastest growing MRI mode is


This timescale at around is slightly shorter than the timescale of increase of the strength of magnetic field in the disk and the mass accretion rate at the event horizon. By the analysis of MRI for Newtonian MHD, the angular frequency of the mode is for the fastest growing mode for (parallel to polar axis) direction in Keprelian accretion disk, where is wave number and is Alfvén speed of component. The volume averaged Alfvén speed near the ISCO () at the equator is the order of . Thus the wavelength of this mode is with a grid size near the ISCO and at the equator . Our simulation shows that the rising timescale in mass accretion rate . Corresponding wavelength of MRI is estimated . This is consistent with the filamentary structure shown in magnetic pressure panel in Fig. 4. The wavelength of the fastest growing mode is not completely resolved in our simulation, though physical mechanism of MRI is already shown well. If we resolve it completely with five times or more higher resolution, the rise time would reduce to about .

Since the episodic period of this quick increase is which is about 10 times longer than the Keprelian orbital period near the ISCO (). This timescale for the recurrence is consistent with the analysis by local shearing box (Stone et al., 1996; Suzuki & Inutsuka, 2009; O’Neill et al., 2011).

Along the polar axis funnel nozzles appear (Fig. 4). The outward going electric-magnetic luminosity, opening angle of this jet are also time variable like as the mass accretion rate. The radial velocity just above the black hole and becomes positive at typically , but sometimes .

Figure 5(b) shows radial electric-magnetic luminosity calculated by the area integration only around the polar region () at the radius . Electric-magnetic luminosity shows similar short time variability with the magnetic field amplification near the disk inner edge and the mass accretion rate at the event horizon. Figure. 5(c) shows enlarged properties of time evolution of the Fig. 5(b). Some timescales of rising flares and repeat cycle of flares are written. We can see typical rising timescale of the flares is as same as that for rising timescale of magnetic field in the disk, i.e., and the typical cycle of flares is also as same as repeat cycle of magnetic field amplification . We have observed some active phases in electric-magnetic luminosity in the jet, including the epoch shown in Fig. 5(b). The time evolution of electric-magnetic luminosity shows the most active phase from to . In this phase, the averaged radial electric-magnetic flux reaches the averaged disk Alfvén flux at the equator as shown in the Fig. 5 (a). The disk Alfvén flux at the equator is evaluated as an average of component of Alfvén energy flux at the equator times half of Alfvén speed inside the disk (). This is consistent with the assumption A, as we discuss in the next section.

Figure 5.— (a): time evolution of area averaged radial electric-magnetic flux at and (red). Time evolution of averaged electric-magnetic flux inside the disk calculated by at the equator and at . (b) Time evolution of electric-magnetic (Poynting) power in the jet at calculated by the area integration of electric-magnetic flux at . (c) same as (b) but only very active phase is shown. Typical rising and repeat timescales of flares, and , respectively. are presented in the figure.

4. Particle Acceleration

As shown in Fig. 5 (b) flares in electric-magnetic power in the jet are observed, where the Alfvén speed is almost speed of light because of the low mass density. Large amplitude Alfvén waves become electric-magnetic waves by mode conversion of strongly relativistic waves (Daniel & Tajima, 1997, 1998; Ebisuzaki & Tajima, 2014a). The interaction of the electric-magnetic waves and the plasma can result in the acceleration of the charged particles by the ponderomotive force i.e., wakefield acceleration (Tajima & Dawson, 1979). The key for the efficient wakefield acceleration is the Lorentz invariant dimensionless strength parameter of the wave (Esarey et al., 2009),


The velocity of the oscillation motion of the charged particles via electric field becomes speed of light, when . If the strength parameter highly exceeds unity, the ponderomotive force works to accelerate the charged particles to relativistic regime to wave propagating direction.

4.1. Comparison with Ebisuzaki & Tajima (2014a)

In order to evaluate , Ebisuzaki & Tajima (2014a) used three assumptions A, B and C. Based on the results of numerical simulation, we intend to confirm three assumptions. First, assumption A tells us that the Alfvén flux in the jet is equal to that in the disk. As shown in Fig. 5(a), these fluxes are almost at same level at the active phase of the electric magnetic luminosity in the jet. Thus Alfvén waves emitted from the disk via Alfvén burst goes to the jet as assumed in Ebisuzaki & Tajima (2014a), in other words, assumption A.

Second, Ebisuzaki & Tajima (2014a) assumed that magnetic field amplification occurs at for the standard disk model (Shakura & Sunyaev, 1973). In Ebisuzaki & Tajima (2014a) the timescale of the magnetic field amplification () and frequency of the Alfvén wave are determined by the MRI growth rate (Eq. (2)). Although they evaluated it around , as shown in Fig. 3, the transition between low and high plasma disk states occurs not at , but around in our numerical calculations. Since magnetic field amplification mainly occurs inside compared with the assumption by Ebisuzaki & Tajima (2014a), the timescales are shorter than those of them due to faster rotation period. If we apply instead of for Ebisuzaki & Tajima (2014a) model, the timescales are in good agreement with both our numerical results as shown in table 1. The reason why the magnetic field amplification at smaller radius may be due to the high spinning of the black hole, i.e., for which both the event horizon and the radius of ISCO is smaller than those for non-spinning black hole case (). This timescale is well consistent to the rising timescales of blazar flares observed for 3C454.3 which will be discussed in next subsection. In other words, assumption B is OK in qualitatively, but Eq. (2) must to be   (see also table  2).

Finally, Ebisuzaki & Tajima (2014a) estimated the repeat timescale as the crossing time of the Alfvén wave in the disk, i.e, (assumption C) for the standard disk model (Shakura & Sunyaev, 1973). When we apply the radius instead of in Eq. (3), we obtain , ignoring factor which is order of unity. This value is good agreement with our typical . The case of is also listed in the table 1.

We can reevaluate the strength parameter as


Here electric field is estimated as . The angular frequency of the Alfvén wave is , where is assumed to be 0.1 . We use the values at and the time averaged mass accretion rate ) at the event horizon is used as a normalization to 10% Eddington luminosity for solar masses. The estimated strength parameter highly exceeds unity as discussed in Ebisuzaki & Tajima (2014a). This suggests efficient particle acceleration via wakefield acceleration can occur in the jet. Since both electrons and protons are accelerated at the large amplitude electric-magnetic flares via pomderomotive forces and these particle move with the waves, high energy non-thermal electrons are concentrated at these waves.

If we apply the radius for the magnetic field amplification at instead of for Eq. (3), the estimated values such as the angular frequency of Alfvén wave in the jet (), recurrence rate of the burst (), acceleration time , maximum energy of accelerated particle , total accretion power , Alfvén luminosity in the jet, gamma-ray luminosity , and UHECR luminosity in Ebisuzaki & Tajima (2014a) are revised. Table 2 is revised version of Table 1 in Ebisuzaki & Tajima (2014a). Figure 6 is the revised version of Fig. 4 in Ebisuzaki & Tajima (2014a) which shows the relation of the maximum energy of accelerated particle as a function of mass of central black hole and accretion power . Since we do not consider any radiative processes, i.e., RIAF ( Eddington luminosity), the gray shadowed region shows the objects for the UHECR ( eV) accelerators.

our results 3C454.3 AO0235+164 ET14 ET14
rising timescale of flares () 50 92
repeat cycle of flares () 354

a. Abdo et al. (2011)
b. Abdo et al. (2010a)

Table 1Comparison the timescales unit in of rising flares and repeat cycle of flares between our numerical results, blazar observations, and Ebisuzaki & Tajima (2014a) (ET14). Black hole masses (Bonnoli et al., 2011), and (Liu et al., 2006), are used. For Ebisuzaki & Tajima (2014a) model, two different radii ( and (their original assumption)) at which magnetic field amplification occurs are considered.
Values Units
erg s
erg s
erg s
erg s

, , , and .

Table 2Revised version of the Table 1 in Ebisuzaki & Tajima (2014a), where angular frequency in the jet is , recurrence rate of the burst is , acceleration length is , maximum energy of accelerated particle is , is the charge of the particle, is bulk Lorentz factor of the jet, total accretion power is , Alfvén luminosity in the jet is , gamma-ray luminosity is , and UHECR luminosity is . instead of is assumed as a radius where the magnetic field amplification occurs.
Figure 6.— Revised version of Fig. 4 in Ebisuzaki & Tajima (2014a) applying the radius instead of as the magnetic field amplification site for the standard disk model (Shakura & Sunyaev, 1973). Solid lines represent maximum energy of accelerated particle for the energies , and eV via pondermotove acceleration on the plane of the central black hole mass and accretion power , assuming charge of the accelerated particle , and bulk Lorentz factor of the jet (). Dashed lines show 10%, 0.1% and 0.001% accretion rate to the Eddington luminosity. Since we do not consider any radiative processes, i.e., RIAF ( Eddington luminosity), the gray shadowed region shows the objects for the UHECR ( eV) accelerators by our models.

4.2. gamma-rays

Blazar is a subclass of AGNs for which the jets are close to our line of sight. A blazar jet is very bright due to relativistic beaming effect, including high energy gamma-ray bands. They are observed in multi-wavelength from radio to TeV gamma rays. Short time variabilities and polarization are observed, including recent AGILE and Fermi observations for high energy gamma ray bands (Abdo et al., 2009a, 2010a, 2010b; Ackermann et al., 2010; Striani et al., 2010; Abdo et al., 2011; Bonnoli et al., 2011; Ackermann et al., 2012; Chen et al., 2013; Ackermann et al., 2016; Britto et al., 2016).

These non-thermal emissions are usually explained by the internal shock model (Rees, 1978), i.e., two shell collisions (Spada et al., 2001; Kino et al., 2004; Mimica & Aloy, 2010; Pe’er et al., 2016) for which a rapid shell catches up a slow shell, forming relativistic shocks. At the shocks particle accelerations, generating non-thermal particles, are expected by Fermi acceleration mechanism. Finally non-thermal emission is produced via synchrotron emission and inverse Compton emission.

Our model can naturally explain properties of observed active gamma-ray flares, i.e., spectrum and timescales of flares. For electrons energy loss via synchrotron radiation causes a cutoff around PeV regime (Ebisuzaki & Tajima, 2014a), although heavier particles, such as protons and heavier nuclei are accelerated up to ultra high energy cosmic ray regime ( eV) and beyond. The accelerated electrons emit radiation from radio to high energy gamma-rays via synchrotron radiation and inverse Compton emission mechanism. The distribution of accelerated non-thermal particles becomes power law with a power law index (Mima et al., 1991), which is consistent with the observed blazar spectrum with the power law index close to . The photon index becomes close to , when the light curve of gamma-rays becomes active phase (Abdo et al., 2010a, 2011; Britto et al., 2016). This anti-correlation between the gamma-ray light curves and the photon index also supports our results.

From our numerical simulations the rising timescale of electric-magnetic flares in the jet is as same as rising timescale of magnetic field amplification in the disk, i.e., typically and timescales of peak to peak in the flares are as same as timescale of repeat cycle of magnetic field amplification in the disk, i.e., typically .

As for comparison with observed gamma-ray flares of blazars the rising timescale of flares and timescales of the cycle of the flares are normalized by the , where is cosmological redshift of the object. The redshifts for two objects are (Lynds, 1967) and . From the observations of line widths of in the broad line region (BLR), the mass of central black hole in 3C454.3 is estimated from (Bonnoli et al., 2011) to (Gu et al., 2001). In this paper we adopt as a mass of central black hole in 3C454.3, since the estimation of the BLR was done by using C line with less contamination by the non-thermal continuum. The mass of central black hole in AO0235+164 is also derived to (Liu et al., 2006) from the line width of in the BLR.

We adopted 7 days as repeat time and 3 days as rise time for 3C454.3, though various timescales with different times are observed from 3C454.3. Among them, the seven days flare observed by Abdo et al. (2011) is the most energetic. The estimated apparent isotropic gamma-ray energy in the seven days flare is 4 times or more higher than those of flares reported in Abdo et al. (2010a); Britto et al. (2016). Sub-energetic and shorter timescale flares reported in Striani et al. (2010); Ackermann et al. (2012); Britto et al. (2016) can be explained the result of the magnetic eruption via reconnection at the smaller region in the accretion flows.

For AO0235+164 we compare the flare reported in Abdo et al. (2010a), since the flare is the most energetic one in the apparent isotropic energy compared with other flares of AO0235+164, for example Ackermann et al. (2012). The rising timescale is three weeks and repeat timescale is four weeks. Table 1 summarize of the comparison between our results, theoretical model by Ebisuzaki & Tajima (2014a) and observations. Both timescales for 3C454.3 are good agreement with our results. For AO0235+164 the timescales are longer than those for our results which suggests the magnetic field amplification may occur outward where the timescale of MRI growth becomes longer.


We have performed 3D GRMHD simulation of accretion flows around a spinning black hole () in order to study the AGN jets from the system of supermassive black hole and surrounding accretion disk as an ultra high energy cosmic ray accelerator via wake field acceleration mechanism.

We start our simulation from a hydrostatic disk, i.e., Fishbone-Moncrief solution with a weak magnetic field and random perturbation in thermal pressure which violate the hydrostatic state and axis-symmetry. We follow the time evolution of the system until . Initially imposed magnetic field is well amplified, due to differential rotation of the disk. Non axis-symmetric mode, i.e., bar mode near the disk edge grows up. The typical timescale of the magnetic field growth near the disk edge is which corresponds to inverse of the growth rate of the MRI. Amplified magnetic field once drops then grows up again. The typical repeat timescale is which corresponds to the analysis by high resolution local shearing box simulations. The transition between low state and high state repeats. This short time variability for the growth of poloidal magnetic field neat the disk edge also can be seen in the mass accretion rate at the event horizon which means the mass accretion seen in our numerical simulation triggered by angular momentum transfer by the growth of the magnetic fields.

We have two types outflows as shown in Fig. 3. One is the low density, magnetized, and collimated outflow, i.e, jets. The other is disk winds which are dense gas flow and between the disk surface and jets. In the jet short time variabilities of the electric-magnetic luminosity are observed. The timescales are similar with those seen in the mass accretion rate at the event horizon and and poloidal Alfvén energy flux in the disk near the ISCO, i.e., typical rising timescales of flares and typical repeat cycle are as same as the rising timescales of the magnetic field amplification and repeat cycle of the magnetic field amplification, and , respectively. Thus short pulsed relativistic Alfvén waves are emitted from the accretion disk, when a part of stored magnetic field energy is released. Since the strength parameter of these waves extremely high as for the solar masses central black hole and 10% Eddington accretion rate accretion flows, the wakefield acceleration proposed by Tajima & Dawson (1979) can be applied in the jet. There are some advantages against Fermi acceleration mechanism (Fermi, 1954). When we apply this mechanism to the cosmic ray acceleration, the highest energy of cosmic ray reaches eV which is enough high to explain the ultra high energy cosmic rays. We observe magnetic field amplification occurs around . If we apply the model by Ebisuzaki & Tajima (2014a) assuming that magnetic field amplification occurs much inside the disk. The two timescales are consistent with our numeral results.

Since both protons and electrons are accelerated via ponderomotive force in the jet. High energy gamma-ray emission are observed if we see the jet almost on-axis, i.e., blazars. The observed gamma-ray flare timescales such as rising timescales of flares and repeat cycle of flares for 3C454.3 by Fermi Gamma-ray Observatory (Abdo et al., 2010a) are well explained by our bow wake acceleration model.

Our numerical resolution used in this paper is still not large enough to resolve the fastest growing mode of MRI in the disk. According to our estimation, at least 5-10 time more better resolution near the equator is necessary to capture the fastest growing mode of MRI Such high resolution calculations by increasing the grid number of polar coordinate and also by increasing the grid number of both radial and azimutial coordinate are possible by using K-computer, PEZY computer 8 which is next generation supercomputer, and post-K computer9.

Lastly, the consequences from our present work include the following implication on the gravitational observation. Since non-axis-symmetric mode grows in the disk, mass accretion onto the black hole causes the emission of gravitational waves (Kiuchi et al., 2011). We estimate the levels of the signal of these gravitational waves, assuming the black hole and accretion disk system. The dimensionless amplitude of the gravitational wave at the coalescing phase can be estimated as Matsubayashi et al. (2004):


where is efficiency and we here assume as , is reduced mass of the unit of solar mass for the black hole and the blob mass of the blob. The mass of the blob is estimated by . Figure 7 (a) shows time evolution of estimated amplitude of the gravitational waves from our mass accretion rate for the blazar 3C454.3.

Figure 7.— Top: Time evolution of amplitude of gravitational waves for 3C454.3 () derived by the mass accretion rate shown in Fig. 2. Bottom : Estimated gravitational wave signal when gas blobs accrete onto the black holes. Sensitivity curves of grand based gravitational detector (KAGRA) and proposed space gravitational detectors (eLISA, preDECIGO, DECIGO, and BBO) are also presented.

Figure 7 (b) shows the estimated amplitude of the gravitational wave for some objects, such as gamma-ray active blazars AO0235+164 () and 3C454.3 (), nearby AGN jet M87, and famous stellar black hole Cygnus X-1. We assume typical frequency of the gravitational wave signal is , where is redshift of the object. Approximated sensitivity curves of the KAGRA (Nakano et al., 2015) proposed space gravitational wave detectors such as eLISA (Klein et al., 2016), preDECIGO Nakamura et al. (2016), DECIGO (Kawamura et al., 2006) and BBO (Yagi & Seto, 2011) are also presented. The signal level is by far small compared to the limit of the presently operating or proposed gravitational antennas.

We thank H. Shinkai for introducing some references for the fitting formula on the sensitivities of gravitational wave detectors. This work was carried out on Hokusai-greatwave system at RIKEN and XC30 at CFCA at NAOJ. This work was supported in part by the Grants-in-Aid of the Ministry of Education, Science, Culture, and Sport 15K17670 (MA), 26287056 (AM & SN), Grant-in-Aid for Scientific Research on Innovative Areas Grant Number 26106006 (TE), the Mitsubishi Foundation (SN), Associate Chief Scientist Program of RIKEN SN, a RIKEN pioneering project ‘Interdisciplinary Theoretical Science (iTHES)’ (SN), and ’Interdisciplinary Theoretical & Mathematical Science Program (iTHEMS)’ of RIKEN (SN).


  1. affiliation: Computational Astrophysics Laboratory, RIKEN, 2-1 Hirosawa, Wako, Saitama, 351-0198, Japan
  2. affiliation: Computational Astrophysics Laboratory, RIKEN, 2-1 Hirosawa, Wako, Saitama, 351-0198, Japan
  3. affiliation: Department of Physics and Astronomy, University of California at Irvine, Irvine, CA 92697, United States
  4. affiliation: Astrophysical BigBang Laboratory, RIKEN, 2-1 Hirosawa, Wako, Saitama, 351-0198, Japan
  5. affiliation: Interdisciplinary Theoretical & Mathematical Science Program (iTHEMS), RIKEN
  6. affiliation: Interdisciplinary Theoretical Science Research Group (iTHES), RIKEN
  7. slugcomment: submitted to ApJ


  1. Abdo, A. A., Ackermann, M., Ajello, M., et al. 2009, ApJ, 700, 597
  2. Abdo, A. A., Ackermann, M., Ajello, M., et al. 2010a, ApJ, 710, 1271
  3. Abdo, A. A., Ackermann, M., Agudo, I., et al. 2010b , ApJ, 716, 30
  4. Abdo, A. A., Ackermann, M., Ajello, M., et al. 2011, ApJ, 733, L26
  5. Ackermann, M., Ajello, M., Baldini, L., et al. 2010, ApJ, 721, 1383
  6. Ackermann, M., Ajello, M., Ballet, J., et al. 2012, ApJ, 751, 159
  7. Ackermann, M., Anantua, R., Asano, K., et al. 2016, ApJ, 824, L20
  8. Asada, K., Nakamura, M., Doi, A., Nagai, H., & Inoue, M. 2014, ApJ, 781, L2
  9. Balbus, S. A., & Hawley, J. F. 1991, ApJ, 376, 214
  10. Begelman, M. C., Blandford, R. D., & Rees, M. J. 1984, Reviews of Modern Physics, 56, 255
  11. Biretta, J. A., Sparks, W. B., & Macchetto, F. 1999, ApJ, 520, 621
  12. Boccardi, B., Krichbaum, T. P., Bach, U., et al. 2016, A&A, 585, A33
  13. Bonnoli, G., Ghisellini, G., Foschini, L., Tavecchio, F., & Ghirlanda, G. 2011, MNRAS, 410, 368
  14. Britto, R. J., Bottacini, E., Lott, B., Razzaque, S., & Buson, S. 2016, ApJ, 830, 162
  15. Burgarella, D., Livio, M., & O’Dea, C. P. 1993, Astrophysical jets : proceedings of the Astrophysical Jets Meeting, Baltimore, 1992 May 12-14 , by Burgarella, Denis.;Livio, Mario;O’Dea, Christopher P.  Cambridge ; New York : Cambridge University Press, 1993. Space Telescope Science Institute symposium series ; 6. Space Telescope Science Institute (U.S.),
  16. Chandrasekhar, S. 1960, Proceedings of the National Academy of Science, 46, 253
  17. Chang, F.-Y., Chen, P., Lin, G.-L., Noble, R., & Sydora, R. 2009, Physical Review Letters, 102, 111101
  18. Chen, P., Tajima, T., & Takahashi, Y. 2002, Physical Review Letters, 89, 161101
  19. Chen, L.-E., Li, H.-Z., Yi, T.-F., Zhou, S.-B., & Li, K.-Y. 2013, Research in Astronomy and Astrophysics, 13, 5
  20. Corde, S., Adli, E., Allen, J. M., et al. 2015, Nature, 524, 442
  21. Daniel, J., & Tajima, T. 1997, Phys. Rev. D, 55, 5193
  22. Daniel, J., & Tajima, T. 1998, ApJ, 498, 296
  23. Dermer, C. D., Razzaque, S., Finke, J. D., & Atoyan, A. 2009, New Journal of Physics, 11, 065016
  24. Duez, M. D., Liu, Y. T., Shapiro, S. L., Shibata, M., & Stephens, B. C. 2006, Phys. Rev. D, 73, 104015
  25. Ebisuzaki, T., & Tajima, T. 2014a, Astroparticle Physics, 56, 9
  26. Ebisuzaki, T., & Tajima, T. 2014b, European Physical Journal Special Topics, 223, 1113
  27. Esarey, E., Schroeder, C. B., & Leemans, W. P. 2009, Reviews of Modern Physics, 81,
  28. Fermi, E. 1954, ApJ, 119, 1
  29. Ferrari, A. 1998, ARA&A, 36, 539
  30. Fossati, G., Maraschi, L., Celotti, A., Comastri, A., & Ghisellini, G. 1998, MNRAS, 299, 433
  31. Gammie, C. F., McKinney, J. C., & Tóth, G. 2003, ApJ, 589, 444
  32. Gu, M., Cao, X., & Jiang, D. R. 2001, MNRAS, 327, 1111
  33. Kawamura, S., Nakamura, T., Ando, M., et al. 2006, Classical and Quantum Gravity, 23, S125
  34. Kino, M., Mizuta, A., & Yamada, S. 2004, ApJ, 611, 1021
  35. Kiuchi, K., Shibata, M., Montero, P. J., & Font, J. A. 2011, Physical Review Letters, 106, 251102
  36. Klein, A., Barausse, E., Sesana, A., et al. 2016, Phys. Rev. D, 93, 024003
  37. Kotera, K., & Olinto, A. V. 2011, ARA&A, 49, 119
  38. Hughes, P. A. (ed.) 1991, Beams and jets in astrophysics, Cambridge Univ. Press
  39. Lada, C. J. 1985, ARA&A, 23, 267
  40. Lau, C. K., Yeh, P. C., Luk, O., et al. 2015, Physical Review Special Topics Accelerators and Beams, 18, 024401
  41. Leemans, W. P., Nagler, B., Gonsalves, A. J., et al. 2006, Nature Physics, 2, 696
  42. Liu, F. K., Zhao, G., & Wu, X.-B. 2006, ApJ, 650, 749
  43. Lynds, C. R. 1967, ApJ, 147, 837
  44. Machida, M., & Matsumoto, R. 2003, ApJ, 585, 429
  45. Matsubayashi, T., Shinkai, H.-a., & Ebisuzaki, T. 2004, ApJ, 614, 864
  46. Matsumoto, R., & Tajima, T. 1995, ApJ, 445, 767
  47. Matsumoto, R., Uchida, Y., Hirose, S., et al. 1996, ApJ, 461, 115
  48. McKinney, J. C. 2006, MNRAS, 368, 1561
  49. McKinney, J. C., & Gammie, C. F. 2004, ApJ, 611, 977
  50. McKinney, J. C., Tchekhovskoy, A., & Blandford, R. D. 2012, MNRAS, 423, 3083
  51. Mima, K., Horton, W., Tajima, T., & Hasegawa, A. 1991, American Institute of Physics Conference Series, 230, 27
  52. Mimica, P., & Aloy, M. A. 2010, MNRAS, 401, 525
  53. Mineshige, S., Kusnose, M., & Matsumoto, R. 1995, ApJ, 445, L43
  54. Mirabel, I. F., & Rodríguez, L. F. 1999, ARA&A, 37, 409
  55. Misner, C. W., Thorne, K. S., & Wheeler, J. A. 1973, San Francisco: W.H. Freeman and Co., 1973,
  56. Nagataki, S. 2009, ApJ, 704, 937
  57. Nagataki, S. 2011, PASJ, 63, 1243
  58. Nakamura, K., Nagler, B., Tóth, C., et al. 2007, Physics of Plasmas, 14, 056708
  59. Nakamura, T., Ando, M., Kinugawa, T., et al. 2016, Progress of Theoretical and Experimental Physics, 2016, 093E01 p
  60. Nakano, H., Tanaka, T., & Nakamura, T. 2015, Phys. Rev. D, 92, 064003
  61. Narayan, R., & Yi, I. 1994, ApJ, 428, L13
  62. Noble, S. C., Gammie, C. F., McKinney, J. C., & Del Zanna, L. 2006, ApJ, 641, 626
  63. O’Neill, S. M., Reynolds, C. S., Miller, M. C., & Sorathia, K. A. 2011, ApJ, 736, 107
  64. Pe’er, A., Long, K., & Casella, P. 2016, arXiv:1610.08712
  65. Penna, R. F., McKinney, J. C., Narayan, R., et al. 2010, MNRAS, 408, 752
  66. Rees, M. J. 1978, MNRAS, 184, 61P
  67. Shakura, N. I., & Sunyaev, R. A. 1973, A&A, 24, 337
  68. Shibata, K., Tajima, T., & Matsumoto, R. 1990a, ApJ, 350, 295
  69. Shiokawa, H., Dolence, J. C., Gammie, C. F., & Noble, S. C. 2012, ApJ, 744, 187
  70. Spada, M., Ghisellini, G., Lazzati, D., & Celotti, A. 2001, MNRAS, 325, 1559
  71. Stone, J. M., Hawley, J. F., Gammie, C. F., & Balbus, S. A. 1996, ApJ, 463, 656
  72. Striani, E., Vercellone, S., Tavani, M., et al. 2010, ApJ, 718, 455
  73. Suzuki, T. K., & Inutsuka, S.-i. 2009, ApJ, 691, L49
  74. Tajima, T. 2010, Proceeding of the Japan Academy, Series B, 86, 147
  75. Tajima, T., Nakajima, K., & Mourou, G. 2017, Nuovo Cimento Rivista Serie, 40, 33
  76. Tajima, T., & Dawson, J. M. 1979, Physical Review Letters, 43, 267
  77. Tajima, T., & Gilden, D. 1987, ApJ, 320, 741
  78. Takahashi Y., Hillman L. W, & Tajima T., 2000, in High Field Science, ed. Tajima T., Mima K., Baldis H. (Klewer), 171
  79. Tsiganos K. (ed.), 1996, Solar and Astrophysical Magnetohydrodynamic Flows, Dordrecht: Kluwer
  80. Yagi, K., & Seto, N. 2011, Phys. Rev. D, 83, 044011
  81. Velikhov E. P. 1959, ZhETF, 36, 1398
Comments 0
Request Comment
You are adding the first comment!
How to quickly get a good reply:
  • Give credit where it’s due by listing out the positive aspects of a paper before getting into which changes should be made.
  • Be specific in your critique, and provide supporting evidence with appropriate references to substantiate general statements.
  • Your comment should inspire ideas to flow and help the author improves the paper.

The better we are at sharing our knowledge with each other, the faster we move forward.
The feedback must be of minimum 40 characters and the title a minimum of 5 characters
Add comment
Loading ...
This is a comment super asjknd jkasnjk adsnkj
The feedback must be of minumum 40 characters
The feedback must be of minumum 40 characters

You are asking your first question!
How to quickly get a good answer:
  • Keep your question short and to the point
  • Check for grammar or spelling errors.
  • Phrase it like a question
Test description