Bottom-up Acceleration of Ultra-High-Energy Cosmic Rays in the Jets of Active Galactic Nuclei
It has been proposed that Ultra-High-Energy Cosmic Rays (UHECRs) up to eV could be produced in the relativistic jets of powerful Active Galactic Nuclei (AGNs) via a one-shot reacceleration of lower-energy CRs produced in supernova remnants (the espresso mechanism). We test this theory by propagating particles in realistic 3D magneto-hydrodynamic simulations of ultra-relativistic jets and find that about 10% of the CRs entering the jet are boosted by at least a factor of in energy, where is the jet’s effective Lorentz factor, in agreement with the analytical expectations. Furthermore, about of the CRs undergo two or more shots and achieve boosts well in excess of . Particles are typically accelerated up to the Hillas limit, suggesting that the espresso mechanism may promote galactic-like CRs to UHECRs even in AGN jets with moderate Lorentz factors, and not in powerful blazars only. Finally, we find that the sign of the toroidal magnetic field in the jet and in the cocoon controls the angular distribution of the reaccelerated particles, leading to a UHECR emission that may be either quasi-isotropic or beamed along the jet axis. These findings strongly support the idea that espresso acceleration in AGN jets can account for the UHECR spectra, chemical composition, and arrival directions measured by Auger and Telescope Array.
The origin of the highest-energy cosmic rays (CRs) is one of the most prominent unresolved questions in astrophysics. Below eV, particles are thought to be accelerated in supernova remnants (SNRs) via diffusive shock acceleration (DSA) (e.g., Bell, 1978; Blandford & Ostriker, 1978; Berezhko & Völk, 2007; Ptuskin et al., 2010; Caprioli et al., 2010; Caprioli & Spitkovsky, 2014), consistent with observations of individual SNRs such as Tycho (Morlino & Caprioli, 2012; Slane et al., 2014), IC443, and W44 (Ackermann et al., 2013).
More precisely, data suggest that the CR flux has a cut-off in rigidity around eV (the CR “knee”), such that nuclei with atomic number are accelerated up to a maximum energy ; this corresponds to CR iron making up to eV and hence to a heavier and heavier composition above the knee (e.g., Bartoli et al., 2015; Hörandel, 2005; Dembinski et al., 2017).
As for Ultra-High-Energy Cosmic Rays (UHECRs) with energies between eV and eV, their sources and acceleration mechanism remain much less clear. The Pierre Auger Observatory has enabled a better understanding of the UHECR spectrum by showing evidence that the highest energy bins should contain nuclei heavier than hydrogen and helium (Aab et al., 2014a, b; Abbasi et al., 2015; Aab et al., 2017a). While at eV the composition is proton-only, at eV the UHECR composition is nitrogen-like, outlining a scenario where the UHECR composition becomes heavier as the energy increases, with a possible contribution from iron close to the highest-energy cut-off (Heinze et al., 2019). When statistics, limitations of current nuclear interaction models, and pipeline analyses are taken into account, this scenario is not at odds with Telescope Array data, which is also consistent with a lighter chemical composition on the whole energy range (Pierog, 2013; Abbasi et al., 2015).
1.1 Possible Acceleration Mechanisms
Based on energetics and luminosity arguments, -ray bursts (GRBs) (e.g., Vietri, 1995; Waxman, 1995), tidal disruption events (TDEs) (e.g., Farrar & Piran, 2014), newly-born millisecond pulsars (e.g., Blasi et al., 2007; Fang et al., 2012), and active galactic nuclei (AGNs) (e.g., Ostrowski, 2000; Murase et al., 2012) have been suggested as possible sources of particles up to eV. However, the actual mechanism(s) through which acceleration should proceed are not well delineated, in the sense that “models” very often are just back-of-the-envelope estimates of the maximum energy achievable in a given system. The Hillas criterion (Hillas, 1984; Cavallo, 1978)
expresses the minimum combination of size and magnetic field (in kpc and in G) necessary, but not sufficient, to accelerate a nucleus of charge up to energy (in units of eV) in a flow with speed . For relativistic flows, this criterion corresponds to a constraint on the particle Larmor radius, , and applies to both stochastic and one-shot acceleration mechanisms. In some cases, well-defined acceleration mechanisms (e.g., DSA, magnetic reconnection, shear acceleration,…) are also considered, but calculations require large extrapolations and parametrization of poorly-constrained ingredients such as particle injection and/or scattering.
DSA at non-relativistic shocks is a very robust acceleration process, but shocks on stellar (e.g., SNRs) and galactic scales (e.g., the wind termination shock), have hard times reaching the highest CR energies (e.g., Bell et al., 2013; Cardillo et al., 2015; Bustard et al., 2017); the interplay of multiple non-relativistic shocks in the backflowing material of AGN lobes (e.g., Matthews et al., 2019) and accretion shocks in galaxy clusters (e.g., Kang et al., 1996) may be more promising, though. DSA at relativistic shocks, which applies to GRBs, TDEs, and AGN internal shocks, has been shown to be generally less efficient and much slower than its non-relativistic counterpart (e.g., Sironi & Spitkovsky, 2011; Sironi et al., 2013; Araudo et al., 2018; Bell et al., 2018). Stochastic shear acceleration at the jet/cocoon interface in relativistic jets has also been suggested (e.g., Ostrowski, 2000; Kimura et al., 2018); any stochastic model, however, requires ad-hoc assumptions on spectrum and amplitude of the background magnetic turbulence, which are poorly constrained both observationally and theoretically. Furthermore, while newly-born millisecond pulsars can generate large voltages between poles and equator, it is not clear if/how particles could manage to cross magnetic field lines and tap the full potential drop.
Two additional general considerations. First, any acceleration mechanism that produces spectra or steeper is strongly disfavored if acceleration has to start from “thermal” particles, since the energetic constraint is already demanding at eV and a spectrum steeper than would have most of the power in low-energy particles. Second, the hint of a UHECR chemical composition significantly heavier than solar is at odds with leptonic environments such as pulsar magnetospheres (though see Fang et al., 2013) and poses tight constraints on sources and acceleration mechanisms. Kinetic simulations have recently shown that DSA at non-relativistic shocks naturally boosts the injection of heavy nuclei, in agreement with the elemental abundances of Galactic CRs (Caprioli et al., 2017); however, test-particle DSA only leads to spectra or steeper, and the standard non-linear theory of DSA (e.g., Jones & Ellison, 1991; Malkov & O’C. Drury, 2001), which predicts flatter spectra as a consequence of the back-reaction of accelerated particles, is at odds with the steep -ray spectra observed in SNRs (Caprioli, 2012, 2011). This list of possible UHECR sources and acceleration mechanisms is far from being comprehensive (for reviews see, e.g., Aharonian et al., 2002; Kotera & Olinto, 2011; Blandford et al., 2014), but it is safe to say that there is no consensus on where and how UHECRs are produced.
The goal of this paper is to build a solid theoretical framework where parameterizations are reduced to a minimum, if not eliminated, and UHECR acceleration is followed bottom-up from injection to the highest energies in an astrophysical source described in the most realistic way.
1.2 Espresso Reacceleration in AGNs
AGNs are excellent candidates as UHECR sources: an AGN jet with radius of hundreds of pc and of a few tens of G satisfies the Hillas criterion up to the highest energies for iron nuclei (equation 1). Also, AGN luminosities are consistent with the energy injection rate required to sustain the flux of UHECRs, erg Mpc yr (e.g., Katz et al., 2013), as extensively discussed, e.g., in Caprioli (2018) and references therein. In fact, assuming a typical density of AGNs Mpc, the luminosity of each AGN in UHECRs has to be
which is smaller or comparable to the bolometric luminosity of typical AGNs, erg s (e.g., Woo & Urry, 2002; Lusso et al., 2012). However, the ultimate source of energy that can be exploited to accelerate UHECRs is the jet power, which is a factor of 10–100 larger than (e.g., Ghisellini et al., 2009).
Caprioli (2015), hereafter C15, suggested that UHECRs may be produced in relativistic AGN jets through a very general mechanism dubbed espresso acceleration. The basic idea is that CR seeds accelerated up to eV in SNRs can penetrate into a relativistic jet and —independently of their exact trajectory— receive a one-shot boost of a factor of in energy, where is the Lorentz factor of the relativistic flow (the steam). With , as inferred from multi-wavelength observations of powerful blazars (e.g., Tavecchio et al., 2010; Zhang et al., 2014), a single espresso shot may be sufficient to boost the energy of galactic CRs by a factor , transforming the highest-energy galactic CRs at eV in the highest-energy UHECRs at eV.
The most appealing feature of the espresso mechanism is its simplicity: CR seeds undergo a Compton-like scattering with a relativistic “wall” (the jet magnetic field) and the resulting energy gain is of order if the initial and final directions of flight differ by more than . A close relative of this process is the first upstream–downstream–upstream cycle for a particle accelerated at a relativistic shock (e.g., Vietri, 1995; Achterberg et al., 2001). In the espresso framework, no assumptions are made on particle pitch-angle scattering (diffusion) or on the properties of the underlying magnetic turbulence, differently from the stochastic models that rely on repeated acceleration at the jet interface (e.g., Ostrowski, 1998, 2000; Kimura et al., 2018; Fang & Murase, 2018) or on multiple DSA in the jet cocoon (Matthews et al., 2019). In any stochastic model, the maximum energy critically depends on the rate at which CRs diffuse back to the acceleration sites, and in turn on the amplitude and spectrum of the magnetic turbulence at Larmor-radius scales, which is hard to constrain observationally.
It is worth estimating whether the energy budget in galactic CR seeds is consistent with a reacceleration scenario. We express the seed luminosity by considering the luminosity of the Milky Way in Galactic CRs, erg s (e.g., Hillas, 2005), and scaling it proportionally to the star formation rate , which for our Galaxy is a few solar masses per year; is a proxy for the rate of SN explosions and hence for the rate of production of CRs. Considering an injected spectrum , which contains the same energy per decade, allows us to work with the total seed luminosity rather than the luminosity in PeV particles.
Then, we estimate the fraction of the seeds that can be reprocessed by the jet. C15 (§3.2) gives the flux of CRs in the galactic halo that goes through the lateral surface of the jet, which reads
where is the minimum between the halo scale-height and the jet length, is the galaxy radius, and is the jet semi-aperture in degrees. With the reasonable assumptions of (it may be if self-confinement in the galactic halo is effective, see Blasi & Amato, 2019) and , we obtain .
Since the energy of the CR knee should not differ much from galaxy to galaxy (C15), galactic CRs need to gain a factor in energy to be promoted to UHECR; hence, the total luminosity in reaccelerated seeds can be estimated as
For our reference parameters, (equation 2) for a Milky-Way–like AGN host, and potentially much larger in starburst galaxies where .
C15 also pointed out that any reacceleration model that uses galactic CRs as seeds naturally predicts a match between the chemical composition at/above the knee and that of UHECRs: given a rigidity cutoff at a few PV, the UHECR spectrum should be proton-dominated at eV and increasingly heavier at higher energies, consistent with experimental data. This argument was put originally forward in the context of the espresso mechanism but has been applied to other frameworks, too (e.g., Kimura et al., 2018; Fang & Murase, 2018).
The espresso scenario has been corroborated with analytical calculations of CR trajectories in idealized jet structures, which confirmed that the vast majority of the trajectories lead to boosts regardless of the radial and longitudinal jet structures (Caprioli, 2018). Nevertheless, analytical calculations cannot answer some fundamental questions such as:
In a realistic jet, what is the fraction of CR seeds that can undergo espresso acceleration?
Do particles typically get a boost of ? Is it possible to undergo more than one shot and thus exceed such an estimate?
Is acceleration up to the Hillas limit generally achievable?
How does the spectrum of reaccelerated particles compare with the injected one?
Are reaccelerated particles released isotropically or are they beamed along the jet?
To address these points, we present here simulations in which test-particle CRs are propagated in more and more realistic jets; we start from idealized cylindrically-symmetric jets where the electromagnetic fields are prescribed analytically (§2) and culminate with self-consistent 3D magneto-hydrodynamic (MHD) simulations of ultra-relativistic jets (§3). Then, in §4 we discuss the energy spectrum and angular distribution of the accelerated particles. In §5 we conclude by summarizing the implications of our results for the origin of UHECRs in AGN jets.
2 particle trajectories and energy gain in cylindrical jets
The basic idea behind the espresso acceleration is that particles with Larmor radii large enough can penetrate into relativistic jets and experience the potential drop associated with the strong motional electric field, as seen from the laboratory frame. In C15 it was argued that particles on average gain a factor of in energy, where is the jet Lorentz factor, provided that the ingoing and outgoing flight directions are uncorrelated, which requires particles to perform at least a quarter of gyration in the flow before being released.
In this section, we use either an analytical Hamiltonian formalism (Caprioli, 2018) or a direct numerical approach to study particle trajectories and energy evolution in idealized jet structures; the main goal is to assess acceleration in non-homogeneous and finite jets.
Throughout the paper, we denote quantities in the laboratory and flow frames respectively with and , and initial/final quantities with the subscripts .
2.1 The Analytical Approach
We use cylindrical coordinates (,,) and consider a cylindrical jet with radius , with a magnetic field that is purely toroidal in the flow frame, , corresponding to a potential vector:
The flow has a velocity in the laboratory frame, where and are the speed of the flow and speed of light respectively; is the flow Lorentz factor. The potential vector transforms as , which means that is larger by a factor of in the laboratory frame.
We consider the Hamiltonian of a particle with mass , charge , and Lorentz factor in the flow frame:
where , , and are the canonical momenta. , , and are conserved quantities because is independent of , , and . Let us then consider a relativistic particle in the laboratory frame with initial energy and momentum
where and define the cosines of the flight angles with respect to and .
If we restrict ourselves to cases with , which correspond to and hence to planar orbits, the initial momentum is and the energy gain in the laboratory frame achieved by truncating the orbit at an arbitrary time reads (Caprioli, 2018):
is the average particle Larmor radius, , being the radially-averaged toroidal magnetic field. Equation 8 illustrates that the energy gain in the laboratory frame is due to the potential energy tapped by particles that penetrate into the jet: the closer they get to (where ), the larger the gain. For and in the limit , one has:
In general, the energy gain of a particle that travels back and forth between and some minimum oscillates between 1 and a value that depends on the initial Larmor radius. The average energy gain along the orbit, , obeys for any increasing profile of , and the maximum energy gain depends on the initial Larmor radius of the particle.
In calculating using equation 8, we distinguish the cases and , corresponding to final Larmor radii smaller or larger than the jet radius, in the flow frame. Note that is constant over the gyration, but it corresponds to different and in the laboratory frame.
: Particles get back to after one gyration; for planar orbits, the maximum energy is achieved at where and the momentum is entirely along the -direction. The maximum energy gain reads:
: Particles are not well magnetized and cross the entire flow; the maximum gain is achieved for and reads:
i.e., , which means that particles gain energy up to the Hillas limit.
The energy gain does not depend much on the pitch angle as long as it is not very close to 1. Finally, for orbits are non planar and have a finite constant ; in this case, particles do not reach because of the centrifugal barrier but the corrections to equation 11 and 12 are small, (Caprioli, 2018).
2.2 The Numerical Approach
In this section we numerically integrate orbits on top of a prescribed electromagnetic configuration for which an analytical solution cannot be easily found. Particles are propagated using the relativistic Boris algorithm (e.g., Birdsall & Langdon, 1991), which ensures long-term stability of the orbits; we have also tried the Vay pusher (Vay, 2008), which is known to perform better in relativistic flows, and it led to consistent results.
We consider the more general case of a cylindrical jet with and length . The Lorentz factor outside the jet is set up with a sharp exponential radial profile such that till it reaches for continuity purposes. We set a toroidal magnetic field inspired by the one generated by a current-carrying homogeneous wire, which goes as inside the jet and as outside, i.e.:
The magnetic field also decreases linearly along the jet, as suggested by radio observations (e.g., O’Sullivan & Gabuzda, 2009), which allows particles to eventually leave the jet.
Throughout the paper, we consider two possible orientations of the jet toroidal magnetic field, , which physically correspond to the jet current along being and call them case A and B, respectively. The sign of depends on the direction of the poloidal field in the material that is accreted on the black hole (e.g., Begelman et al., 1984), as recently shown also by first-principle PIC simulations (Parfrey et al., 2018), which in general does not have a preferential direction. One possible model that breaks such a symmetry is the “cosmic battery” mechanism (Contopoulos & Kazanas, 1998); in this picture radiative losses induce a drag between electrons and protons in the accretion disk, which generates a toroidal current and hence a fixed sign of the poloidal field that is accreted onto the black hole. The cosmic battery would lead to a counter-clockwise toroidal magnetic field (, , case A), provided that the accreted material does not have its own magnetic field in excess of the battery-generated one. In the absence of a definitive motivation to fix the signs of and , we leave both possibilities open and discuss how our results differ between case A and B.
Figure 1 shows the projection of the orbit on the plane (top panel) and the energy gain (bottom panel) for four particles with , , and different initial Larmor radii parameterized by , as in the legend; we also consider both jet polarizations, which lead to concave/convex trajectories. Espresso acceleration occurs for both toroidal magnetic field directions, as attested by the energy gain oscillating between 1 and during one orbit. Particles with reach a maximum energy gain , in agreement with equation 11, while particles with gain energy until their Larmor radius becomes comparable with (see equation 12). The dashed lines in the bottom panel of Figure 1 show the orbit-averaged energy gain , which is always a fraction of order one of the maximum energy gain .
To further test the effect of the magnetic field on the trajectories, we have added a poloidal component to the field, comparable in strength to the toroidal one, and found that the induced drift has no appreciable effect on the particle energy gain.
2.3 Non Homogeneous Jets
In general, relativistic jets do not have a uniform Lorentz factor as we have assumed so far. The interaction of the jet with the ambient medium creates regions that are pinch unstable (e.g., Hardee, 2000; Mignone et al., 2010; Tchekhovskoy & Bromberg, 2016), which leads to the formation of clumps where the Lorentz factor is larger separated by slower regions along the jet spine (e.g., Agudo et al., 2001; Hardcastle et al., 2016). This feature is commonly found also in high-resolution MHD simulations, so it is worth singling out the effects of such inhomogeneities on particle orbits.
We consider particles initialized with and different values of in a region of finite extent that has ; beyond , we set . The top panel of Figure 2 illustrates such a flow profile (greyscale), along with the trajectories of the propagated particles. Their maximum energy gain (bottom panel) is the expected for and increasingly smaller for larger Larmor radii. The critical value can be understood in the following way: in order to achieve the maximum energy boost a particle needs to stay in the high- flow for at least a quarter of a gyration, i.e., for a time , where is the gyration frequency in the flow frame. Since during this time the particle travels a distance along z, corresponds to the energy for which .
In summary, the conditions
express the constraints on the transverse () and longitudinal () sizes of a region with Lorentz factor that allow particles with initial Larmor radius to achieve the maximum theoretical energy boost . Both conditions can physically be interpreted as the Hillas criterion in the transverse and longitudinal directions, where in the longitudinal direction particles see a relativistically-contracted region. Note that both conditions correspond to requiring the final Larmor radius to be at most comparable with the size of the system.
3 Propagation in a full MHD simulation
In order to properly capture all the properties of a realistic astrophysical jet, we performed 3D MHD simulations with the PLUTO code (Mignone et al., 2012), which includes adaptive mesh refinement (AMR). A relativistic magnetized jet is launched in an unmagnetized uniform ambient medium, with a setup similar to the simulations presented in (Mignone et al., 2010). The jet is launched along the -direction through a cylindrical nozzle with a magnetization radius in a box that measures 48 in the and directions and 100 in the -direction. The grid has cells with 4 AMR levels, which allows us to resolve the instabilities at the jet–ambient medium interface in great detail. The characteristics of the system are specified by the jet/ambient density contrast , the launching Lorentz factor of the jet , and the jet sonic and Alfvénic Mach numbers and , where and are the sound and Alfvén speed, respectively; our fiducial parameters are , , 3, and . The jet is initialized with a purely toroidal magnetic field component such that:
We scale the magnetic field with respect to its initial value at the magnetization radius , which is defined by the chosen Alfvénic Mach number.
The simulations considered in this work cannot cover all the possible realizations of an AGN jet. For instance, varying the galactic density/temperature profile is known to affect the jet shape (e.g., Tchekhovskoy & Bromberg, 2016), and it is likely that varying the jet Lorentz factor, luminosity, and magnetization may lead to diverse jet morphologies. Nevertheless, our fiducial simulation and the analytical formalism of §2 allow us to characterize the general features of the acceleration mechanism in relativistic jets.
Figure 3 shows the volume rendering of the tracer variable for our benchmark run; represents the local fraction of jet material, such that corresponds to pure ambient medium and to pure jet medium. Jet and ambient material are well mixed, attesting that the Kelvin–Helmholtz instability in the cocoon and the wobbling of the jet lead to an effective mass entrainment (e.g., Mignone et al., 2010; Tchekhovskoy & Bromberg, 2016). This means that galactic-like CRs embedded in the ambient medium can be easily entrained, too, and be brought in contact with the most relativistic jet layers even without diffusing. Figure 4 shows 2D maps of and , which in the considered plane () correspond in modulus to the radial and toroidal components of the electric and magnetic field, respectively. These variables trace the morphology of the different jet regions (the spine, for , and the cocoon, for . We can see that the motional electric field reverses its sign between the two regions. This ultimately affects trajectories and fate of the reaccelerated particles, as we will explain more in detail in the following sections.
3.1 Particle Trajectories
For both case A and B, we propagate 100,000 test-particle protons with a broad range of initial Larmor radii and positions; the trajectories of nuclei with charge can be derived simply by considering protons with the same rigidity , whose Larmor radius is . We account for the two possible orientations of the toroidal magnetic field by flipping the sign of the propagated particles, so we can use the same MHD background and consider both case A () and B (). We initialize protons with normalized Larmor radii logarithmically spaced in the interval in order to cover an extended range of initial rigidities up to the Hillas limit, where particles should traverse the entire flow without significant acceleration. We notice that the actual strength of the magnetic field in the spine is larger than (see Figure 4): averaging over the regions where returns , so that the effective Hillas condition (equation 16) is met for .
Particles are initialized at linearly-spaced positions around the spine of the jet such that , and . The initial pitch angles are also linearly spaced with and (see definitions in §2).
Particle trajectories in the MHD jet show many features common with those discussed above for simplified jets. Figure 5 illustrates two examples of espresso-accelerated particles for case A and B (left and right panels, respectively). The top panels show the particle trajectories overplotted on a 2D slice of the component of the flow velocity, while the bottom panels show their energy gain as a function of ; the color code indicates the instantaneous Lorentz factor that they probe, .
The case-A particle (left panels) gains a factor of in energy during its first gyrations, where , and then encounters a major jet kink at , which fosters another energy boost, though with a smaller . Eventually, the particle loses some energy in crossing the cocoon and escapes the system with the canonical energy gain. The case-B particle (right panels) first gains a factor of , then undergoes another boost around , and is finally released with a total energy gain , well in excess of . Both particles gain energy up to the Hillas limit, though: the case-A particle has an initial Larmor radius and escapes the jet with , while the case-B particle is initialized with and escapes the jet with . It is important to stress that acceleration occurs when particles plunge into the relativistic flows, and neither at the interface between the spine and the cocoon, where the velocity shear is the largest (the boundary between the red and blue regions in Figure 5), nor in wake shocks.
These paradigmatic particle trajectories show that, regardless of the sign of the motional electric field: i) energization occurs because of espresso acceleration, and not because of stochastic or diffusive processes; ii) energy gains significantly larger than are possible and are favored by a non-uniform, turbulent jet that allows multiple acceleration cycles; iii) crossing the cocoon, which has an electric potential opposite to the jet’s, induces some energy losses that are generally smaller than the gain due to espresso acceleration; iv) particles tend to gain energy up to the Hillas limit.
3.2 Energy Gain Dependence on
It is crucial to point out that, even if the jet is launched with , most of the jet material moves with typical Lorentz factors smaller than . The histogram in Figure 6 shows the distribution of Lorentz factors in the relativistic flow (defined as regions where ). Only a few percent of the jet (in volume) moves with , most of the material having ; more precisely, the mean value of in the relativistic regions, which we use to define an effective , is . Regions with are rare and a consequence of the jet pinching.
Figure 6 also shows the average as a function of the largest Lorentz factor that particles probe along their trajectories. We can distinguish two regimes of acceleration depending on whether regions with a given have a large/small filling factor. On average, particles that experience regions with achieve the full energy boost, while the energy gain saturates to for particles that probe regions with larger Lorentz factors. The explanation for this trend hinges on the longitudinal condition in equation 16, which may be easily violated for the most relativistic regions. For instance, in our MHD simulation a typical particle with would need a region with that is as large as to be fully boosted, while in reality these regions are rare and/or much smaller.
In summary, particle energy gains are consistent with the espresso prediction of , where is the largest Lorentz factor that has a non-negligible filling factor; by the same token, should also be the characteristic Lorentz factor inferred from multi-wavelength observations of AGN jets.
3.3 Energy Gain Dependence on Initial Positions
Strictly speaking, CR seeds live in the ambient medium, but the vigorous mass entrainment effectively convects them inside the cocoon (Figure 3). On top of such a convective entrainment, seeds should diffuse through the ambient/jet interface thanks to both large-scale and microscopic turbulence. The role of large-scale turbulence is captured by the high resolution of our MHD simulation, while sub-grid magnetic fluctuations are not accounted for in the present work. In order for particles to be espresso accelerated, they need to reach the relativistic regions of the jet and the rate of percolation of seed particles into the jet spine may be enhanced by pitch-angle scattering due to small-scale fluctuations; therefore, our results can be viewed as a lower limit on the jet effectiveness in injecting seeds into the espresso mechanism.
Figure 7 illustrates the correlation between the particle initial position and final energy gain for seeds with , , and . The top panel shows the distribution of in a slice at , while the middle and bottom panels show the maps of for case A and B. In both cases particles typically gain at least a factor of in energy, but particles initialized in the highly-relativistic regions generally gain even more; the same trend is observed regardless of the initial parameters. The motivation for this correlation is that after the first acceleration shot the particle momentum is preferentially along (), which tends to reduce the energy gain of the following cycles (see equation 8); therefore, the possibility of having multiple shots relies on large-scale jet perturbations, in particular spine kinks and jet wobbling. Also, pitch-angle scattering due to small-scale turbulence generally helps breaking the correlation between in- and out-going angles, thereby fostering multiple acceleration cycles.
3.4 Energy Distribution of the Accelerated Particles
Figure 8 shows the spectrum of the escaping particles (solid black line), divided in the spectra produced by particles with a given initial (color lines). Four features are noteworthy: i) low-energy particles can be espresso-accelerated multiple times; for them energy gains are common; ii) the fraction of particles that are not reaccelerated increases at both extremes of the range; the optimal rigidity range for reacceleration is ; iii) the final spectrum is truncated at , a manifestation of the Hillas criterion; iv) the final spectrum tends to be flatter than the injected one because particles pile up close to the Hillas limit.
We have already discussed how seeds with do not gain much energy because they cannot complete one gyration in the flow; here we notice that for smaller and smaller the peaks in the color histograms in Figure 8 are more and more marked, attesting that a larger and larger fraction of the particles is not reaccelerated, the reason being that it is hard for low-energy particles to penetrate deep into the jet spine. The low-energy particles that manage to be reaccelerated often exhibit a double espresso shot with , but do not necessarily make it up to the Hillas limit. This may be due to the limited dynamical range achievable in 3D MHD simulations: a realistic AGN may show multiple kinks and accommodate even three or more shots and pile even more particles up at the Hillas limit. In general, we expect that the spectrum of reaccelerated particles should become rather flat at low energies, consistent with modeling the UHECR flux and composition (e.g., Gaisser et al., 2013; Aloisio et al., 2014; Taylor, 2014).
Multiple acceleration also shots ease the requirement of having for producing UHECRs starting from Galactic CRs that have a knee at a few PeV (C15), in the sense that the required energy gain could be achieved even in AGNs with smaller . Since the maximum energy gain scales as , where is the number of espresso cycles, even a moderate and 3-4 shots could lead to UHECR production. By the same token, multiple shots could provide and thus allow nitrogen and oxygen nuclei to reach the canonical eV, relaxing the need to have iron nuclei at the highest energies.
Finally, it is worth stressing that —even if multiple espresso shots were involved— the energization mechanism would still be distinct from stochastic shear acceleration or DSA in the jet cocoon (e.g., Kimura et al., 2018; Matthews et al., 2019), in that the energy gain is so effective that a few acceleration events are generally sufficient to reach the Hillas limit and no pitch-angle diffusion is required between two successive shots.
4 Espresso Acceleration in Astrophysical Jets
Let us now consider the results above in the context of typical AGN jets and for different species in seed CRs. In the MHD simulation, Larmor radii are normalized to the jet radius and the initial magnetic field; if we set pc and G, we can associate physical rigidities to the propagated particles, and the rigidity of the knee, GV, would correspond to . In this section we focus only on particles that have initial rigidities GV.
Following C15, we parametrize the energy flux of galactic CR seeds below the knee as:
where we consider the CR species [H, He, C/N/O, Fe] as grouped according to their effective atomic number and mass . The normalizations are chosen according to the abundance ratios at eV such that [1, 0.46, 0.30, 0.14]. The spectral slope observed at Earth is for protons and for heavier ions, a manifestation of the so-called discrepant hardening (Ahn et al., 2010; Caprioli et al., 2011). However, a conservation argument suggests that seeds in the galactic halo must have a spectrum parallel to the injection one, which should be significantly harder and closer to the universal spectrum produced at strong SNR shocks (, e.g., Blandford & Ostriker, 1978; Bell, 1978) or at most as suggested by -ray observations of young SNRs (Caprioli, 2011, 2012), anisotropy constraints (Blasi & Amato, 2012a, b), and secondary/primary abundances in Galactic CRs (e.g., Aguilar et al., 2016). A value of closer to the injection slope may also be realized in galaxies denser or bigger than the Milky Way, where spallation losses may dominate over diffusive escape (see the discussion in Caprioli, 2018). For these reasons, we initialize the seeds with the abundances above, , and introduce an abrupt rigidity cut-off at GV to facilitate the discrimination between seeds and reaccelerated particles.
4.1 Injection Efficiency
Figure 9 shows the cumulative distribution of the final energy gains of particles with GV, for both case A and B. If we consider the whole injection domain defined in §3.1 (upper curves), we find that () of case A(B) particles gain at least a factor of 2 in energy, and about () of case A(B) particles gain a factor of . Also, () of the particles achieve an energy gain of 100 or more in case A(B), corresponding to , i.e., two full espresso shots. The particles that do not gain much energy belong to the extremes of the range of initial rigidities: they either traverse the spine of the jet without gyrating because of their large Larmor radius, or had a Lorentz factor too small to enter the flow.
Including all the particles in the domain means assuming that particles can diffuse from ambient to jet medium, filling the domain in a uniform fashion. Instead, assuming that seeds are confined to the entrained ambient material only (regions with ) returns a lower limit on the fraction of seeds that can be reprocessed by the jet (bottom curves in Figure 9). In this case we find that () of case A(B) particles gain at least a factor of 2 in energy, and about () of case A(B) particles gain a factor of , as shown in Figure 9. We also find that of the particles in case B achieve an energy gain of 100, while particles in case A do not make it up to such a threshold. While it is indeed possible that the sign of the toroidal magnetic field may affect the fraction of particles that get multiple shots, in realistic systems diffusion is expected to be effective at some level. The fractions quoted above do not depend much on the chosen range of rigidity, but it is clear from Figure 8 that the injection efficiency would be quite larger if we only considered particles particles with . In summary, we conclude that a fraction of the seeds has to be reprocessed by at least one shot, and that about 0.05–0.1% can undergo two or more shots. Since a fraction is needed to account for the UHECR flux in the one-shot scenario (C15), we conclude that realistic AGN jets should be able to reaccelerate enough CR seeds to produce UHECRs.
4.2 The Spectrum of Released Particles
The espresso model predicts that the chemical composition observed in galactic CRs, which is increasingly heavy above eV and dominated by iron nuclei around eV (e.g., Hörandel et al., 2006; Kampert & Unger, 2012), should be mapped into UHECRs. Auger observations suggest a proton-dominated flux at eV and a heavier composition at higher energies, which is consistent with such a scenario. Figure 10 shows the energy spectrum produced in our fiducial MHD jet once galactic CRs are initialized as discussed above. We only consider particles that gained to eliminate reacclerated particles that should not contribute to the UHECR flux.
The cutoffs of seed species are boosted up in energy by a factor of , consistently with the results in §3, and spectra exhibit a high-energy tail because of the particles that underwent multiple acceleration cycles. Below the cutoffs spectra are still power-laws: parallel to the seed ones in case A, and slightly harder than the seed ones in case B (see also Figure 8). Modeling the UHECR flux and composition above eV (e.g., Gaisser et al., 2013; Aloisio et al., 2014; Taylor, 2014) favors rather hard injection spectra ( or flatter), which may seem hard to obtain with espresso reacceleration. Nevertheless, since particles tend to pile up close to the Hillas limit and since particles with are hardly reaccelerated, the spectrum of accelerated particles naturally flattens at low energies (see Figure 8), in agreement with the models above, which assume single power-law injection spectra.
Unfortunately the limited effective Lorentz factor and the intrinsic Hillas limit of our simulations do not allow particles to be accelerated beyond eV; since the proton spectrum cuts off quite close to the seed iron spectrum, we cannot directly observe the typical light-heavy-light-heavy modulation predicted in C15, either. In other words, the resolution of our fiducial MHD simulation forces us to initialize knee CRs with a Larmor radius that is too close to the jet scales for accommodating the energy gain required to produce eV particles. The extrapolation of our results to jets with larger Lorentz factors and/or with larger dynamical range between the seed rigidities and the jet extent is straightforward and is discussed in the following sections.
4.3 The Potential Role of Nuclei Photodisintegration
Photodisintegration inside or around the sources may destroy heavy ions and, in general, affect the spectra of escaping UHECRs in a different way for protons and other species. For ballistic propagation, even the intense radiation fields in the broad line region of the most luminous AGNs should not affect the UHECR spectrum appreciably (e.g., Dermer, 2007). Nevertheless, Unger et al. (2015) pointed out that —in the presence of magnetic irregularities that scatter and increase the UHECR confinement time in/around their sources — heavy ions may easily be photodisintegrated. They also argue that the products of such a disintegration naturally account for the light composition observed below eV, with the Galactic to extra-galactic transition populated of secondary protons, mainly. Fang & Murase (2018) recently suggested that a similar process may also account for the flux of high-energy neutrinos measured by Ice Cube and for the GeV -ray background measured by Fermi.
Relaxing the assumption that all the UHECR species are produced with the same slope also relaxes the requirement that the UHECR injection spectrum has to be quite flat (e.g., Gaisser et al., 2013; Aloisio et al., 2014; Taylor, 2014). In fact, it may be possible to reproduce the transition between Galactic and extra-galactic CRs starting with relatively-steep seed spectra (e.g., , with a low-rigidity cutoff for heavy nuclei and the contribution of secondary protons around eV. Note that this still does not explain the reason why the Galactic and extra-galactic components connect so smoothly, despite their different sources and transport properties. On the other hand, a single-source scenario for all of the CRs is strongly disfavored by the observed modulation in the chemical composition.
A detailed account for photodisintegration inside or in proximity of the jet is beyond the goal of this paper, but it is reasonable that heavy ions shot towards the black hole will be preferentially disintegrated by being exposed to the intense radiation fields in the AGN broad line region and/or in the dusty torus. In general, we see that the trajectories of lower energy UHECRs are far from being ballistic, especially for the lower-energy particles that propagate towards the jet base, so we expect photodisintegration to act as a high-pass filter that may make the UHECR injection spectrum flatter at low energy.
4.4 Dependence on the Jet Lorentz Factor
Although we initialized our fiducial jet with , the Lorentz factor that we would infer from a multi-wavelength analysis of its emission would be smaller (see Figure 6). Typical Lorentz factors inferred from multi-wavelength analysis of blazars are as large as 30–50, (Tavecchio et al., 2010; Zhang et al., 2014), but also more ordinary Seyfert galaxies may contain a similar ultra-relativistic jet spine despite being ordinarily associated with (e.g., Chiaberge et al., 2001). Very generally, when stratified jets are considered, there arises a unified scheme where Seyfert galaxies and blazars only differ by their viewing angle (e.g. Urry & Padovani, 1995; Ghisellini et al., 2005). Modelling hyper-relativistic flows in 3D MHD simulations is extremely challenging, so we resort to extrapolate our findings to realistic jets by scaling the measured energy gain with the jet Lorentz factor.
From figure 6 we see that the typical energy gain increases as up to and then levels off because the filling factor of regions with is small. By the same token, it is reasonable that also a multi-wavelength analysis would preferentially highlight regions with ; the most relativistic regions of the jet may still be the source of the jet fast variability and produce interesting relativistic effects (e.g., Giannios et al., 2009, and references therein). If the results above were extrapolated to the typical Lorentz factors inferred in powerful AGNs, i.e., for , about of the CR seeds should consistently achieve energy gains as large as (Figure 9), which would saturate the Hillas limit and produce the highest-energy CRs.
As discussed in C15, it is reasonable that the knee rigidity PeV does not depend much on the host galaxy. Assuming that the highest-energy galactic CRs are produced in SNRs around the Sedov stage, and that in the interstellar medium there is roughly equipartition between magnetic energy and thermal gas motions, one finds , where and are the number density and total mass of the galaxy considered. If any, in dense starburst galaxies the knee might be at energies a factor of a few larger because of such a scaling with , but also because of more effective CR diffusive shock reacceleration (Caprioli et al., 2018). Moreover, since the elemental composition is determined by the intrinsic dependence of ion injection into DSA (Caprioli et al., 2017), the elemental composition of seeds in other galaxies should resemble the one in the Milky Way, too. Finally, it is also possible that multiple episodes of AGN activity produced circumgalactic halos much more turbulent than the Galactic one, which results is a more effective confinement of both CRs accelerated in SNRs and low-energy reaccelerated ones. Such an effect may change slope and normalization of the seed spectrum, but it would systematically go in the direction of making it flatter and enhanced.
4.5 Angular Distribution of Reaccelerated Particles
A natural question is whether we should expect a correlation between the UHECR directions of arrival and the local AGN population; such a question is intimately connected to how particles are released from their sources, i.e., whether reaccelerated particles are strongly beamed along the jet or not. In the first case, UHECRs would preferentially come from AGNs that point at us, such as blazars or flat-spectra-radio quasars, while in the second case also more standard Seyfert galaxies may contribute, generally producing a more isotropic signal.
Intergalactic (and possibly galactic) magnetic fields may scramble the UHECR trajectories enough to break any correlation, but magnetic turbulence and hence scattering rates are quite uncertain in this respect. With our simulations we can however quantify if reaccelerated particles are preferentially released quasi isotropically or beamed along the jet axis.
Figure 11 shows how reaccelerated particles with are released. The vectors in the top panels illustrate the final directions of flight and the bottom panels show the distribution of , the cosine of the angle between the final particle velocity and the -axis; colors correspond to different energy bins. The lowest-energy particles (yellow arrows) tend to remain confined close to the jet spine, while higher-energy particles are released with significantly different angular distributions in case A and B. Case-A particles escape quasi isotropically, while case-B ones are strongly beamed along the jet axis. The very reason for this discrepancy hinges on the sign of the radial electric field in the cocoon, , which reverses its sign in the two cases. In fact, particles preferentially move along the -direction while in the spine, but their final direction is determined by the deflection they experience after having left the relativistic regions and entered the cocoon, where the flow is in the -direction (top panel of Figure 5). A different sign of in the jet produces a radial electric field that either disperses (, case A) or collimates (, case B) the flux of particles that escape the jet (Figure 4). Note that also in case B about half of the highest energy particles escape within an angle of , a bit larger than the canonical relativistic beaming of . Also considering that the sign of might as well change during the AGN lifetime, we conclude that the UHECR emission is not necessarily beamed along the jet as strongly as its -ray emission.
4.6 Correlation between UHECR Arrival Directions and Local AGNs
The results of §4.5 suggest that the jet cocoon is crucial for scattering escaping particles and potentially isotropizing them. More evolved jets with larger lobes and/or larger separation of scales between the relativistic spine and the cocoon may be even more effective in this respect. The typical size of AGN lobes (tens to hundreds of kpc) is in fact much larger than the particle Larmor radius even at the highest rigidities (kpc at V in a G field). Note that this does not mean that the Hillas limit is larger in the lobes, since there flows are non-relativistic and the maximum energy scales with the flow .
The astrophysical implication of a quasi-isotropic particle release is that we may receive UHECRs also from nearby AGNs whose jets do not point in our direction, i.e., non-blazar AGNs. For instance, Centaurus A is usually classified as a FR-I AGN with , despite likely hosting a very relativistic jet with (e.g., Chiaberge et al., 2001) and could be a source of UHECRs (e.g., Wykes et al., 2013).
It is still possible that a few powerful AGN may contribute in a more prominent way to the UHECR flux on top of a background that is quasi-isotropic, a noteworthy feature being the dipole measured by Auger (Aab et al., 2017b, 2018a). The only marginal ( post-trial) evidence for a hotspot has been reported by the Telescope Array Collaboration for events above eV Abbasi et al. (2014); such a hotspot correlates with the position of Mrk 421, one of the most powerful blazars in the local universe with a bolometric luminosity of and a luminosity distance of about Mpc (see Caprioli, 2018, for a more extended discussion). More recently, the Auger Collaboration reported on intermediate-scale anisotropy (Aab et al., 2018b) and analyzed the correlation of UHECR arrival directions with local starburst galaxies and GeV-bright AGNs, finding a correlation with the former and a correlation with the latter. These non-conclusive results are not at odds with the scenario outlined here because it is reasonable that starburst galaxies produce and retain more CR seeds (e.g., Yoast-Hull et al., 2014, 2015), but also because espresso acceleration may have occurred during previous AGN activities in the same galaxies. In fact, it is important to remember that UHECRs can self-confine themselves around their sources (Blasi et al., 2015) and be effectively released after the typical lifetime of an AGN jet, of the order of a few Myr. Finally, it is also possible that seed spectra in dense starburst galaxies may be flatter than in the Milky Way due to the balance of injection with spallation rather than diffusive escape (C15).
In summary, several different considerations suggest that, even if espresso acceleration in AGN jets were the main mechanism for generating UHECRs, their direction of arrival should not necessarily correlate with the most powerful nearby blazars: i) reaccelerated particles may be released almost isotropically, hence also non-blazar AGNs may contribute; ii) when just a few shots are allowed, even AGN jets with moderate produce UHECRs (note that the Hillas criterion does not depend on ); iii) self-confinement and propagation delay may offset the UHECR arrival with respect to the signature of a prominent jet activity.
We have extensively analyzed the espresso paradigm for the acceleration of UHECRs in relativistic AGN jets by propagating test-particles in both synthetic jet structures and full 3D relativistic MHD simulations. Our bottom-up approach accounts for all of the fundamental ingredients of a universal acceleration theory: it characterizes the conditions needed for injecting particles and follows their evolution in realistic environments up to the expected maximum energy; also, it does not include any sub-grid modeling or fine-tuning, and is consistent with the current UHECR phenomenology in terms of spectral slope, chemical composition, and anisotropy. To our knowledge, no other current model for the origin of UHECRs shares all of these properties.
The main results of our analysis corroborate the picture in which galactic CR seeds are reaccelerated through a one/two-shot process that taps the motional electric field in the ultra-relativistic jet spine. To be as general as possible, we consider two cases (A and B), where the initial toroidal magnetic field in the jet has different signs ( and , respectively).
By using analytical/numerical integration of particle trajectories in simplified jet structures (§2), we find that:
Particles are typically boosted by a factor of in energy regardless of the magnetic structure of the jet, in the sense that depends neither on the sign and the radial profile of toroidal magnetic field, nor on its poloidal component.
Achieving an energy gain in a jet region of transverse and longitudinal extent and requires the initial particle Larmor radius to satisfy (equation 16).
About 10% of the seeds are boosted by a factor of in energy and about of them gain more than , which correspond to one and two espresso shots, respectively (Figure 9). This efficiency is quite larger that the required to sustain the UHECR flux (C15).
The number of shots that a seed may undergo generally depends on the jet size and Lorentz factor distribution; nevertheless, since the energy gain scales as , a few shots are typically sufficient to reach the Hillas limit, i.e., UHECR energies, even in AGNs with moderate .
The sign of determines the sign of the (motional) radial electric field in the cocoon, which controls the angular distribution of the reaccelerated particles. In case A () particles are mainly released quasi isotropically, while in case B () particles are prefentially beamed along the jet axis (Figure 11).
Since reaccelerated particles may be released quasi-isotropically, and since also AGNs with apparent low Lorentz factors can produce UHECRs, the proposed model is consistent with the quasi-isotropic UHECR flux observed.
The possible correlation of the UHECR arrival direction with local starburst galaxies (Aab et al., 2018b) may be explained by the fact that CR seeds are overabundant in starburst galaxies. Hot-spots (Abbasi et al., 2014) may still appear in the direction of powerful blazar such as Mrk 421, especially for case-B jets.
In principle, espresso reacceleration may be at work also in other jetted astrophysical objects (microquasars, T-Tauri stars, TDEs, GRBs, kilonovae,…) and in relativistic flows (e.g., pulsar wind nebulae) but given the small spacial extent of these sources it is unlikely that they can reprocess enough CR seeds to contribute to the UHECR flux. On the other hand, espresso acceleration in these systems may produce a local population of energetic particles and be important for sculpting their non-thermal emission, especially if also electron reacceleration is considered.
In future works we plan to address in greater details the role of attenuation losses (especially photo-disintegration of the heavy nuclei shot towards the AGN, which may lead to a different spectra for protons and other ions) and the production of high-energy neutrinos in AGN jets.
- Aab et al. (2014a) Aab, A., Abreu, P., Aglietta, M., et al. 2014a, Phys. Rev. D, 90, 122005, doi: 10.1103/PhysRevD.90.122005
- Aab et al. (2014b) —. 2014b, Phys. Rev. D, 90, 122006, doi: 10.1103/PhysRevD.90.122006
- Aab et al. (2017a) Aab, A., Abreu, P., Aglietta, M., et al. 2017a, Phys. Rev. D, 96, 122003, doi: 10.1103/PhysRevD.96.122003
- Aab et al. (2017b) —. 2017b, Science, 357, 1266, doi: 10.1126/science.aan4338
- Aab et al. (2018a) —. 2018a, ApJ, 868, 4, doi: 10.3847/1538-4357/aae689
- Aab et al. (2018b) —. 2018b, ApJ, 853, L29, doi: 10.3847/2041-8213/aaa66d
- Abbasi et al. (2015) Abbasi, R., Bellido, J., Belz, J., et al. 2015, ArXiv e-prints. https://arxiv.org/abs/1503.07540
- Abbasi et al. (2014) Abbasi et al., R. U. 2014, ApJ Lett., 790, L21, doi: 10.1088/2041-8205/790/2/L21
- Achterberg et al. (2001) Achterberg, A., Gallant, Y. A., Kirk, J. G., & Guthmann, A. W. 2001, MNRAS, 328, 393, doi: 10.1046/j.1365-8711.2001.04851.x
- Ackermann et al. (2013) Ackermann et al., M. 2013, Science, 339, 807, doi: 10.1126/science.1231160
- Agudo et al. (2001) Agudo, I., Gómez, J.-L., Martí, J.-M., et al. 2001, ApJ, 549, L183, doi: 10.1086/319158
- Aguilar et al. (2016) Aguilar, M., Ali Cavasonza, L., Ambrosi, G., et al. 2016, Phys. Rev. Lett., 117, 231102, doi: 10.1103/PhysRevLett.117.231102
- Aharonian et al. (2002) Aharonian, F. A., Belyanin, A. A., Derishev, E. V., Kocharovsky, V. V., & Kocharovsky, V. V. 2002, Phys. Rev. D, 66, 023005, doi: 10.1103/PhysRevD.66.023005
- Ahn et al. (2010) Ahn, H. S., Allison, P., Bagliesi, M. G., et al. 2010, The Astrophysical Journal, 714, L89, doi: 10.1088/2041-8205/714/1/l89
- Aloisio et al. (2014) Aloisio, R., Berezinsky, V., & Blasi, P. 2014, J. Cosmology Astropart. Phys, 10, 20, doi: 10.1088/1475-7516/2014/10/020
- Araudo et al. (2018) Araudo, A. T., Bell, A. R., Blundell, K. M., & Matthews, J. H. 2018, MNRAS, 473, 3500, doi: 10.1093/mnras/stx2552
- Bartoli et al. (2015) Bartoli, B., Bernardini, P., Bi, X. J., et al. 2015, Phys. Rev. D, 92, 092005, doi: 10.1103/PhysRevD.92.092005
- Begelman et al. (1984) Begelman, M. C., Blandford, R. D., & Rees, M. J. 1984, Reviews of Modern Physics, 56, 255, doi: 10.1103/RevModPhys.56.255
- Bell (1978) Bell, A. R. 1978, MNRAS, 182, 147
- Bell et al. (2018) Bell, A. R., Araudo, A. T., Matthews, J. H., & Blundell, K. M. 2018, MNRAS, 473, 2364, doi: 10.1093/mnras/stx2485
- Bell et al. (2013) Bell, A. R., Schure, K. M., Reville, B., & Giacinti, G. 2013, MNRAS, 431, 415, doi: 10.1093/mnras/stt179
- Berezhko & Völk (2007) Berezhko, E. G., & Völk, H. J. 2007, ApJ, 661, L175, doi: 10.1086/518737
- Birdsall & Langdon (1991) Birdsall, C. K., & Langdon, A. B. 1991, Plasma Physics via Computer Simulation (CRC Press). http://adsabs.harvard.edu/abs/1991ppcs.book.....B
- Blandford et al. (2014) Blandford, R., Simeon, P., & Yuan, Y. 2014, Nuclear Physics B Proceedings Supplements, 256, 9, doi: 10.1016/j.nuclphysbps.2014.10.002
- Blandford & Ostriker (1978) Blandford, R. D., & Ostriker, J. P. 1978, ApJL, 221, L29, doi: 10.1086/182658
- Blasi & Amato (2012a) Blasi, P., & Amato, E. 2012a, J. Cosmology Astropart. Phys, 1, 10, doi: 10.1088/1475-7516/2012/01/010
- Blasi & Amato (2012b) —. 2012b, J. Cosmology Astropart. Phys, 1, 11, doi: 10.1088/1475-7516/2012/01/011
- Blasi & Amato (2019) Blasi, P., & Amato, E. 2019, arXiv e-prints. https://arxiv.org/abs/1901.03609
- Blasi et al. (2007) Blasi, P., Amato, E., & Caprioli, D. 2007, Monthly Notices of the Royal Astronomical Society, 375, 1471, doi: 10.1111/j.1365-2966.2006.11412.x
- Blasi et al. (2015) Blasi, P., Amato, E., & D’Angelo, M. 2015, Physical Review Letters, 115, 121101, doi: 10.1103/PhysRevLett.115.121101
- Bustard et al. (2017) Bustard, C., Zweibel, E. G., & Cotter, C. 2017, ApJ, 835, 72, doi: 10.3847/1538-4357/835/1/72
- Caprioli (2011) Caprioli, D. 2011, J. Cosmology Astropart. Phys, 5, 26, doi: 10.1088/1475-7516/2011/05/026
- Caprioli (2012) —. 2012, J. Cosmology Astropart. Phys, 7, 38, doi: 10.1088/1475-7516/2012/07/038
- Caprioli (2015) —. 2015, ApJ, 811, L38, doi: 10.1088/2041-8205/811/2/L38
- Caprioli (2018) Caprioli, D. 2018, NPPP. https://arxiv.org/abs/1705.00017
- Caprioli et al. (2010) Caprioli, D., Amato, E., & Blasi, P. 2010, APh, 33, 160, doi: 10.1016/j.astropartphys.2010.01.002
- Caprioli et al. (2011) Caprioli, D., Blasi, P., & Amato, E. 2011, APh, 34, 447, doi: 10.1016/j.astropartphys.2010.10.011
- Caprioli & Spitkovsky (2014) Caprioli, D., & Spitkovsky, A. 2014, ApJ, 783, 91, doi: 10.1088/0004-637X/783/2/91
- Caprioli et al. (2017) Caprioli, D., Yi, D. T., & Spitkovsky, A. 2017, Physical Review Letters, 119, doi: 10.1103/PhysRevLett.119.171101
- Caprioli et al. (2018) Caprioli, D., Zhang, H., & Spitkovsky, A. 2018, JPP. https://arxiv.org/abs/1801.01510
- Cardillo et al. (2015) Cardillo, M., Amato, E., & Blasi, P. 2015, APh, 69, 1, doi: 10.1016/j.astropartphys.2015.03.002
- Cavallo (1978) Cavallo, G. 1978, A&A, 65, 415
- Chiaberge et al. (2001) Chiaberge, M., Capetti, A., & Celotti, A. 2001, MNRAS, 324, L33, doi: 10.1046/j.1365-8711.2001.04642.x
- Contopoulos & Kazanas (1998) Contopoulos, I., & Kazanas, D. 1998, ApJ, 508, 859, doi: 10.1086/306426
- Dembinski et al. (2017) Dembinski, H., Engel, R., Fedynitch, A., et al. 2017, International Cosmic Ray Conference, 35, 533. https://arxiv.org/abs/1711.11432
- Dermer (2007) Dermer, C. D. 2007, ArXiv e-prints. https://arxiv.org/abs/0711.2804
- Fang et al. (2012) Fang, K., Kotera, K., & Olinto, A. V. 2012, ApJ, 750, 118, doi: 10.1088/0004-637X/750/2/118
- Fang et al. (2013) —. 2013, Journal of Cosmology and Astro-Particle Physics, 2013, 010, doi: 10.1088/1475-7516/2013/03/010
- Fang & Murase (2018) Fang, K., & Murase, K. 2018, Nature Physics, 14, 396, doi: 10.1038/s41567-017-0025-4
- Farrar & Piran (2014) Farrar, G. R., & Piran, T. 2014, arXiv e-prints. https://arxiv.org/abs/1411.0704
- Gaisser et al. (2013) Gaisser, T. K., Stanev, T., & Tilav, S. 2013, Frontiers of Physics, 8, 748, doi: 10.1007/s11467-013-0319-7
- Ghisellini et al. (2005) Ghisellini, G., Tavecchio, F., & Chiaberge, M. 2005, A&A, 432, 401, doi: 10.1051/0004-6361:20041404
- Ghisellini et al. (2009) Ghisellini, G., Tavecchio, F., & Ghirlanda, G. 2009, MNRAS, 399, 2041, doi: 10.1111/j.1365-2966.2009.15397.x
- Giannios et al. (2009) Giannios, D., Uzdensky, D. A., & Begelman, M. C. 2009, MNRAS, 395, L29, doi: 10.1111/j.1745-3933.2009.00635.x
- Hardcastle et al. (2016) Hardcastle, M. J., Lenc, E., Birkinshaw, M., et al. 2016, MNRAS, 455, 3526, doi: 10.1093/mnras/stv2553
- Hardee (2000) Hardee, P. E. 2000, ApJ, 533, 176, doi: 10.1086/308656
- Heinze et al. (2019) Heinze, J., Fedynitch, A., Boncioli, D., & Winter, W. 2019, arXiv e-prints. https://arxiv.org/abs/1901.03338
- Hillas (1984) Hillas, A. M. 1984, Ann. Rev. of A&A, 22, 425, doi: 10.1146/annurev.aa.22.090184.002233
- Hillas (2005) —. 2005, Journal of Physics G Nuclear Physics, 31, 95, doi: 10.1088/0954-3899/31/5/R02
- Hörandel (2005) Hörandel, J. R. 2005, astro-ph/0508014
- Hörandel et al. (2006) Hörandel et al., J. R. 2006, Journal of Physics Conference Series, 39, 463, doi: 10.1088/1742-6596/39/1/122
- Jones & Ellison (1991) Jones, F. C., & Ellison, D. C. 1991, Space Science Reviews, 58, 259
- Kampert & Unger (2012) Kampert, K.-H., & Unger, M. 2012, Astroparticle Physics, 35, 660, doi: 10.1016/j.astropartphys.2012.02.004
- Kang et al. (1996) Kang, H., Ryu, D., & Jones, T. W. 1996, ApJ, 456, 422, doi: 10.1086/176666
- Katz et al. (2013) Katz, B., Waxman, E., Thompson, T., & Loeb, A. 2013, ArXiv e-prints. https://arxiv.org/abs/1311.0287
- Kimura et al. (2018) Kimura, S. S., Murase, K., & Zhang, B. T. 2018, Phys. Rev. D, 97, 023026, doi: 10.1103/PhysRevD.97.023026
- Kotera & Olinto (2011) Kotera, K., & Olinto, A. V. 2011, ARA&A, 49, 119, doi: 10.1146/annurev-astro-081710-102620
- Lusso et al. (2012) Lusso, E., Comastri, A., Simmons, B. D., et al. 2012, MNRAS, 425, 623, doi: 10.1111/j.1365-2966.2012.21513.x
- Malkov & O’C. Drury (2001) Malkov, M. A., & O’C. Drury, L. 2001, Rep. Prog. Phys., 64, 429
- Matthews et al. (2019) Matthews, J. H., Bell, A. R., Blundell, K. M., & Araudo, A. T. 2019, MNRAS, 482, 4303, doi: 10.1093/mnras/sty2936
- Mignone et al. (2010) Mignone, A., Rossi, P., Bodo, G., Ferrari, A., & Massaglia, S. 2010, MNRAS, 402, 7, doi: 10.1111/j.1365-2966.2009.15642.x
- Mignone et al. (2012) Mignone, A., Zanni, C., Tzeferacos, P., et al. 2012, ApJS, 198, 7, doi: 10.1088/0067-0049/198/1/7
- Morlino & Caprioli (2012) Morlino, G., & Caprioli, D. 2012, A&A, 538, A81, doi: 10.1051/0004-6361/201117855
- Murase et al. (2012) Murase, K., Dermer, C. D., Takami, H., & Migliori, G. 2012, ApJ, 749, 63, doi: 10.1088/0004-637X/749/1/63
- Ostrowski (1998) Ostrowski, M. 1998, A&A, 335, 134
- Ostrowski (2000) —. 2000, MNRAS, 312, 579, doi: 10.1046/j.1365-8711.2000.03146.x
- O’Sullivan & Gabuzda (2009) O’Sullivan, S. P., & Gabuzda, D. C. 2009, MNRAS, 400, 26, doi: 10.1111/j.1365-2966.2009.15428.x
- Parfrey et al. (2018) Parfrey, K., Philippov, A., & Cerutti, B. 2018, arXiv e-prints. https://arxiv.org/abs/1810.03613
- Pierog (2013) Pierog, T. 2013, Journal of Physics Conference Series, 409, 012008, doi: 10.1088/1742-6596/409/1/012008
- Ptuskin et al. (2010) Ptuskin, V., Zirakashvili, V., & Seo, E.-S. 2010, ApJ, 718, 31, doi: 10.1088/0004-637X/718/1/31
- Sironi & Spitkovsky (2011) Sironi, L., & Spitkovsky, A. 2011, ApJ, 726, 75, doi: 10.1088/0004-637X/726/2/75
- Sironi et al. (2013) Sironi, L., Spitkovsky, A., & Arons, J. 2013, ApJ, 771, 54, doi: 10.1088/0004-637X/771/1/54
- Slane et al. (2014) Slane, P., Lee, S.-H., Ellison, D. C., et al. 2014, ApJ, 783, 33, doi: 10.1088/0004-637X/783/1/33
- Tavecchio et al. (2010) Tavecchio, F., Ghisellini, G., Ghirlanda, G., Foschini, L., & Maraschi, L. 2010, MNRAS, 401, 1570, doi: 10.1111/j.1365-2966.2009.15784.x
- Taylor (2014) Taylor, A. M. 2014, Astroparticle Physics, 54, 48 , doi: https://doi.org/10.1016/j.astropartphys.2013.11.006
- Tchekhovskoy & Bromberg (2016) Tchekhovskoy, A., & Bromberg, O. 2016, MNRAS, 461, L46, doi: 10.1093/mnrasl/slw064
- Unger et al. (2015) Unger, M., Farrar, G. R., & Anchordoqui, L. A. 2015, ArXiv e-prints. https://arxiv.org/abs/1505.02153
- Urry & Padovani (1995) Urry, C. M., & Padovani, P. 1995, PASP, 107, 803, doi: 10.1086/133630
- Vay (2008) Vay, J.-L. 2008, Physics of Plasmas, 15, 056701, doi: 10.1063/1.2837054
- Vietri (1995) Vietri, M. 1995, ApJ, 453, 883, doi: 10.1086/176448
- Waxman (1995) Waxman, E. 1995, Physical Review Letters, 75, 386, doi: 10.1103/PhysRevLett.75.386
- Woo & Urry (2002) Woo, J.-H., & Urry, C. M. 2002, ApJ, 579, 530, doi: 10.1086/342878
- Wykes et al. (2013) Wykes, S., Croston, J. H., Hardcastle, M. J., et al. 2013, A&A, 558, A19, doi: 10.1051/0004-6361/201321622
- Yoast-Hull et al. (2015) Yoast-Hull, T. M., Gallagher, J. S., & Zweibel, E. G. 2015, MNRAS, 453, 222, doi: 10.1093/mnras/stv1525
- Yoast-Hull et al. (2014) Yoast-Hull, T. M., Gallagher, III, J. S., Zweibel, E. G., & Everett, J. E. 2014, ApJ, 780, 137, doi: 10.1088/0004-637X/780/2/137
- Zhang et al. (2014) Zhang, B., Zhao, X., & Cao, Z. 2014, IJAA, 4, 499, doi: 10.4236/ijaa.2014.43046