Characterization of the production of intense Alfvén pulses : GRMHD simulation of black hole accretion disks
Abstract
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 threedimensional generalrelativistic 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 electricmagnetic 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 ultrahigh energy cosmic rays and electrons which finally emit gammarays. 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 microquasars, and down to young stellar objects. Properties such as time variabilities of blazar gammaray flares and spectrum observed by Fermi Gammaray Observatory are well explained by linear acceleration of electrons by the bow wake.
Subject headings:
magnetohydrodynamics (MHD) — accretion disks — quasars: supermassive black holes — galaxies: jets1. Introduction
Active galactic Nuclei (AGNs) are high energy astronomical objects, so that they emit nonthermal radiation in any frequency ranges of radiation, in other words, radio, infrared, visible lights, ultraviolet, Xrays, and gammarays (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 headon collisions, which the particles gain the energy, are more frequent than rearend 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 electricmagnetic 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:

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

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

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

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:
(1) 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:
(2) 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:
(3) where is the episodedependent parameter of the order of unity.
They found that the nondimensional 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 colinearly accelerate to the jet particle up to the maximum energy:
(4) 
where is the charge of the particle and is the bulk Lorentz factor of the jet. Recent onedimensional 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 BoyerLindquist coordinate or KerrSchild coordinate as follows. The line element in BL coordinate is
(5) 
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 KerrSchild coordinate is
(6)  
where , , , , , , and .
We also use socalled modified KerrSchild coordinate so that the numerical grids are fine near the event horizon. The transformation between KerrSchild coordinate and modified KerrSchild 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 BoyerLindquist coordinate and KerrSchild coordinate are related with , , , and . Contravariant vectors in KerrSchild coordinate and modified KerrSchild coordinate are related with , , , and .
Mass and energymomentum conservation laws are,
(7)  
(8) 
where is energy momentum tensor, is rest mass density, is fluid 4velocity, 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
(9)  
(10)  
(11) 
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
(12) 
where and is the completely antisymmetric LeviCivita symbol (). The magnetic field observed by normal observer is
(13) 
where is the normal observer’s fourvelocity 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
(14) 
The time component of is also zero. We also introduce four magnetic field which is measured by an observer at rest in the fluid,
(15) 
and are related with
(16)  
(17) 
is satisfied. By using this magnetic four vector electricmagnetic component of energy momentum tensor and the dual of Faraday tensor can be written as
(18)  
(19) 
where the magnetic pressure is . The electricmagnetic part of the energy momentum tensor can be written by and as
(20) 
The Maxwell equations are written as
(21) 
By using Eqs. (19) these equations give
(22)  
(23) 
These are nomonopole 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 selfgravity 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) RungeKutta 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 FishboneMoncrief 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 hydrostatic 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 4velocity is firstly given at BoyerâLindquist coordinate then transformed to KerrSchild one and to modified KerrSchild 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.
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 fourvector potential of the electromagnetic field. The Faraday tensor is defined by this vector potential as follows.
(24) 
In this study we assume initially closed poloidal magnetic field, i.e., the toroidal component of the initial vector potential is,
(25) 
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 axissymmetry maximally 5% amplitude random perturbation is imposed in the thermal pressure i.e., thermal pressure is , where is equilibrium thermal pressure derived by FishboneMoncrief solution and is random number in the range of . This perturbation violates axissymmetry of the system and triggers generation of non axisymmetric mode.
3. Episodic Eruption of Disks and Jets
3.1. Bfield 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
(26) 
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.
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 antiparallel 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 nonaxissymmetric 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 
(advectiondominated)  (coolingdominated)  
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 
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
(27) 
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
(28) 
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 electricmagnetic 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 electricmagnetic luminosity calculated by the area integration only around the polar region () at the radius . Electricmagnetic 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 electricmagnetic luminosity in the jet, including the epoch shown in Fig. 5(b). The time evolution of electricmagnetic luminosity shows the most active phase from to . In this phase, the averaged radial electricmagnetic 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.
4. Particle Acceleration
As shown in Fig. 5 (b) flares in electricmagnetic 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 electricmagnetic waves by mode conversion of strongly relativistic waves (Daniel & Tajima, 1997, 1998; Ebisuzaki & Tajima, 2014a). The interaction of the electricmagnetic 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),
(29) 
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 nonspinning 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
(30) 
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 electricmagnetic flares via pomderomotive forces and these particle move with the waves, high energy nonthermal 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, gammaray 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)
Values  Units  

s  
s  
s  
eV  
erg s  
erg s  
erg s  
erg s  
–  
– 
, , , and .
4.2. gammarays
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 gammaray bands. They are observed in multiwavelength 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 nonthermal 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 nonthermal particles, are expected by Fermi acceleration mechanism. Finally nonthermal emission is produced via synchrotron emission and inverse Compton emission.
Our model can naturally explain properties of observed active gammaray 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 gammarays via synchrotron radiation and inverse Compton emission mechanism. The distribution of accelerated nonthermal 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 gammarays becomes active phase (Abdo et al., 2010a, 2011; Britto et al., 2016). This anticorrelation between the gammaray light curves and the photon index also supports our results.
From our numerical simulations the rising timescale of electricmagnetic 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 gammaray 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 nonthermal 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 gammaray 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). Subenergetic 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.
5. DISCUSSIONS and SUMMARY
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., FishboneMoncrief solution with a weak magnetic field and random perturbation in thermal pressure which violate the hydrostatic state and axissymmetry. We follow the time evolution of the system until . Initially imposed magnetic field is well amplified, due to differential rotation of the disk. Non axissymmetric 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 electricmagnetic 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 gammaray emission are observed if we see the jet almost onaxis, i.e., blazars. The observed gammaray flare timescales such as rising timescales of flares and repeat cycle of flares for 3C454.3 by Fermi Gammaray 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 510 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
Kcomputer, PEZY computer
Lastly, the consequences from our present work include the following implication on the gravitational observation. Since nonaxissymmetric 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):
(31) 
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 (b) shows the estimated amplitude of the gravitational wave for some objects, such as gammaray active blazars AO0235+164 () and 3C454.3 (), nearby AGN jet M87, and famous stellar black hole Cygnus X1. 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 Hokusaigreatwave system at RIKEN and XC30 at CFCA at NAOJ. This work was supported in part by the GrantsinAid of the Ministry of Education, Science, Culture, and Sport 15K17670 (MA), 26287056 (AM & SN), GrantinAid 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).
Footnotes
 affiliation: Computational Astrophysics Laboratory, RIKEN, 21 Hirosawa, Wako, Saitama, 3510198, Japan
 affiliation: Computational Astrophysics Laboratory, RIKEN, 21 Hirosawa, Wako, Saitama, 3510198, Japan
 affiliation: Department of Physics and Astronomy, University of California at Irvine, Irvine, CA 92697, United States
 affiliation: Astrophysical BigBang Laboratory, RIKEN, 21 Hirosawa, Wako, Saitama, 3510198, Japan
 affiliation: Interdisciplinary Theoretical & Mathematical Science Program (iTHEMS), RIKEN
 affiliation: Interdisciplinary Theoretical Science Research Group (iTHES), RIKEN
 slugcomment: submitted to ApJ
 http://www.pezy.co.jp/en/index.html
 http://www.aics.riken.jp/en/postk/project
References
 Abdo, A. A., Ackermann, M., Ajello, M., et al. 2009, ApJ, 700, 597
 Abdo, A. A., Ackermann, M., Ajello, M., et al. 2010a, ApJ, 710, 1271
 Abdo, A. A., Ackermann, M., Agudo, I., et al. 2010b , ApJ, 716, 30
 Abdo, A. A., Ackermann, M., Ajello, M., et al. 2011, ApJ, 733, L26
 Ackermann, M., Ajello, M., Baldini, L., et al. 2010, ApJ, 721, 1383
 Ackermann, M., Ajello, M., Ballet, J., et al. 2012, ApJ, 751, 159
 Ackermann, M., Anantua, R., Asano, K., et al. 2016, ApJ, 824, L20
 Asada, K., Nakamura, M., Doi, A., Nagai, H., & Inoue, M. 2014, ApJ, 781, L2
 Balbus, S. A., & Hawley, J. F. 1991, ApJ, 376, 214
 Begelman, M. C., Blandford, R. D., & Rees, M. J. 1984, Reviews of Modern Physics, 56, 255
 Biretta, J. A., Sparks, W. B., & Macchetto, F. 1999, ApJ, 520, 621
 Boccardi, B., Krichbaum, T. P., Bach, U., et al. 2016, A&A, 585, A33
 Bonnoli, G., Ghisellini, G., Foschini, L., Tavecchio, F., & Ghirlanda, G. 2011, MNRAS, 410, 368
 Britto, R. J., Bottacini, E., Lott, B., Razzaque, S., & Buson, S. 2016, ApJ, 830, 162
 Burgarella, D., Livio, M., & O’Dea, C. P. 1993, Astrophysical jets : proceedings of the Astrophysical Jets Meeting, Baltimore, 1992 May 1214 , 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.),
 Chandrasekhar, S. 1960, Proceedings of the National Academy of Science, 46, 253
 Chang, F.Y., Chen, P., Lin, G.L., Noble, R., & Sydora, R. 2009, Physical Review Letters, 102, 111101
 Chen, P., Tajima, T., & Takahashi, Y. 2002, Physical Review Letters, 89, 161101
 Chen, L.E., Li, H.Z., Yi, T.F., Zhou, S.B., & Li, K.Y. 2013, Research in Astronomy and Astrophysics, 13, 5
 Corde, S., Adli, E., Allen, J. M., et al. 2015, Nature, 524, 442
 Daniel, J., & Tajima, T. 1997, Phys. Rev. D, 55, 5193
 Daniel, J., & Tajima, T. 1998, ApJ, 498, 296
 Dermer, C. D., Razzaque, S., Finke, J. D., & Atoyan, A. 2009, New Journal of Physics, 11, 065016
 Duez, M. D., Liu, Y. T., Shapiro, S. L., Shibata, M., & Stephens, B. C. 2006, Phys. Rev. D, 73, 104015
 Ebisuzaki, T., & Tajima, T. 2014a, Astroparticle Physics, 56, 9
 Ebisuzaki, T., & Tajima, T. 2014b, European Physical Journal Special Topics, 223, 1113
 Esarey, E., Schroeder, C. B., & Leemans, W. P. 2009, Reviews of Modern Physics, 81,
 Fermi, E. 1954, ApJ, 119, 1
 Ferrari, A. 1998, ARA&A, 36, 539
 Fossati, G., Maraschi, L., Celotti, A., Comastri, A., & Ghisellini, G. 1998, MNRAS, 299, 433
 Gammie, C. F., McKinney, J. C., & Tóth, G. 2003, ApJ, 589, 444
 Gu, M., Cao, X., & Jiang, D. R. 2001, MNRAS, 327, 1111
 Kawamura, S., Nakamura, T., Ando, M., et al. 2006, Classical and Quantum Gravity, 23, S125
 Kino, M., Mizuta, A., & Yamada, S. 2004, ApJ, 611, 1021
 Kiuchi, K., Shibata, M., Montero, P. J., & Font, J. A. 2011, Physical Review Letters, 106, 251102
 Klein, A., Barausse, E., Sesana, A., et al. 2016, Phys. Rev. D, 93, 024003
 Kotera, K., & Olinto, A. V. 2011, ARA&A, 49, 119
 Hughes, P. A. (ed.) 1991, Beams and jets in astrophysics, Cambridge Univ. Press
 Lada, C. J. 1985, ARA&A, 23, 267
 Lau, C. K., Yeh, P. C., Luk, O., et al. 2015, Physical Review Special Topics Accelerators and Beams, 18, 024401
 Leemans, W. P., Nagler, B., Gonsalves, A. J., et al. 2006, Nature Physics, 2, 696
 Liu, F. K., Zhao, G., & Wu, X.B. 2006, ApJ, 650, 749
 Lynds, C. R. 1967, ApJ, 147, 837
 Machida, M., & Matsumoto, R. 2003, ApJ, 585, 429
 Matsubayashi, T., Shinkai, H.a., & Ebisuzaki, T. 2004, ApJ, 614, 864
 Matsumoto, R., & Tajima, T. 1995, ApJ, 445, 767
 Matsumoto, R., Uchida, Y., Hirose, S., et al. 1996, ApJ, 461, 115
 McKinney, J. C. 2006, MNRAS, 368, 1561
 McKinney, J. C., & Gammie, C. F. 2004, ApJ, 611, 977
 McKinney, J. C., Tchekhovskoy, A., & Blandford, R. D. 2012, MNRAS, 423, 3083
 Mima, K., Horton, W., Tajima, T., & Hasegawa, A. 1991, American Institute of Physics Conference Series, 230, 27
 Mimica, P., & Aloy, M. A. 2010, MNRAS, 401, 525
 Mineshige, S., Kusnose, M., & Matsumoto, R. 1995, ApJ, 445, L43
 Mirabel, I. F., & Rodríguez, L. F. 1999, ARA&A, 37, 409
 Misner, C. W., Thorne, K. S., & Wheeler, J. A. 1973, San Francisco: W.H. Freeman and Co., 1973,
 Nagataki, S. 2009, ApJ, 704, 937
 Nagataki, S. 2011, PASJ, 63, 1243
 Nakamura, K., Nagler, B., Tóth, C., et al. 2007, Physics of Plasmas, 14, 056708
 Nakamura, T., Ando, M., Kinugawa, T., et al. 2016, Progress of Theoretical and Experimental Physics, 2016, 093E01 p
 Nakano, H., Tanaka, T., & Nakamura, T. 2015, Phys. Rev. D, 92, 064003
 Narayan, R., & Yi, I. 1994, ApJ, 428, L13
 Noble, S. C., Gammie, C. F., McKinney, J. C., & Del Zanna, L. 2006, ApJ, 641, 626
 O’Neill, S. M., Reynolds, C. S., Miller, M. C., & Sorathia, K. A. 2011, ApJ, 736, 107
 Pe’er, A., Long, K., & Casella, P. 2016, arXiv:1610.08712
 Penna, R. F., McKinney, J. C., Narayan, R., et al. 2010, MNRAS, 408, 752
 Rees, M. J. 1978, MNRAS, 184, 61P
 Shakura, N. I., & Sunyaev, R. A. 1973, A&A, 24, 337
 Shibata, K., Tajima, T., & Matsumoto, R. 1990a, ApJ, 350, 295
 Shiokawa, H., Dolence, J. C., Gammie, C. F., & Noble, S. C. 2012, ApJ, 744, 187
 Spada, M., Ghisellini, G., Lazzati, D., & Celotti, A. 2001, MNRAS, 325, 1559
 Stone, J. M., Hawley, J. F., Gammie, C. F., & Balbus, S. A. 1996, ApJ, 463, 656
 Striani, E., Vercellone, S., Tavani, M., et al. 2010, ApJ, 718, 455
 Suzuki, T. K., & Inutsuka, S.i. 2009, ApJ, 691, L49
 Tajima, T. 2010, Proceeding of the Japan Academy, Series B, 86, 147
 Tajima, T., Nakajima, K., & Mourou, G. 2017, Nuovo Cimento Rivista Serie, 40, 33
 Tajima, T., & Dawson, J. M. 1979, Physical Review Letters, 43, 267
 Tajima, T., & Gilden, D. 1987, ApJ, 320, 741
 Takahashi Y., Hillman L. W, & Tajima T., 2000, in High Field Science, ed. Tajima T., Mima K., Baldis H. (Klewer), 171
 Tsiganos K. (ed.), 1996, Solar and Astrophysical Magnetohydrodynamic Flows, Dordrecht: Kluwer
 Yagi, K., & Seto, N. 2011, Phys. Rev. D, 83, 044011
 Velikhov E. P. 1959, ZhETF, 36, 1398