Bottom-up Acceleration of Ultra-high-energy Cosmic Rays in the Jets of Active Galactic Nuclei

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; Matthews et al. 2019) have been suggested as possible sources of particles up to ∼10²⁰ 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 (Cavallo 1978; Hillas 1984),

$$\begin{eqnarray}&&{D}_{\mathrm{kpc}}{B}_{\mu {\rm{G}}}\beta \gtrsim \displaystyle \frac{100}{Z}\displaystyle \frac{E}{{10}^{20}\,\mathrm{eV}},\end{eqnarray} \tag{ 1 }$$

expresses the minimum combination of size D and magnetic field B (in kpc and in μG) necessary, but not sufficient, to accelerate a nucleus of charge Z up to energy E (in units of 10²⁰ eV) in a flow with speed βc. For relativistic flows, this criterion corresponds to a constraint on the particle Larmor radius, ${ \mathcal R }(E)\lesssim D$, 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 parameterization of poorly constrained ingredients such as particle injection and/or scattering.

DSA at nonrelativistic 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 nonrelativistic 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 nonrelativistic counterpart (e.g., Sironi & Spitkovsky 2011; Sironi et al. 2013; Araudo et al. 2018; Bell et al. 2018). Magnetic reconnection in newly born millisecond pulsars is another popular CR acceleration mechanism; pulsars can generate large voltages between poles and equator, but it is not clear if/how particles could manage to cross magnetic field lines and tap the full potential drop. Stochastic turbulent acceleration and shear acceleration at the jet/cocoon interface in relativistic jets have also been suggested (e.g., Hardcastle et al. 2009; O'Sullivan et al. 2009; Ostrowski 2000; Kimura et al. 2018).

Any acceleration mechanism that produces spectra E⁻² or steeper is strongly disfavored if acceleration has to start from "thermal" particles, since the energetic constraint is already demanding at E ≳ 10¹⁸ eV and a spectrum steeper than E⁻² would have most of the power in low-energy particles. Kinetic simulations have recently shown that DSA at nonrelativistic 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 ∝E⁻² or steeper, and the standard nonlinear theory of DSA (e.g., Jones & Ellison 1991; Malkov & Dury 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 2011, 2012). 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.

AGNs are excellent candidates as UHECR sources: an AGN jet with a radius of hundreds of parsecs and B 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, Q_UHECR ∼ 5 × 10⁴³ 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 n_AGN ≈ 10⁻⁴ Mpc⁻³, 10%–20% of which is radio-loud (e.g., Jiang et al. 2007), the luminosity of each AGN in UHECRs has to be

$$\begin{eqnarray}&&\tilde{{ \mathcal L }}\approx \displaystyle \frac{{Q}_{\mathrm{UHECR}}}{{n}_{\mathrm{AGN}}}\approx {10}^{40}\,\mathrm{erg}\,{{\rm{s}}}^{-1}.\end{eqnarray} \tag{ 2 }$$

Such a luminosity is smaller than the bolometric luminosity of typical AGNs, ${{ \mathcal L }}_{\mathrm{AGN}}^{\mathrm{bol}}\approx {10}^{42}\,-\,{10}^{45}$ erg s⁻¹ (e.g., Woo & Urry 2002; Lusso et al. 2012). Note that the upper end of this luminosity distribution is populated by powerful Fanaroff-Riley II (FR II) galaxies, which have number densities n_FRII ≈ 10⁻⁷ Mpc⁻³ and hence $\tilde{{ \mathcal L }}\approx {10}^{43}\,\mathrm{erg}\,{{\rm{s}}}^{-1}$. 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 ${{ \mathcal L }}_{\mathrm{AGN}}^{\mathrm{bol}}$ (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 10¹⁷ 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 Γ ≳ 30, as inferred from multiwavelength 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 Γ² ≳ 10³, transforming the highest-energy galactic CRs at 10¹⁷ eV in the highest-energy UHECRs at 10²⁰ 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 π/2. 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); the idea was also applied in blazar shocks assuming DSA and large-angle scattering, which in principle leads to multiple Γ² boosts (e.g., Stecker et al. 2007).

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; Fang & Murase 2018; Kimura et al. 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.

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 10¹⁸ eV and increasingly heavier at higher energies, consistent with experimental data. This argument was originally put forward in the context of the espresso mechanism but has been applied to other frameworks too (e.g., Fang & Murase 2018; Kimura et al. 2018).

Before diving into more detailed calculations, it is worth going through a simple estimate of the energetics of the galactic CR reacceleration. We express the seed luminosity by considering the luminosity of the Milky Way in galactic CRs, ${{ \mathcal L }}_{\mathrm{MW}}^{\mathrm{GCR}}\approx 5\times {10}^{40}$ erg s⁻¹ (e.g., Hillas 2005), and scaling it proportionally to the SN rate, ζ_MW ∼ 2 SNe century⁻¹ for our Galaxy. Considering an injected spectrum ∝E⁻², which contains the same energy per decade, allows us to work with the total seed luminosity rather than the luminosity in PeV particles; if the seed injection spectrum were steeper, ∝E^−2.3 as in the Milky Way (Blasi et al. 2012; Aguilar et al. 2016), the energy density in seeds at the knee would be about one order of magnitude smaller.

Then, we estimate the fraction of the seeds that can be reprocessed by the jet. C15 (Section 3.2) gives the flux of CRs in the galactic halo that goes through the lateral surface of the jet, which reads

$$\begin{eqnarray}&&{\rm{\Phi }}\approx 0.02{\rm{\Delta }}\theta {\left(\displaystyle \frac{H}{{R}_{\mathrm{gal}}}\right)}^{2},\end{eqnarray} \tag{ 3 }$$

where H is the minimum between the halo scale height and the jet length L_jet, R_gal is the galaxy radius, and Δθ is the jet semi-aperture in degrees. If the CR scale height is comparable to R_gal, we can put H ∼ L_jet; therefore, for L_jet/R_gal ≳ 1 and Δθ ≈ 2°, we obtain Φ ≳ 0.04. We comment on such assumptions in the next section.

Since the energy of the CR knee should not differ much from galaxy to galaxy (C15), galactic CRs need to gain a factor ${ \mathcal E }\approx {10}^{3}$ in energy to be promoted to UHECR; hence, the total luminosity in reaccelerated seeds can be estimated as

$$\begin{eqnarray}&&{{ \mathcal L }}_{\mathrm{AGN}}^{\mathrm{UHECR}}\approx { \mathcal E }{{ \mathcal L }}_{\mathrm{MW}}^{\mathrm{GCR}}\displaystyle \frac{\zeta {\rm{\Phi }}}{{\zeta }_{\mathrm{MW}}}\approx 2\times {10}^{42}\,\displaystyle \frac{\zeta }{{\zeta }_{\mathrm{MW}}}\displaystyle \frac{{\rm{\Phi }}}{0.04}\,\mathrm{erg}\,{{\rm{s}}}^{-1}.\end{eqnarray} \tag{ 4 }$$

For our reference parameters, ${{ \mathcal L }}_{\mathrm{AGN}}^{\mathrm{UHECR}}\gg \tilde{{ \mathcal L }}$ (Equation (2)) for a Milky Way-like AGN host. Due to the contribution of Type Ia SNe, the SN rate per unit of stellar mass is roughly the same in spiral and elliptical galaxies (e.g., Turatto et al. 1994); therefore, one can have ζ ≫ ζ_MW both in starburst and in massive elliptical galaxies, the latter being the typical AGN hosts. Moreover, it is plausible that the larger the SN/seed production rate, the stronger the interstellar magnetic turbulence and hence the longer the seed confinement time.

A natural question is what kinds of AGNs are the best candidates for espresso acceleration. Let us consider acceleration to ∼10²⁰eV, first. If such energies had to be achieved with one shot only, the jet Lorentz factor would need to be Γ ≳ 30, and such large Lorentz factors are typically inferred in blazars and radio-loud quasars (e.g., Tavecchio et al. 2010; Zhang et al. 2014). If a few espresso cycles were allowed (say, N), significantly lower values of Γ would suffice because the energy gain scales as Γ^2N. Since AGN jets typically have bulk flows Γ ≳ 5, and possibly even spines with Γ ≳ 10 (Chiaberge et al. 2001; Ghisellini et al. 2005; Lister et al. 2019), even ordinary Seyfert galaxies could in principle accelerate UHECRs.

In terms of required UHECR luminosity per single AGN (Equation (4)), it is likely that radio-quiet AGNs may not be powerful enough to contribute substantially to the UHECR flux (e.g., Kimura et al. 2018). When considering radio-loud AGNs, one has to distinguish between FR I jets, which are typically decelerated to nonrelativistic bulk flows within 1 kpc (e.g., Wardle & Aaron 1997; Arshakian & Longair 2004; Mullin & Hardcastle 2009), and FR II jets, which rather show Γ ≳ 10 at kiloparsec scales and beyond (e.g., Sambruna et al. 2002; Siemiginowska et al. 2002; Tavecchio et al. 2004; Harris & Krawczynski 2006).

FR I jets may have H ≲ 0.1R_gal and reprocess a fraction of the galactic seeds smaller than FR II jets (Equation (3)) but still be able to satisfy the condition in Equation (4). In such galaxies, the jet Lorentz factor may be of order of a few, and UHECR production would critically rely on the cocoon magnetic turbulence to be strong enough to allow a few acceleration cycles; for Γ as low as ≳3, N ≳ 3 espresso shots would be sufficient to promote galactic CRs to UHECRs. Instead, if the jet velocity becomes subrelativistic, espresso acceleration morphs into stochastic acceleration (Fang & Murase 2018; Kimura et al. 2018). For instance, N ≳ 12 for a jet bulk flow of ≃0.5c and Γ ≃ 1.33 (see Equation (10) of Caprioli 2018).

On the other hand, radio-loud FR II AGNs, despite being quite rare, are both powerful and extended enough to easily satisfy the energetic constraints for production of UHECRs via reacceleration of galactic CR seeds. Also, the large Lorentz factors persisting over kiloparsec scales provide the ideal conditions for espresso acceleration to occur, even via one/two shots only.

One final note: the boost of 10³ in energy and the corresponding requirements on Γ are based on the assumption that the knee energy of galactic CRs is the same in every galaxy (Caprioli 2015); it is indeed possible that the knee energy may be a factor of a few to 10 larger in some AGN hosts (e.g., in galaxies with prominent winds; see, e.g., Bustard et al. 2017), which would significantly reduce the requirements on the jet Lorentz factors in such environments.

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:

• 1.  
  In a realistic jet, what is the fraction of CR seeds that can undergo espresso acceleration?
• 2.  
  Do particles typically get a boost of Γ²? Is it possible to undergo more than one shot and thus exceed such an estimate?
• 3.  
  Is acceleration up to the Hillas limit generally achievable?
• 4.  
  How does the spectrum of reaccelerated particles compare with the injected one?
• 5.  
  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 (Section 2) and culminate with self-consistent 3D magnetohydrodynamic (MHD) simulations of ultrarelativistic jets (Section 3). Then, in Section 4 we discuss the energy spectrum and angular distribution of the accelerated particles. We conclude by summarizing the implications of our results for the origin of UHECRs in AGN jets (Section 5).

We use cylindrical coordinates (r, ϕ, z) and consider a cylindrical jet with radius ${R}_{\mathrm{jet}}$, with a magnetic field that is purely toroidal in the flow frame, ${\boldsymbol{B}}^{\prime} =-B^{\prime} (r){\boldsymbol{\phi }}$, corresponding to a potential vector:

$$\begin{eqnarray}&&{\boldsymbol{A}}^{\prime} ({\boldsymbol{r}})=A^{\prime} (r){\boldsymbol{z}},\quad \mathrm{with}\quad A^{\prime} (r)=-{\int }_{0}^{r}B^{\prime} (r^{\prime} ){dr}^{\prime} .\end{eqnarray} \tag{ 5 }$$

The flow has a velocity $\beta {\boldsymbol{z}}=v/c{\boldsymbol{z}}$ in the laboratory frame, where v and c are the speed of the flow and speed of light, respectively; Γ ≡ (1 − β²)^−1/2 is the flow Lorentz factor. The potential vector transforms as ${\boldsymbol{A}}={\rm{\Gamma }}{\boldsymbol{A}}^{\prime}$, which means that B_ϕ is larger by a factor of Γ in the laboratory frame.

We consider the Hamiltonian of a particle with mass m, charge q, and Lorentz factor γ' in the flow frame:

$$\begin{eqnarray}&&{ \mathcal H }^{\prime} =\sqrt{{P}_{r}{{\prime} }^{2}+\displaystyle \frac{{P}_{\phi }{{\prime} }^{2}}{{r}^{2}}+{[{P}_{z}^{\prime} -{qA}^{\prime} (r)]}^{2}+m}=\gamma ^{\prime} m,\end{eqnarray} \tag{ 6 }$$

where ${P}_{r}^{\prime} =\gamma ^{\prime} m\dot{r}$, ${P}_{\phi }^{\prime} =\gamma ^{\prime} {{mr}}^{2}\dot{\phi }$, and ${P}_{z}^{\prime} =\gamma ^{\prime} m\dot{z}+{qA}^{\prime}$ are the canonical momenta. ${ \mathcal H }^{\prime}$, P_ϕ', and P_z' are conserved quantities because ${ \mathcal H }^{\prime}$ is independent of t, ϕ, and z. Let us then consider a relativistic particle in the laboratory frame with initial energy E_i and momentum

$$\begin{eqnarray}&&{{\boldsymbol{p}}}_{{\rm{i}}}\simeq {E}_{{\rm{i}}}(-{\xi }_{{\rm{i}}}\sqrt{1-{\mu }_{{\rm{i}}}^{2}};\sqrt{(1-{\mu }_{{\rm{i}}}^{2})(1-{\xi }_{{\rm{i}}}^{2})};{\mu }_{{\rm{i}}}),\end{eqnarray} \tag{ 7 }$$

where ${\mu }_{{\rm{i}}}\equiv {p}_{z}/| {\boldsymbol{p}}|$ and ${\xi }_{{\rm{i}}}\equiv {p}_{r}/| {\boldsymbol{p}}|$ define the cosines of the flight angles with respect to ${\boldsymbol{z}}$ and $-{\boldsymbol{r}}$.

If we restrict ourselves to cases with ${\xi }_{{\rm{i}}}$ = 1, which correspond to P_ϕ'= 0 and hence to planar orbits, the initial momentum is ${{\boldsymbol{p}}}_{{\rm{i}}}\simeq {E}_{{\rm{i}}}(-\sqrt{1-{\mu }_{{\rm{i}}}^{2}};0;{\mu }_{{\rm{i}}})$ and the energy gain in the laboratory frame achieved by truncating the orbit at an arbitrary time t_f reads as (Caprioli 2018)

$$\begin{eqnarray}\begin{array}{rcl}{ \mathcal E } & \equiv & \displaystyle \frac{{E}_{{\rm{f}}}}{{E}_{{\rm{i}}}}={{\rm{\Gamma }}}^{2}(1-\beta {\mu }_{{\rm{i}}})(1+\beta {\mu }_{{\rm{f}}}^{{\prime} })\\ & = & (1-\beta {\mu }_{{\rm{i}}})\left\{1+\displaystyle \frac{{{\rm{\Gamma }}}^{2}\beta }{\alpha ^{\prime} }\left[1-\displaystyle \frac{A(r)}{A({R}_{\mathrm{jet}})}\right]\right\},\end{array}\end{eqnarray} \tag{ 8 }$$

where

$$\begin{eqnarray}&&\alpha ^{\prime} \equiv \displaystyle \frac{{ \mathcal R }^{\prime} }{{R}_{\mathrm{jet}}}=\displaystyle \frac{E^{\prime} }{q\langle B^{\prime} {\rangle }_{r}{R}_{\mathrm{jet}}}={{\rm{\Gamma }}}^{2}{\alpha }_{{\rm{i}}}={\alpha }_{{\rm{f}}}.\end{eqnarray} \tag{ 9 }$$

${ \mathcal R }^{\prime}$ is the average particle Larmor radius, $A^{\prime} (r)\equiv r\langle B^{\prime} {\rangle }_{r}$, $\langle B^{\prime} {\rangle }_{r}$ 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 r = 0 (where A' = 0), the larger the gain. For ${\mu }_{{\rm{i}}}\ne 1$ and in the limit β → 1, Γ² ≫ α' one has

$$\begin{eqnarray}&&{ \mathcal E }\approx \displaystyle \frac{{{\rm{\Gamma }}}^{2}}{\alpha ^{\prime} }\left[1-\displaystyle \frac{A(r)}{A({R}_{\mathrm{jet}})}\right].\end{eqnarray} \tag{ 10 }$$

In general, the energy gain of a particle that travels back and forth between ${R}_{\mathrm{jet}}$ and some minimum r oscillates between 1 and a value ${{ \mathcal E }}_{\max }$ that depends on the initial Larmor radius. The average energy gain along the orbit, $\langle { \mathcal E }\rangle$, obeys ${{ \mathcal E }}_{\max }\geqslant \langle { \mathcal E }\rangle \geqslant {{ \mathcal E }}_{\max }/2$ for any increasing profile of B(r), and the maximum energy gain depends on the initial Larmor radius of the particle.

In calculating ${{ \mathcal E }}_{\max }$ using Equation (8), we distinguish the cases α' > 1/2 and α' < 1/2, 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 α_i and α_f in the laboratory frame.

• 1.  
  α' < 1/2: Particles get back to r = ${R}_{\mathrm{jet}}$ after one gyration; for planar orbits, the maximum energy is achieved at r = r_min, where P_r = 0 and the momentum is entirely along the z-direction. The maximum energy gain reads as

  $$\begin{eqnarray}&&{{ \mathcal E }}_{\max }={{\rm{\Gamma }}}^{2}(1-\beta {\mu }_{{\rm{i}}})(1+\beta )\simeq 2{{\rm{\Gamma }}}^{2}.\end{eqnarray} \tag{ 11 }$$
• 2.  
  α' > 1/2: Particles are not well magnetized and cross the entire flow; the maximum gain is achieved for A (r = 0) = 0 and reads as

  $$\begin{eqnarray}&&{{ \mathcal E }}_{\max }=(1-\beta {\mu }_{{\rm{i}}})\left[1+\displaystyle \frac{{{\rm{\Gamma }}}^{2}\beta }{\alpha ^{\prime} }\right]\simeq \displaystyle \frac{{{\rm{\Gamma }}}^{2}}{\alpha ^{\prime} }=\displaystyle \frac{1}{{\alpha }_{{\rm{i}}}},\end{eqnarray} \tag{ 12 }$$

  i.e., α_f ≃ 1, which means that particles gain energy up to the Hillas limit.

The energy gain does not depend much on the pitch angle ${\mu }_{{\rm{i}}}$ as long as it is not very close to 1. Finally, for ${\xi }_{{\rm{i}}}$ < 1 orbits are nonplanar and have a finite constant P_ϕ; in this case, particles do not reach r = 0 because of the centrifugal barrier, but the corrections to Equation (11) and (12) are small, ${ \mathcal O }(1/{{\rm{\Gamma }}}^{2})$ (Caprioli 2018).

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 Γ₀ = 10 and length L_jet = 200${R}_{\mathrm{jet}}$. The Lorentz factor outside the jet is set up with a sharp exponential radial profile such that ${\rm{\Gamma }}(r\gt {R}_{\mathrm{jet}})={{\rm{\Gamma }}}_{0}{e}^{-10(r/{R}_{\mathrm{jet}}-1)}$ until it reaches Γ = 1 for continuity purposes. We set a toroidal magnetic field inspired by the one generated by a current-carrying homogeneous wire, which goes as r inside the jet and as 1/r outside, i.e.,

$$\begin{eqnarray}&&{\boldsymbol{B}}(r\leqslant {R}_{\mathrm{jet}},z)=B({R}_{\mathrm{jet}},z)\displaystyle \frac{r}{{R}_{\mathrm{jet}}}{\boldsymbol{\phi }},\end{eqnarray} \tag{ 13 }$$

$$\begin{eqnarray}&&{\boldsymbol{B}}(r\geqslant {R}_{\mathrm{jet}},z)=B({R}_{\mathrm{jet}},z)\displaystyle \frac{{R}_{\mathrm{jet}}}{r}{\boldsymbol{\phi }}\end{eqnarray} \tag{ 14 }$$

$$\begin{eqnarray}&&B({R}_{\mathrm{jet}},z)={B}_{0}\displaystyle \frac{{L}_{\mathrm{jet}}-z}{{L}_{\mathrm{jet}}}.\end{eqnarray} \tag{ 15 }$$

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, ${B}_{\phi }\lessgtr 0$, which physically correspond to the jet current along z being ${J}_{z}\lessgtr 0$, and call them case A and case B, respectively. The sign of B_ϕ 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 particle-in-cell simulations (Parfrey et al. 2019), 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 counterclockwise toroidal magnetic field (J_z < 0, B_ϕ < 0, 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 J_z and B_ϕ, we leave both possibilities open and discuss how our results differ between case A and case B.

Figure 1 shows the projection of the orbit on the x–z plane (top panel) and the energy gain ${ \mathcal E }$ (bottom panel) for four particles with μ_i = 0, ξ_i = 1, and different initial Larmor radii parameterized by $\alpha ^{\prime} ={ \mathcal R }^{\prime} /{R}_{\mathrm{jet}}$, 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 ${{ \mathcal E }}_{\max }$ during one orbit. Particles with α'< 1/2 reach a maximum energy gain ${{ \mathcal E }}_{\max }\sim 2{{\rm{\Gamma }}}^{2}$, in agreement with Equation (11), while particles with α'> 1/2 gain energy until their Larmor radius becomes comparable to ${R}_{\mathrm{jet}}$ (see Equation (12)). The dashed lines in the bottom panel of Figure 1 show the orbit-averaged energy gain $\langle { \mathcal E }\rangle$, which is always a fraction of order one of the maximum energy gain ${{ \mathcal E }}_{\max }$.

Zoom In Zoom Out Reset image size

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 ${\boldsymbol{B}}\times {\boldsymbol{\nabla }}B$ drift has no appreciable effect on the particle energy gain.

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 ${\mu }_{{\rm{i}}}$ = 0 and different values of $\alpha ^{\prime} \in [0.05,0.4]$ in a region of finite extent H = 2.5${R}_{\mathrm{jet}}$ that has Γ = 10; beyond x = H, we set Γ = 5. The top panel of Figure 2 illustrates such a flow profile (gray scale), along with the trajectories of the propagated particles. Their maximum energy gain (bottom panel) is the expected 2Γ² ≃ 200 for $\alpha ^{\prime} \lesssim \bar{\alpha }^{\prime} \simeq 0.16$ and increasingly smaller for larger Larmor radii. The critical value $\bar{\alpha }^{\prime}$ 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 one-quarter of a gyration, i.e., for a time T_acc ≳ Γπ/(2Ω'), where ${\rm{\Omega }}^{\prime} \equiv { \mathcal R }^{\prime} /c$ is the gyration frequency in the flow frame. Since, during this time, the particle travels a distance Δz = T_accc = πα' ${R}_{\mathrm{jet}}$Γ/2 along z, $\bar{\alpha }^{\prime}$ corresponds to the energy for which T_accc = H.

Zoom In Zoom Out Reset image size

In summary, the conditions

$$\begin{eqnarray}&&{{ \mathcal R }}_{{\rm{i}}}\lesssim \displaystyle \frac{{R}_{\mathrm{jet}}}{2{{\rm{\Gamma }}}^{2}};\quad {{ \mathcal R }}_{{\rm{i}}}\lesssim \displaystyle \frac{2}{\pi }\displaystyle \frac{H}{{{\rm{\Gamma }}}^{3}}\approx \displaystyle \frac{H}{{{\rm{\Gamma }}}^{3}}\end{eqnarray} \tag{ 16 }$$

express the constraints on the transverse (${R}_{\mathrm{jet}}$) and longitudinal (H) sizes of a region with Lorentz factor Γ that allow particles with initial Larmor radius ${{ \mathcal R }}_{{\rm{i}}}$ 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 ${{ \mathcal R }}_{{\rm{f}}}$ to be at most comparable to the size of the system.

For both case A and case B, we propagate ∼100,000 test-particle protons with a broad range of initial Larmor radii ${ \mathcal R }$ and positions; the trajectories of nuclei with charge q = Ze can be derived simply by considering protons with the same rigidity ρ ≡ E/q, whose Larmor radius is ${ \mathcal R }=\tfrac{E}{{ZeB}}=\tfrac{\rho }{B}$. 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 (B_ϕ < 0) and case B (B_ϕ > 0). We initialize protons with normalized Larmor radii logarithmically spaced in the interval ${\alpha }_{{\rm{i}}}=\tfrac{\rho }{{B}_{0}{R}_{\mathrm{jet}}}\in [{10}^{-3.6},8]$ 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 B₀ (see Figure 4): averaging over the regions where Γ ≥ 2 returns B_eff ∼ 7.2B₀, so that the effective Hillas condition (Equation (16)) is met for α_H ∼ 7.2.

Particles are initialized at linearly spaced positions $({r}_{{\rm{i}}},{\phi }_{{\rm{i}}},{z}_{{\rm{i}}})$ around the spine of the jet such that ${r}_{{\rm{i}}}/{R}_{\mathrm{jet}}\in [0.2,5]$, ${z}_{{\rm{i}}}/{R}_{\mathrm{jet}}\in [2,60]$, and ${\phi }_{{\rm{i}}}\in [0,2\pi ]$. The initial pitch angles are also linearly spaced with ${\mu }_{{\rm{i}}}\in [-0.5,0.5]$ and ${\xi }_{{\rm{i}}}\in [-1,1]$ (see definitions in Section 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 case B (left and right panels, respectively). The top panels show the particle trajectories overplotted on a 2D slice of the z component of the flow velocity, while the bottom panels show their energy gain ${ \mathcal E }$ as a function of z; the color code indicates the instantaneous Lorentz factor that they probe, Γ_pr.

Zoom In Zoom Out Reset image size

The case A particle (left panels) gains a factor of $\sim {{\rm{\Gamma }}}_{\mathrm{pr}}^{2}$ in energy during its first gyrations, where Γ_pr ≲ 3.6, and then encounters a major jet kink at z ≳ 54${R}_{\mathrm{jet}}$, 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 $\sim {{\rm{\Gamma }}}_{\mathrm{pr}}^{2}$ energy gain. The case B particle (right panels) first gains a factor of $\sim 2{{\rm{\Gamma }}}_{\mathrm{pr}}^{2}$, then undergoes another boost around 55–60 ${R}_{\mathrm{jet}}$, and is finally released with a total energy gain ${ \mathcal E }\sim 50$, well in excess of ${{\rm{\Gamma }}}_{\mathrm{pr}}^{2}$. Both particles gain energy up to the Hillas limit, though: the case A particle has an initial Larmor radius α_i = 0.42 and escapes the jet with α_f ≃ 5.9 ≈ α_H, while the case B particle is initialized with α_i = 0.13 and escapes the jet with α_f ≃ 6.34 ≈ α_H. 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 ${{\rm{\Gamma }}}_{\mathrm{pr}}^{2}$ are possible and are favored by a nonuniform, 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; and (iv) particles tend to gain energy up to the Hillas limit.

It is crucial to point out that, even if the jet is launched with Γ₀ = 7, 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 Γ ≥ 2). Only a few percent of the jet (in volume) moves with Γ ∼ Γ₀, most of the material having Γ ≲ 5; more precisely, the mean value of Γ in the relativistic regions, which we use to define an effective Γ_eff, is ∼3.2. Regions with Γ > 7 are rare and a consequence of the jet pinching.

Zoom In Zoom Out Reset image size

Figure 6 also shows the average ${{ \mathcal E }}_{\max }$ 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 Γ ≲ ${{\rm{\Gamma }}}_{\mathrm{eff}}$ achieve the full Γ² energy boost, while the energy gain saturates to ${ \mathcal E }\approx 15\mbox{--}20$ 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 α ∼ 0.01 would need a region with Γ ≳ 7 that is as large as ≳6${R}_{\mathrm{jet}}$ 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 ${ \mathcal E }\gtrsim {{\rm{\Gamma }}}_{\mathrm{eff}}^{2}$, where ${{\rm{\Gamma }}}_{\mathrm{eff}}$ is the largest Lorentz factor that has a nonnegligible filling factor; by the same token, ${{\rm{\Gamma }}}_{\mathrm{eff}}$ should also be the characteristic Lorentz factor inferred from multiwavelength observations of AGN jets.

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 subgrid 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 α_i = 0.2, ${\mu }_{{\rm{i}}}$ = 0, and ${\xi }_{{\rm{i}}}$ = 1. The top panel shows the distribution of Γ(x, z) in a slice at y = 0, while the middle and bottom panels show the maps of ${ \mathcal E }(x,z)$ for case A and case B, respectively. In both cases particles typically gain at least a factor of ∼${{\rm{\Gamma }}}_{\mathrm{eff}}$² 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 +z (μ ≈ 1), 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 in breaking the correlation between in- and outgoing angles, thereby fostering multiple acceleration cycles.

Zoom In Zoom Out Reset image size

Figure 8 shows the spectrum of the escaping particles (solid black line), divided into the spectra produced by particles with a given initial α_i (color lines). Four features are noteworthy: (i) low-energy particles can be espresso-accelerated multiple times; for them energy gains ${ \mathcal E }\gtrsim 100$ are common; (ii) the fraction of particles that are not reaccelerated increases at both extremes of the α_i range; the optimal rigidity range for reacceleration is ${\alpha }_{{\rm{i}}}\in \sim [0.01,1];$ (iii) the final spectrum is truncated at α_H ∼ 7.2, 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.

Zoom In Zoom Out Reset image size

We have already discussed how seeds with α ≳ 1 do not gain much energy because they cannot complete one gyration in the flow; here we notice that for smaller and smaller α_i the peaks in the color histograms in Figure 8 are more and more marked, attesting that a larger and larger fraction of the particles are 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 ${ \mathcal E }\gtrsim 100$ 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 shots also ease the requirement of having Γ = 30 for producing UHECRs starting from galactic CRs that have a knee at a few PeV (C15), in the sense that the required energy gain ${ \mathcal E }\sim {10}^{3}$ could be achieved even in AGNs with smaller ${{\rm{\Gamma }}}_{\mathrm{eff}}$. Since the maximum energy gain scales as ${{\rm{\Gamma }}}_{\mathrm{eff}}$^2N, where N is the number of espresso cycles, even a moderate ${{\rm{\Gamma }}}_{\mathrm{eff}}$ ≳ 3 and 3–4 shots could lead to UHECR production. By the same token, multiple shots could provide ${ \mathcal E }\gtrsim 3\times {10}^{3}$ and thus allow nitrogen and oxygen nuclei to reach the canonical 10²⁰ eV, relaxing the need to have iron nuclei at the highest energies.

It is worth discussing how the energization process that we observe is distinct from other processes suggested in the literature, such as stochastic shear acceleration, DSA in the cocoon, and turbulent acceleration (e.g., Ostrowski 2000; Hardcastle et al. 2009; O'Sullivan et al. 2009; Kimura et al. 2018; Matthews et al. 2019).

First of all, in our simulations all of the accelerated particles gain energy through one or at most two/three shots. Their trajectories closely resemble the ordered gyrations outlined in Section 2 (Figure 5), with no evidence of spatial random walk and pitch-angle scattering between two successive shots. Multiple shots are typically due to large-scale kinks in the jet, rather than to pitch-angle diffusion.

In stochastic shear acceleration, instead, particles are expected to be accelerated at the interface between the cocoon and the jet (Ostrowski 1998, 2000; Kimura et al. 2018) and to diffuse throughout the cocoon; the same random walk in the cocoon should also be present if acceleration were due to second-order Fermi processes (e.g., Hardcastle et al. 2009; O'Sullivan et al. 2009) or DSA in the backflowing material (Matthews et al. 2019). However, in our simulations, acceleration always occurs within the jet (i.e., in relativistic regions where τ = 1, i.e., pure jet material), consistently with the espresso scenario.

These results suggest that when an ultrarelativistic jet spine is present, such as in FR II and in subkiloparsec scales of FR I AGNs, acceleration proceeds via the espresso mechanism, typically up to the Hillas limit. Our simulations cannot exclude that stochastic acceleration might play some role when the jet becomes trans/nonrelativistic, although such processes have never been quantified via bottom-up calculations in environments with self-consistent MHD configurations, but only with Monte Carlo simulations where jet/cocoon structures are prescribed and pitch-angle scattering is imposed by hand.

Our simulations, while showing that espresso acceleration in relativistic jets can be more prominent than stochastic processes for UHECR production, of course cannot prove that stochastic acceleration cannot ever matter. For instance, it is possible that stochastic acceleration might reenergize the seed population and/or produce lower-energy particles responsible for the observed radio emission (e.g., Stawarz & Ostrowski 2002).

Figure 9 shows the cumulative distribution of the final energy gains of particles with ρ_i ∈ 3 × [10³, 10⁶] GV, for both case A and case B. If we consider the whole injection domain defined in Section 3.1 (upper curves), we find that ∼38% (∼53%) of case A (B) particles gain at least a factor of 2 in energy, and about ∼14% (∼18%) of case A (B) particles gain a factor of ${{\rm{\Gamma }}}_{\mathrm{eff}}$². Also, ∼0.06% (∼0.18%) of the particles achieve an energy gain of 100 or more in case A (B), corresponding to ${ \mathcal E }\gtrsim {{\rm{\Gamma }}}_{\mathrm{eff}}^{4}$, 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.

Zoom In Zoom Out Reset image size

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 τ < 0.1) 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 ∼14% (∼27%) of case A (B) particles gain at least a factor of 2 in energy, and about ∼2% (∼4%) of case A (B) particles gain a factor of ${{\rm{\Gamma }}}_{\mathrm{eff}}$², as shown in Figure 9. We also find that ∼0.06% 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, diffusion is expected to be effective at some level in realistic systems. 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 with α_i ≳ 0.01. In summary, we conclude that a fraction η ≈ 3%–20% 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 η≳10⁻⁴ 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.

FR I jets tend to wobble and exhibit multiple kinks (e.g., Mignone et al. 2010; Tchekhovskoy & Bromberg 2016), which may favor multiple acceleration cycles and allow the acceleration of UHECRs in sources such as Centaurus A (see the discussion in Section 3.5).

The espresso model predicts that the chemical composition observed in galactic CRs, which is increasingly heavy above 10¹³ eV and dominated by iron nuclei around 10¹⁷ eV (e.g., Hörandel et al. 2006; Kampert & Unger 2012), should be mapped into UHECRs. Auger observations suggest a proton-dominated flux at 10¹⁸ 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 ${ \mathcal E }\gtrsim {{\rm{\Gamma }}}_{\mathrm{eff}}^{2}$ to eliminate reaccelerated particles that should not contribute to the UHECR flux.

Zoom In Zoom Out Reset image size

The cutoffs of seed species are boosted up in energy by a factor of ${{\rm{\Gamma }}}_{\mathrm{eff}}$² ∼ 10, consistently with the results in Section 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 10¹⁸ eV (e.g., Gaisser et al. 2013; Aloisio et al. 2014; Taylor 2014) favors rather hard injection spectra (∝E^−1.5 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 α_i ≪ 1 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 10¹⁸ 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 ${ \mathcal E }\sim {10}^{3}$ energy gain required to produce 10²⁰ 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.

Although we initialized our fiducial jet with Γ₀ = 7, the Lorentz factor that we would infer from a multiwavelength analysis of its emission would be smaller (see Figure 6). Radio observations do not rule out the existence of an extended, ultrarelativistic jet spine even in ordinary AGNs (e.g., Chiaberge et al. 2001; Ghisellini et al. 2005); moreover, they suggest that FR I and FR II spines should have no appreciable differences in their initial Lorentz factors (Giovannini et al. 2001; Casandjian & Grenier 2008). In general, high-resolution relativistic MHD simulations (e.g., Mignone et al. 2010) show that the relativistic spine may extend for hundreds of jet radii up to the termination shock at the jet's head; even if deceleration is observed in the outer layers of jets, the Lorentz factor of the spine typically remains unchanged from the injection value (Rossi et al. 2008). This picture should mainly apply to FR II jets since the relativistic component exceeds many kiloparsecs in length as we discuss in Section 1.4. Modeling hyper-relativistic flows in 3D MHD simulations is extremely challenging, so we resort to extrapolating 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 ${ \mathcal E }\propto {{\rm{\Gamma }}}^{2}$ up to ${{\rm{\Gamma }}}_{\mathrm{eff}}\sim 3.2$ and then levels off because the filling factor of regions with ${\rm{\Gamma }}\gtrsim {{\rm{\Gamma }}}_{\mathrm{eff}}$ is small. By the same token, it is reasonable that also a multiwavelength analysis would preferentially highlight regions with Γ ≲ ${{\rm{\Gamma }}}_{\mathrm{eff}}$; 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 ${{\rm{\Gamma }}}_{\mathrm{eff}}$ ≈ Γ_AGN, about 10% of the CR seeds should consistently achieve energy gains as large as ${{\rm{\Gamma }}}_{\mathrm{AGN}}^{2}$ (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 ρ_knee ≈ 3 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 ${\rho }_{\mathrm{knee}}\propto {n}_{g}^{1/6}{M}_{g}^{1/2}$, where n_g and M_g are the number density and total mass of the galaxy considered. 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 in a more effective confinement of both CRs accelerated in SNRs and low-energy reaccelerated ones. Such an effect may change the slope and normalization of the seed spectrum, but it would systematically go in the direction of making it flatter and enhanced.

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 might preferentially come from AGN jets that point at us (from BL Lac and flat-spectrum radio quasars), while in the second case all radio-loud AGNs 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, determine whether reaccelerated particles are preferentially released quasi-isotropically or beamed along the jet axis.

Figure 11 shows how reaccelerated particles with ${ \mathcal E }\gt {{\rm{\Gamma }}}_{\mathrm{eff}}^{2}$ are released. The vectors in the top panels illustrate the final directions of flight, and the bottom panels show the distribution of ${\mu }_{{\rm{f}}}$, the cosine of the angle between the final particle velocity and the z-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 case 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, E_r = v_zB_ϕ, which reverses its sign in the two cases. In fact, particles preferentially move along the +z-direction while in the spine, but their final direction is determined by the deflection they experience after leaving the relativistic regions and entering the cocoon, where the flow is in the −z-direction (top panel of Figure 5). A different sign of B_ϕ in the jet produces a radial electric field that either disperses (E_r > 0, case A) or collimates (E_r < 0, case B) the flux of particles that escape the jet (Figure 4). Therefore, the final anisotropy of the particles is not determined by pitch-angle scattering in the cocoon turbulence, but rather by the global sign of the motional electric field. Note that also in case B about half of the highest-energy particles escape within an angle of ∼30°, a bit larger than the canonical relativistic beaming of ∼${{\rm{\Gamma }}}_{\mathrm{eff}}$⁻¹ ∼ 20°. Also considering that the sign of B_ϕ 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.

Zoom In Zoom Out Reset image size

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 10¹⁸ eV, with the galactic to extragalactic transition populated with secondary protons, mainly. Fang & Murase (2018) recently suggested that photodisintegration, photomeson production by secondary nucleons, and Bethe–Heitler pair production processes may also account for the flux of high-energy neutrinos measured by IceCube and for the GeV gamma-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 extragalactic CRs starting with relatively steep seed spectra (e.g., ∝E⁻²), with a low-rigidity cutoff for heavy nuclei and the contribution of secondary protons around 10¹⁸ eV. Note that this still does not explain the reason why the galactic and extragalactic 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 toward 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 toward the jet base, so we expect photodisintegration to act as a high-pass filter that may make the UHECR injection spectrum flatter at low energies in the central regions of radio-loud quasars.

The results of Section 4.4 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 kiloparsecs) is in fact much larger than the particle Larmor radius even at the highest rigidities (${ \mathcal R }\sim 1$ kpc at 10¹⁸ V in a μG field). Note that this does not mean that the Hillas limit is larger in the lobes, since there flows are nonrelativistic 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., nonblazar AGNs. For instance, Centaurus A is usually classified as an FR I AGN with Γ ≲ 7, despite likely hosting a very relativistic jet with Γ ≈ 15–20 (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 AGNs 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, 2018). The only marginal (3.4σ post-trial) evidence for a hotspot has been reported by the Telescope Array Collaboration for events above 5.7 × 10¹⁹ 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 L_bol ≈ 2 × 10⁴⁴ erg s⁻¹ and a luminosity distance of about 134 Mpc (see Caprioli 2018, for a more extended discussion).

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, and hence nonblazar AGNs may also contribute; (ii) when just a few shots are allowed, even AGN jets with moderate Γ ≳ 3 produce UHECRs (note that the Hillas criterion does not depend on Γ); and (iii) self-confinement and propagation delay may offset the UHECR arrival with respect to the signature of a prominent jet activity.