Ultrahigh-energy Cosmic-ray Nuclei from Black Hole Jets: Recycling Galactic Cosmic Rays through Shear Acceleration

TL;DR

This paper models how shear acceleration in active galactic nuclei jets can reaccelerate galactic cosmic rays to ultrahigh energies, producing a hard spectrum and explaining observed composition and anisotropy.

Contribution

It introduces a Monte Carlo simulation of shear acceleration in AGN jets, demonstrating its effectiveness in producing ultrahigh-energy cosmic rays with specific spectral and compositional features.

Findings

Reacceleration yields a hard cosmic ray spectrum ($E^{-1}$ to $E^0$).

High-energy nuclei abundance is enhanced due to TeV-PeV injections.

Model aligns with observed cosmic ray spectrum, composition, and anisotropy data.

Abstract

We perform Monte Carlo simulations of transrelativistic shear acceleration dedicated to a jet-cocoon system of active galactic nuclei. A certain fraction of galactic cosmic rays in a halo is entrained, and sufficiently high-energy particles can be injected to the reacceleration process and further accelerated up to 100 EeV. We show that the shear reacceleration mechanism leads to a hard spectrum of escaping cosmic rays, $dL_E/dE\propto E^{-1}-E^0$, distinct from a conventional $E^{-2}$ spectrum. The supersolar abundance of ultrahigh-energy nuclei is achieved due to injections at TeV-PeV energies. As a result, we find that the highest-energy spectrum and mass composition can be reasonably explained by our model without contradictions with the anisotropy data.

Tables

Table 1. Table 1: The parameter sets for the models shown in Figs. 2 , 5 , and 6

models	$R_{\mathrm{jet}}$111In unit of kpc.	${\xi }_{c/j}$	$l_{\mathrm{coh}}/R_{\mathrm{coc}}$	$Z_i$	$B_{\mathrm{coc}}$222in unit of $\mu$G.	${\beta }_{\mathrm{jet}}$	$E_{\mathrm{coh},i}$333in unit of EeV.
Reference	0.5	10	0.03	1	3	0.7	0.42
A-1	0.5	10	0.03	26	3	0.7	11
A-2	0.5	10	0.03	1	15	0.7	2.1
A-3	0.5	10	0.03	1	3	0.5	0.42
B-1	0.5	4	0.03	1	3	0.7	0.17
B-2	1.5	10	0.03	1	3	0.7	1.2
C-1	0.5	10	0.003	1	3	0.7	0.042
C-2	0.5	100	0.003	1	3	0.7	0.42

Equations

$$t_{\mathrm{acc}}=\frac{\Delta t}{(\Delta E/E)}\sim {\zeta }_a\frac{{\lambda }_{i,\mathrm{coc}}}{c{\Gamma }_j^2{\beta }_{\mathrm{jet}}^2},$$

$$t_{\mathrm{conf}}={\zeta }_c\frac{R_{\mathrm{jet}}}{c},$$

$$E_{i,\max }\approx \zeta e{Z}_i{B}_{\mathrm{coc}}{l}_{\mathrm{coc}}^{1/2}{R}_{\mathrm{jet}}^{1/2}{\Gamma }_{\mathrm{jet}}{\beta }_{\mathrm{jet}},$$

$$\frac{\partial f}{\partial t}=\frac{1}{p^2}\frac{\partial }{\partial p}\left(p^2{D}_p\frac{\partial f}{\partial p}\right)-\frac{f}{t_{\mathrm{esc},\mathrm{sl}}}+Q_0\delta (p-p_{\mathrm{inj},\mathrm{sl}}),$$

$$s_{\rm sl}=\left\{\begin{array}[]{ll}\frac{\delta-1}{2}-q_{\rm sl}&(p<p_{\rm inj,sl})\\ \frac{\delta-1}{2}+q_{\rm sl}&(p>p_{\rm inj,sl})\end{array}\right.,$$

$$q_{\mathrm{sl}}=\sqrt{\frac{(\delta +3{)}^2}{4}+\frac{t_{\mathrm{acc},\mathrm{sl}}}{t_{\mathrm{esc},\mathrm{sl}}}}.$$

$$n_{i,d}=K_i{\left(\frac{E_{i,\mathrm{inj}}}{\mathrm{TeV}}\right)}^{-{\alpha }_i+1}\exp \left(-\frac{E_{i,\mathrm{inj}}}{Z_i\mathrm{PeV}}\right).$$

$${\dot{N}}_{i,\mathrm{inj}}\approx \frac{N_{i,\mathrm{inj}}}{t_{\mathrm{ad}}}\approx \frac{2\pi R_{\mathrm{jet}}^2{H}_h{n}_{i,d}}{t_{\mathrm{ad}}}.$$

$$E_{i,\max }\approx \zeta e{Z}_i{B}_{\mathrm{coc}}{l}_{\mathrm{coh}}^{1-1/\delta }{R}_{\mathrm{jet}}^{1/\delta }{\Gamma }_{\mathrm{jet}}^{2/\delta }{\beta }_{\mathrm{jet}}^{2/\delta }.$$

Full text

I Introduction

The origin of ultrahigh-energy cosmic rays (UHECRs) has been under intense debate for more than half a century Linsley (1963). Observationally, remarkable developments was made by High-Resolution Fly’s Eye, Pierre Auger Observatory (PAO), and Telescope Array (TA) Kotera and Olinto (2011). The spectrum of UHECRs has a cutoff around 60 EeV Abbasi et al. (2008); Abraham et al. (2008), consistent with the energy of the Gresien-Zatspin-Kuzmin cutoff for protons Greisen (1966); Zatsepin and Kuzmin (1966) or a photodisintegration cutoff for irons. The absence of small-scale anisotropy in the arrival direction of the highest-energy cosmic rays (CRs) places a lower limit on the number density of UHECR sources, $n_{s}\gtrsim{10}^{-6}-{10}^{-5}{\rm~{}Mpc}^{-3}$ Kashti and Waxman (2008); Takami and Sato (2009); Takami et al. (2012); Aab et al. (2015), which may be stronger at 10 EeV Takami et al. (2016). The CR composition is estimated from the depth of the shower maximum, $X_{\rm{max}}$ Abraham et al. (2010); Aab et al. (2014); Aab (2015); Fukushima (2015). The results of PAO and TA on the mean depth $\langle{X_{\rm{max}}}\rangle$ seem compatible Abbasi et al. (2016). However, the interpretation of the data is controversial, partly due to uncertainty in hadron physics implemented in extensive air-shower simulations. The latest hadron interaction models imply that the composition gradually becomes heavier for $\gtrsim 3$ EeV. On the other hand, UHECRs are believed to be dominated by light elements around $1-3$ EeV Abbasi et al. (2010, 2015).

The heavy mass composition, if true, challenges astrophysical models for UHECR sources. The simultaneous fittings of the spectrum and composition suggest (i) a mass composition heavier than the solar abundance and (ii) a spectrum harder than a conventional $E^{-2}$ spectrum for the plausible redshift evolution Allard (2012); Taylor et al. (2015); Aab et al. (2016). The former difficulty is more serious if any anisotropy is established for the highest-energy CR nuclei, since the similar level of the anisotropy is expected at the same rigidity of protons ($\sim 1-10$ EeV) Lemoine and Waxman (2009); Abreu et al. (2011); Liu et al. (2013) (see also Takami and Sato (2009)). Gamma-ray bursts Murase et al. (2008); Horiuchi et al. (2012); Wang et al. (2008); Globus et al. (2015) and newborn pulsars Fang et al. (2012, 2014) can provide a metal-rich composition, but UHECR nuclei must survive against photodisintegration in the sources and the intrinsic abundance ratio is essentially treated as a free parameter. Also, the luminosity argument Blandford (2000) and the nondetection of ultrahigh-energy neutrinos Aartsen et al. (2016) disfavor steady UHECR proton sources Murase and Takami (2009); Takami et al. (2016); Fang and Kotera (2016). Steady sources, such as active galactic nuclei (AGNs), can accelerate CR nuclei up to 100 EeV (e.g., Murase et al. (2012, 2014)), but the origin of heavy composition with a hard spectrum has been unclear.

In this work, we provide a new scenario that simultaneously explains the spectrum and composition, overcoming the above difficulties. First, in Section II, we consider the shear acceleration of CRs around transrelativistic shear layers, where both discrete and continuous shear acceleration mechanisms are discussed. For high-energy CRs, we perform detailed numerical simulations and show that each composition species of the CRs leaving the accelerators have a hard spectrum with a rigidity dependent energy cutoff at $E_{i,\rm{max}}=Z_{i}E_{p,\rm{max}}$ (i.e. Peters cycle Peters (1961); Gaisser et al. (2013)). Here $E_{i,\rm{max}}$ is the ion maximum energy at the CR accelerators, $E_{p,\rm{max}}$ is the proton maximum energy, and $Z_{i}$ is the particle charge for a particle species with $i$. For low-energy CRs, we discuss the analytical CR spectrum, and show that the high-energy CRs can be accelerated mainly via the discrete shear acceleration mechanism. Then, in Section III, we apply the mechanism to the system that is composed of an AGN jet and a cocoon inflated by the jet. We find that TeV-PeV CRs injected from a galactic halo are naturally accelerated by the shear acceleration, which can generate UHECRs with energies up to 100 EeV (Fig. 1). We also calculate the UHECR propagation in intergalactic space, and demonstrate that our model accounts for the observed Auger data well. In Section IV, we summarize our results and discuss implications.

II Shear Acceleration

Shear acceleration is a class of Fermi acceleration mechanisms Berezhko and Krymskii (1981); Earl et al. (1988); Webb (1989); Ostrowski (1990); Rieger and Duffy (2006). The shear acceleration occurs when the relativistic particles are inside an ordered shear velocity field, which is commonly expected in the astrophysical jets Ostrowski (1998); Rieger and Duffy (2004, 2016) and accretion flows Katz (1991); Subramanian et al. (1999); Kimura et al. (2016). When the shear is continuous in the scale of the mean free path for scatterings with magnetic fields, the acceleration mechanism is basically the same as the stochastic acceleration in a turbulence. A particle that has a head-on (tail-on) collision gains (loses) energy, and the particles are statistically accelerated because the head-on collision is more probable than the tail-on collision Fermi (1949); Subramanian et al. (1999). When the scattering mean free path is longer than the scale of the shear velocity gradient, the acceleration is regarded as the Fermi process in the discrete shear Ostrowski (1990, 1998). In our scenario, UHECR production proceeds in this regime due to their large Larmor radii. The spatial diffusion is important, so that we take a numerical approach to properly consider the geometry. Note that the continuous shear acceleration and discrete shear acceleration are different in terms of the properties of CR acceleration, which leads to the important difference in their time scales such as the CR escape time and CR acceleration time. This may result in distinct predictions for CR spectra.

II.1 Discrete shear acceleration

II.1.1 Setup for Monte Carlo simulations

We consider a jet-cocoon system (see, e.g., Begelman and Cioffi (1989); Kawakatu and Kino (2006)). To mimic the geometry of interest (see Fig. 1), we consider two cylinders with radii of $R_{\rm{jet}}$ and $R_{\rm{coc}}$. We parameterize the cocoon radius as $R_{\rm{coc}}\equiv\xi_{c/j}R_{\rm{jet}}$. The shear between the jet and cocoon is given by the jet velocity, $c\beta_{\rm{jet}}$. The cocoon is quasi-spherical in general. For simplicity, we assume the vertical length of the jet and the cocoon to be equal to the cocoon radius: $l_{\rm{jet}}=l_{\rm{coc}}=R_{\rm{coc}}$, which is sufficient for the purpose of this work.

We expect that both of the jet and cocoon have turbulent magnetic fields that scatter the particles. We can parameterize the mean free path inside the cocoon as $\lambda_{i,\rm{coc}}=(E/E_{i,\rm{coh}})^{\delta}l_{\rm{coh}}$, where $l_{\rm{coh}}$ is the coherence length and $E_{i,\rm{coh}}=Z_{i}eB_{\rm{coc}}l_{\rm{coh}}$ ($B_{\rm{coc}}$ is the magnetic field strength in the cocoon). The particles are resonantly scattered by turbulence for $E<E_{i,\rm{coh}}$, which leads to $\delta=1/3$ if we assume the Kolmogorov turbulence inside the cocoon (e.g. Stawarz and Petrosian, 2008). On the other hand, particles are scattered in non-resonant manner with small-scale turbulence for $E>E_{i,\rm{coh}}$, resulting in $\delta=2$ (e.g. Sironi et al., 2013). Both the turbulence and magnetic field are likely to be strong in the jet, and the diffusion process in the strong turbulence is likely to be the Bohm limit Roh et al. (2016); Kimura et al. (2016). Thus, we use the Bohm limit there, $\lambda_{i,\rm{jet}}=E/(Z_{i}eB_{\rm{jet}})$, where $B_{\rm{jet}}$ is the magnetic field strength in the jet. The particles move in a manner of the random walk by these interactions, and undergo multiple passage through the shear layer. This results in the discrete shear acceleration.

For a given nuclear species, we inject 262,144 particles with an injection energy of $E_{i,\rm{inj}}$ (see Section III.1) at the jet-cocoon boundary at $t=0$, and track them by a time of $t=t_{\rm{ad}}\approx{R_{\rm{coc}}}/v_{\rm{exp}}$, where $v_{\rm{exp}}$ is the expansion velocity of the cocoon. After this time scale, we expect that the particles lose their energies due to the adiabatic expansion. Since injected particles are reaccelerated to ultrahigh energies, more than 89 % of the particles escape from the system by the end of simulation runs. The number of the injected particles is normalized by the injection rate $\dot{N}_{\rm{inj}}$ (see Section III.1). The particles travel straightly until they are scattered by a magnetic field. The scattering angle distribution is assumed to be isotropic in the rest frame of each fluid, which is a simplified but reasonable approximation in our problem, given that almost all the particles experience many scatterings during their residence time (cf. Kato and Takahara (2001, 2003); Aoi et al. (2008) and references therein). When the particles diffuse out beyond the cocoon radius, $R_{\rm{coc}}$, or the jet length, $l_{\rm{jet}}$, they are recorded as “escaping” particles.

Hereafter, we consider radio-loud AGNs, in particular Fanaroff-Riley I radio galaxies (FR Is), to demonstrate our results (see Appendix B for an application to radio-quiet AGNs). Powerful kiloparsec-scale jets are commonly seen in radio galaxies, and they are often accompanied by radio lobes or bubbles. The jets sweep up the circumgalactic materials in galactic halos, and eventually propagate into the intergalactic medium. The plasma inflated by the jet forms a cocoon, which is attributed to a radio lobe or bubble. The length of the jet depends on the age of AGN, and we consider the time when the jet finishes sweeping the halo in which galactic CRs (GCRs) are confined, i.e., $l_{\rm{jet}}=H_{h}$, where the scale height of the CR halo $H_{h}$ is set to 5 kpc Strong et al. (2007). Kiloparsec-scale jets of FR Is are only mildly relativistic, so that the jet velocity is set to $\beta_{\rm{jet}}=0.7$ (e.g., Bowman et al. (1996); Laing and Bridle (2002); Canvin et al. (2005)). The ratio of the cocoon to jet radii is given by $\xi_{c/j}=10$ as a reference value Mizuta et al. (2004); Rossi et al. (2008), which leads to $R_{\rm{jet}}=0.5$ kpc. The magnetic fields are assumed to be $B_{\rm{jet}}=0.3\rm~{}mG$ (e.g., Stawarz et al. (2006)) and $B_{\rm{coc}}=3~{}\mu$G Dunn et al. (2005); Kataoka and Stawarz (2005). The expansion velocity of the cocoon is set to $v_{\rm{exp}}=3000\rm~{}km~{}s^{-1}$ Bordas et al. (2011). For the coherence length, we use $l_{\rm{coh}}=0.03R_{\rm{coc}}$ as a reference value. These fiducial parameters are consistent with the observations of radio galaxies, and the spectral shape is largely unaffected by the change of the parameters.

II.1.2 Maximum energy in discrete shear acceleration

In our cases, the jet confines the CRs more efficiently than the cocoon, i.e. $R_{\rm{jet}}/\lambda_{i,\rm{jet}}>R_{\rm{coc}}/\lambda_{i,\rm{coc}}$, which means that the diffusion in the cocoon determines the maximum energy. If the particle diffusing inside the cocoon returns to the jet, the acceleration cycle continues. Otherwise, the particle escapes from the system. Fig. 2 shows the probability of a diffusing particle returning to the jet as a function of $\lambda_{i,\rm coc}/R_{\rm jet}$ for various parameter sets tabulated in Table 1. The lines completely overlap each other except for B-1 and C-2 that have different values of $\xi_{c/j}$. This indicates that the return probability depends only on $\lambda_{i,\rm coc}/R_{\rm jet}$ and $\xi_{c/j}$ under the assumption of $l_{\rm jet}=R_{\rm coc}$. Note that the return probability depends on the scattering mean free path not in the jet but in the cocoon.

When $\lambda_{i,\rm{coc}}\lesssim{R}_{\rm{jet}}$, the majority of the diffusing particles return to the jet after a few random steps. This feature does not change regardless of the physical parameters of the jet-cocoon system, as long as $\lambda_{i,\rm coc}<R_{\rm jet}$. Then, the acceleration time is expressed as

[TABLE]

where $\Delta{t}=\zeta_{a}\lambda_{i,\rm{coc}}/c$ is the typical residence time in the cocoon per cycle and $\Delta{E}/E\sim\Gamma_{\rm{jet}}^{2}\beta_{\rm{jet}}^{2}$ is the mean energy gain per cycle Rieger and Duffy (2004). Here $\zeta_{a}$ is a correction factor that accounts for the average number of steps over the accelerated particles. The residence time in the jet per cycle is much shorter than that in the cocoon, because of the shorter mean free path in the jet.

On the other hand, when $\lambda_{i,\rm{coc}}\gtrsim{R}_{\rm{jet}}$, the majority of the CR particles escape from the cocoon without returning to the jet, as seen in Fig. 2. Only the particles that go back to the jet continue to gain energies by the shear. Thus, the size of efficient CR acceleration region is limited by the jet size. The effective confinement time in the acceleration region can be represented as

[TABLE]

where $\zeta_{c}=\zeta_{c}(\xi_{c/j})$ is a geometrical correction factor that takes into account the weak dependence on $\xi_{c/j}$. As seen in Fig. 2, CRs have a larger chance to return to the jet for a larger $\xi_{c/j}$. Since the CRs escape through the cocoon, the confinement time itself is not directly related to the mean free path in the jet.

The condition $t_{\rm{acc}}\approx t_{\rm{conf}}$ leads to the maximum energy in the energy spectrum of escaping CRs:

[TABLE]

where $\zeta\equiv(\zeta_{c}/\zeta_{a})^{1/2}$ and $\lambda_{i,\rm{coc}}\propto{E^{2}}$ is used. From the simulation results, we found $\zeta\simeq 2.2(\xi_{c/j}/10)^{0.2}$ (see Appendix A for the consistency of this estimate and the simulation results), leading to $E_{i,\rm{max}}\sim 1.6Z_{i}\rm~{}EeV$ for our reference parameter set (see Fig. 3). We confirm this scaling relation for mildly relativistic cases of $\Gamma_{\rm{jet}}\beta_{\rm{jet}}\sim 1$ 444Note that for highly relativistic cases, the anisotropy in the momentum distribution would considerably affect $E_{i,\rm{max}}$ in our estimate..

The discrete shear acceleration process is one of the Fermi acceleration mechanisms, so the accelerated CRs have a power-law spectrum. Almost all the accelerated particles can escape. For $E<E_{i,\rm{max}}$, the escaping CRs show a hard power-law spectrum, $dL_{E}/dE\propto{E^{-1}-E^{0}}$ (see Fig. 3). It has a spectral break at $E\sim{E_{i,\rm{coh}}}$ due to the change of energy dependence of the mean free path. For $E>E_{i,\rm{max}}$, the spectrum has a cutoff that is slower than the exponential (see Appendix A for the detailed results of Monte Carlo simulations, including the parameter dependence of the spectral shape and cases for the Bohm limit). Since we consider kiloparsec-scale jets, we can neglect energy losses due to proton synchrotron, hadronuclear, photohadronic, and photodisintegration processes.

II.2 Continuous shear acceleration

There is a shear layer between the jet and the cocoon where the jet velocity may change linearly Aloy and Mimica (2008). This layer affects the spectrum of CRs if the size of shear layer is larger than the Larmor radius or the scattering mean free path of the CRs Rieger and Duffy (2006). Here, we make a brief discussion about effects of the shear layer, which may have a crucial influence on the injection process to the discrete shear acceleration (see Section III.1).

Inside the shear layer, the evolution of distribution function is described by the diffusion equation in momentum space. Adding the escape term and injection term, which are important in our setup, we can write the CR transport equation as Earl et al. (1988); Rieger and Duffy (2006)

[TABLE]

where $D_{p}\approx p^{2}\lambda_{i,\rm sl}c(dv_{j}/dr)^{2}/15$ is the diffusion coefficient in momentum space ($\lambda_{i,\rm sl}$ is the mean free path and $dv_{j}/dr$ is the velocity gradient in the shear layer), $t_{\rm esc,sl}=R_{\rm sl}^{2}/(2\lambda_{\rm sl}c)$ is the escape time from the shear layer ($R_{\rm{sl}}$ is the size of shear layer), $Q_{0}$ is the injection rate, and $p_{\rm inj,sl}$ is the injection momentum. The acceleration time is estimated to be $t_{\rm acc,sl}=p^{2}/D_{p}\propto p^{-\delta}$, where we write the mean free path as $\lambda_{i,\rm sl}\approx\lambda_{0}(p/p_{0})^{\delta}$. This dependence is the same as that of $t_{\rm esc}$, which means that the acceleration time is shorter for higher energy for $\delta>0$ Rieger and Duffy (2006). Assuming a power-law distribution function $dN/dE=4\pi p^{2}f\propto p^{-s_{\rm sl}}$, we can obtain the steady state solution as

[TABLE]

[TABLE]

We confirm this power-law solution by numerically solving Eq. (4). The spectrum of escaping particles is written as $dL_{E}/dE\approx(dN/dE)/t_{\rm esc,sl}\propto E^{-s_{\rm esc,sl}}$, where $s_{\rm esc,sl}=\mp q_{\rm sl}-(1+\delta)/2$. Considering the linear velocity gradient, $(dv_{j}/dr)\approx c\beta_{\rm jet}/R_{\rm sl}$, we obtain $t_{\rm acc,sl}/t_{\rm esc,sl}=30c^{2}/(R_{\rm sl}^{2}(dv_{j}/dr)^{2})\approx 30\beta_{\rm jet}^{-2}$. Then, the index of the escape spectrum is $s_{\rm esc,sl}\sim 7.3~{}(-8.7)$ for $p>p_{\rm inj}$ ($p<p_{\rm inj}$). This spectrum is so steep that it cannot match the observed UHECR spectrum. Most of the injected particles escape from the shear layer before being accelerated to higher energies. In other words, only few low-energy GCRs that are injected to the continuous shear acceleration can reach the injection energy, above which the discrete shear acceleration operates (see Section III.1). Thus, the low-energy GCRs are unlikely to be accelerated to UHECRs. Here, we assume that the particles are injected at the center of the shear layer for simplicity. In reality, the particles are injected at the edge of the shear boundary. Although this could affect the spectral shape, it is unlikely that the injection position drastically changes the acceleration efficiency. More detailed discussions for the continuous shear acceleration are beyond the scope of this work, and remains as a future work.

III Recycling Galactic CRs as UHECRs

III.1 Injection rate and composition ratio

In our shear reacceleration scenario, we have shown that the spectrum of escaping CRs is generically hard, and $E_{i,\rm{max}}$ is determined by the five parameters ($\beta_{\rm{jet}}$, $R_{\rm{jet}}$, $\xi_{c/j}$, $l_{\rm{coh}}$, $B_{\rm{coc}}$). Next, we estimate the UHECR luminosity and their composition ratio.

CR densities in radio galaxies are highly uncertain. Here, we assume that the proton CR densities are comparable to that in our Galaxy. While the star-formation rate of elliptical galaxies may be lower than that of star-forming galaxies by a factor of 3–10 Salmi et al. (2012); Martig et al. (2013), this uncertainty is easily absorbed by uncertainties in the other parameters. The GCR density inside the CR halo of $H_{h}\sim 5$ kpc Strong et al. (2007) can be expressed as

[TABLE]

Here, CR species are grouped as $i=$ H, He, C–O, Ne–Al, Si–K, Ca–Mn, Fe. Their effective charge $Z_{i}$ and atomic mass $A_{i}$ are $Z_{i}=$ 1, 2, 7, 11, 15, 23, 26 and $A_{i}=$ 1, 4, 14, 23, 30, 49, 56, respectively. We use the observed values at $E\sim 1$ TeV for the normalization of each component: $K_{\rm{H}}=3.6\times 10^{-15}\rm~{}cm^{-3}$ and $K_{i}/K_{\rm{H}}\simeq$ 1, 0.65, 0.33, 0.17, 0.14 0.072, 0.23 Wiebel-Sooth et al. (1998); Hoerandel (2003). In the galactic disk, the proton has softer index than the others Wiebel-Sooth et al. (1998); Hoerandel (2003); Caprioli et al. (2011); Yoon et al. (2011), $\alpha_{\rm{H}}\simeq 2.7$ and $\alpha_{i\neq\rm{H}}\simeq 2.6$ 555This is conservative and would be more appropriate for older galaxies. Harder spectra assumed in Caprioli (2015) relaxes the energetics requirement. In addition, we increase the abundance of nuclei heavier than He by factor of 3 from the value above because most of radio galaxies have more metals than the Galaxy due to their past star formation activities Henry and Worthey (1999); Tang et al. (2009).

The number of swept-up particles of species $i$ by the time when $l_{\rm{jet}}=H_{h}$ is simply given by $2\pi R_{\rm{coc}}^{2}H_{h}n_{i,d}$, and we assume that only the fraction, $R_{\rm{jet}}^{2}/R_{\rm{coc}}^{2}$, is injected into shear acceleration. Thus, the time-integrated number of injected GCRs are written as $N_{i,\rm inj}\approx 2\pi R_{\rm jet}^{2}H_{h}n_{i,d}$. The swept-up particles of $\lambda_{i,\rm sl}<R_{\rm{sl}}$ are accelerated by the continuous shear that is ineffective to produce high-energy CRs (see Section II.2). Only the particles of $\lambda_{i,\rm sl}>R_{\rm{sl}}$ can be injected to the discrete shear acceleration process. Setting $\lambda_{i,\rm sl}=R_{\rm sl}$, the injection energy is given by $E_{i,\rm{inj}}\approx E_{\rm coh}(R_{\rm{sl}}/l_{\rm{coh}})^{3}\sim 15Z_{i}\rm~{}TeV$. Here, we use $\lambda_{i,\rm sl}\sim\lambda_{i,\rm coc}$ and $R_{\rm{sl}}\sim 0.01R_{\rm{jet}}\sim 5$ pc. The injected CRs are accelerated until the adiabatic cooling is effective, $t_{\rm{ad}}\approx{R}_{\rm{coc}}/v_{\rm{exp}}\sim 1.6$ Myr (where $v_{\rm{exp}}\sim 3000\rm~{}km~{}s^{-1}$ Bordas et al. (2011)). The time-averaged injection rate of GCRs of species $i$ to shear acceleration is estimated to be

[TABLE]

Renormalizing the simulation input by the injection rate, we obtain the differential luminosity of UHECRs, $L_{\rm{UHECR}}$. The CR luminosity density at ${10}^{19.5}$ eV is $0.6\times{10}^{44}{\rm~{}erg~{}Mpc^{-3}~{}yr}^{-1}$ (e.g., Murase and Takami (2009)), and the number density of FR Is is roughly $\sim{10}^{-5}-{10}^{-4}\rm~{}Mpc^{-3}$ Padovani et al. (2011); Prescott et al. (2016). Thus, $L_{\rm{UHECR}}\sim 2\times 10^{40}-2\times{10}^{41}~{}\rm~{}erg~{}s^{-1}$ is required. Our model can satisfy this requirement, as shown in Fig. 3. Also, our model can avoid anisotropy constraints at $E\sim 10$ EeV Abu-Zayyad et al. (2012) owing to the high source number density with the heavy composition. The relative abundance ratio at the same rigidity is estimated to be $(f_{\rm{H}},~{}f_{\rm{He}},~{}f_{\rm{C\mathchar 45\relax O}},~{}f_{\rm{Ne\mathchar 45\relax Al}},~{}f_{\rm{Si\mathchar 45\relax K}},~{}f_{\rm{Ca\mathchar 45\relax Mn}},~{}f_{\rm{Fe}})=$ (0.73, 0.21, 0.042, 0.011, 0.0053, 0.0014, 0.0037). Note that we cannot freely change the abundance ratio among heavy nuclei as well as the intrinsic spectral index, because they are determined by the shear acceleration mechanism and observed abundance of Galactic CRs.

III.2 Comparison with observations

We calculate the propagation of the UHECRs from the sources to the Earth using CRPropa 3 Armengaud et al. (2007); Alves Batista et al. (2016). The code includes the photomeson production, the photodisintegration, and the electron-positron pair production through the cosmic microwave background and extragalactic background light (EBL). The nuclear decay process is also included. We use the EBL model of Finke et al. (2010), and assume that all FR Is produce the UHECRs shown in Fig. 3 for simplicity. The luminosity density of bright AGNs positively evolves with redshift Ajello et al. (2014); Ueda et al. (2014), while that of low-luminosity AGNs may have a weaker redshift evolution Ajello et al. (2014); Padovani et al. (2011). In this work, we assume no redshift evolution but stronger evolution models can also fit the data.

We show the spectrum of the UHECRs at the Earth in the upper panel of Fig. 4. The intermediate and heavy nuclei decrease while protons increase during the propagation process due to the photodisintegration. The cutoff at $E\gtrsim 100$ EeV is produced due to the maximum energy of the shear acceleration at the source, which is consistent with the PAO data. We need an additional component to fit the spectrum at $E\sim$ EeV (e.g., Katz et al. (2009); Aloisio et al. (2014); Thoudam et al. (2016)). The middle panel and the lower panel show the mean depth of the shower maximum, $\langle{X}_{\rm{max}}\rangle$, and variance of the shower depth, $\sigma(X_{\rm{max}})$, respectively. These values are calculated using $X_{\rm{max}}$ probability distribution parametrized by De Domenico et al. (2013). Within systematic errors, our model reasonably explains the observed feature of the chemical composition that changes from light to heavy as CR energy increases, without tuning the abundance ratio by hand.

IV Summary and discussion

We have shown that the shear acceleration by black-hole jets provides a promising mechanism of UHECR production. Based on the setup for the jet-cocoon system that is ubiquitous in radio galaxies, we have performed detailed numerical simulations of UHECR acceleration, escape, and propagation in intergalactic space. The radio galaxies can accelerate protons up to a few EeV and irons up to 100 EeV, whose spectra are intrinsically hard as required by the PAO data. TeV–PeV CRs in a galactic halo are injected to the shear acceleration, leading to the enhanced metal abundance suggested by the $\langle{X_{\rm{max}}}\rangle$ and $\sigma(X_{\rm{max}})$ data.

We stress that the spectrum and composition are essentially determined by theoretical calculations and observations of Galactic CRs, respectively. Although the calculation of propagation is slightly affected by the redshift evolution of the sources and the EBL model Aab et al. (2016), this cannot change our conclusion. It is possible to alter the source spectral index by superposing contributions from radio galaxies that have different $E_{i,\rm{max}}$. While more luminous radio galaxies could accelerate UHECRs to higher energies, $L_{\rm{UHECR}}$ is independent of the jet luminosity if all the radio galaxies have the same size of the halos. Then, fainter radio galaxies such as FR Is may give the most important contribution to the observed UHECR flux. On the other hand, the source parameters, such as $l_{\rm{coh}}$, $\beta_{\rm{jet}}$, and $R_{\rm{sl}}$, are uncertain. Phenomenologically, all the uncertainties are absorbed by treating $E_{p,\rm{max}}$ and $L_{\rm{UHECR}}$ as free parameters. The source models with similar values of $E_{p,\rm{max}}$ give the similar shape of the spectra at the Earth, $\langle{X}_{\rm{max}}\rangle$, and $\sigma(X_{\rm{max}})$. According to observations and simulations of the jet propagation Laing and Bridle (2002); Aloy and Mimica (2008), $R_{\rm{sl}}/R_{\rm{jet}}\sim 0.1$ and $\beta_{\rm{jet}}\sim 0.9$ are also possible, where we would need smaller $B_{\rm{coc}}$ and larger $l_{\rm{coh}}$ to obtain the required $E_{p,\rm{max}}$ and $L_{\rm{UHECR}}$.

We have considered shear acceleration in large-scale jets, which is different from the scenario by Murase et al. (2012) for UHECR acceleration in blazar jets. Our model is also different from Caprioli (2015), which relies on the first encounter boost in the relativistic jet of $\Gamma\sim 30$ Gallant and Achterberg (1999), whereas both consider the injection of Galactic CRs. While such jets could exist in sub-parsec scales as suggested in blazars or even kiloparsec scales for the most powerful FR II galaxies, jets of FR Is are significantly decelerated in such large scales, and mildly relativistic jets are considered in this work Bowman et al. (1996); Laing and Bridle (2002); Canvin et al. (2005).

FR Is and their blazar counterparts, BL Lac objects, are observed at different wavelengths from radio to gamma-rays. The charged particles that emit the observed electromagnetic signals are likely to be produced at different locations in the shear layer, e.g., by internal shocks Rees (1978) or turbulence Kakuwa et al. (2015). In the leptonic scenario, the electrons are difficult to get accelerated solely by the discrete shear acceleration mechanism, since their typical energy is lower than $E_{i,\rm{inj}}$ Inoue and Takahara (1996).

Our model is consistent with the convergence picture of UHECRs, neutrinos, and gamma rays Murase and Waxman (2016); Fang and Murase (2017), in which all three messengers are explained simultaneously. In the galaxy cluster and group model, UHECRs can be provided by AGNs Fang and Murase (2017). CRs that do not reach ultrahigh energies can be accelerated by the AGN jet without the shear reacceleration, and the CR spectrum can be effectively extended to ultrahigh energies with a hard spectrum via the shear acceleration mechanism. Also, the corresponding cosmogenic neutrino flux is expected to be $\sim{10}^{-10}{\rm~{}GeV~{}cm^{-2}~{}s^{-1}~{}sr^{-1}}$. Gamma rays and neutrinos associated with large scale jets may not be easy to detect due to long energy-loss time scales (cf. Ostrowski (2000)). Whereas electrons may be difficult to be injected into the shear acceleration process, it is important to study indirect signatures through radio and/or X-ray observations Stawarz and Ostrowski (2002); Laing and Bridle (2002); Stawarz et al. (2006) for testing our model.

Acknowledgements.

We acknowledge Damiano Caprioli, Martin Lemoine, Miguel Mostafa, Michael Ostrowski, and Günter Sigl for useful discussion. The work is supported by Alfred P. Sloan Foundation, NSF Grant No. PHY-1620777 (K.M.), NASA NNX13AH50G, and the IGC post-doctoral fellowship program (S.S.K.). B.T.Z. is supported by China Scholarship Council (CSC) to conduct research at Penn State University. Monte Carlo simulations in this work were carried out on Cray XC30 at Center for Computational Astrophysics, National Astronomical Observatory of Japan and on the cluster Draco in Tohoku University.

Appendix A Details of Monte Carlo simulations

In this appendix, we describe the results of Monte Carlo simulations, focused on the situation that $R_{\rm coc}/\lambda_{i,\rm coc}<R_{\rm jet}/\lambda_{i,\rm jet}$ and $l_{\rm coh}<R_{\rm jet}$.

A.1 Parameter dependence

As discussed in the main text, we obtain the maximum energy by setting $t_{\rm acc}=t_{\rm esc}$, which results in

[TABLE]

We perform Monte Carlo simulations with various parameter sets tabulated in Table 1 to see the values of $\zeta$. The results are shown in Fig. 5, where the lines represent the escape spectra and the corresponding arrows show the peak energy estimated by Eq. (9) with $\delta=2$ and $\zeta\simeq 2.2(\xi_{c/j}/10)^{0.2}$. We can see that the simulation results agree with the estimates well.

According to our simulation results, the spectral shape is not sensitive to the parameters for $E\gtrsim E_{\rm coh}$, as seen in Fig. 5. We try to fit the spectral shape there using a combination of a power law growth and a cutoff. We consider $EL_{E}\propto E^{a}\exp(-(E/E_{0})^{b})$, and find that $a\sim$ 5–9 and $b\sim$ 0.1–0.3 for the parameter range that we explored. Note that the fitting requires $a>1$, $b<1$, and $E_{0}\ll E_{i,\rm max}$ because of slower cutoff than the exponential.

A.2 Bohm diffusion model

We also perform Monte Carlo simulations using the Bohm limit in the cocoon, $\lambda_{i,\rm coc}=E/(Z_{i}eB_{\rm coc})$ for $E<E_{i,\rm coh}$ and $\lambda_{i,\rm{coc}}=(E/E_{i,\rm{coh}})^{2}l_{\rm{coh}}$ for $E>E_{i,\rm coh}$. Figure 6 shows the escape spectra for the cases with the Bohm limit. For these cases, the maximum energy is represented by Eq. (9), while the spectra for $E<E_{i,\rm coh}$ is harder than those with the Kolmogorov turbulence. This difference arises from the difference of energy dependence of the mean free path. For the particles of $E<E_{i,\rm coh}$, the mean free path for the Bohm limit is shorter than that for the Kolmogorov turbulence. The shorter mean free path leads to the higher return probability, which results in the harder escape spectrum for the Bohm limit cases.

The shorter mean free path in the cocoon also increases the value of $E_{i,\rm inj}$, leading to the lower $\dot{N}_{i,\rm inj}$ and $L_{\rm UHECR}$. To obtain the required $L_{\rm UHECR}$ and $E_{p,\rm max}$, we would need lower $B_{\rm coc}$ and higher $\beta_{\rm jet}$.

Appendix B Application to small-scale jets of radio-quiet AGNs

Our model can be applied to the small-scale jets in radio-quiet AGNs (e.g., Mundell et al. (2003); Ho (2008)), which are energetically possible to create UHE CRs Pe’er et al. (2009). The jet size there is much smaller than that in the radio galaxy. We use $R_{\rm jet}\sim 10$ pc, $l_{\rm jet}\sim 200$ pc, and $R_{\rm coc}\sim 100$ pc. The magnetic field can be stronger than that in the radio galaxy owing to its smaller size. We use $B_{\rm coc}\sim 160\rm~{}\mu$G. This jet is embedded in galactic center, so the CR density and metalicity can be enhanced, compared to those around the Earth. We use 20 times higher CR density Abramowski et al. (2016) and 2 times higher metallicity Rudolph et al. (2006) than those in the Galaxy described in the main text. At the center of the radio-quiet AGNs, the outflows of velocity $\sim$100–1000 km s*-1* are observed Steffen et al. (1997), and we use $v_{\rm exp}\sim 1000$ km s*-1*. We set the other parameters to be the same as those for the radio gaalxy; $\beta_{\rm jet}\sim 0.7$, $R_{\rm sl}\sim 0.01R_{\rm coc}$, and $l_{\rm coh}\sim 0.03R_{\rm coc}$.

Using above parameters, we obtain the source spectrum as shown in Fig. 7. The radio-quiet AGNs can accelerate protons up to a few EeV and irons up to several tens of EeV. The number density of radio-quiet AGNs is around $10^{-3}\rm~{}Mpc^{-3}$ Ho (2008), so the required differential luminosity per source is $3\times 10^{39}\rm~{}erg~{}s^{-1}$. However, we find that the radio-quiet AGN model does not reach the required luminosity. It is difficult for the radio-quiet AGN model to achieve both the required values of $E_{p,\rm max}$ and $L_{\rm UHECR}$.