Ultra-High-energy Cosmic Rays from beyond the Greisen–Zatsepin–Kuz'min Horizon

2.1.1. Protons and Neutrons

Energetic protons interact with ambient radiation fields to produce leptons and pions. These photo-hadronic interactions are dominated by Bethe-Heitler (BH) photo-pair production (Bethe & Heitler 1934) and photopion production processes. Photo-pair processes proceed as

$$\begin{eqnarray}&&A^{\prime} +\gamma \to A+{l}^{+}+{l}^{-},\end{eqnarray} \tag{ 1 }$$

where $A^{\prime}$ and A are nucleons before and after the interaction respectively, while l⁺ and l⁻ are the leptons and anti-leptons produced in this process. As photoelectron pair production contributes most to photo-pair CR energy losses (see e.g., Blumenthal 1970; Klein 2006) over the range of energies of interest in this work, ¹⁰ we consider only the production of electrons (e⁻) and positrons (e⁺). Without loss of generality, hereafter, we do not distinguish between electrons and positrons unless otherwise stated, and both are referred to as "electrons."

Photopion production operates mostly through the following channels: (i) resonant single-pion production, (ii) direct single-pion production, and (iii) multiple-pion production (Mücke et al. 1999). ¹¹ Resonant single-pion production arises through the production of Δ⁺ intermediaries, which decay predominantly through the channels:

$$\begin{eqnarray}&&p+\gamma \to {{\rm{\Delta }}}^{+}\to \left\{\begin{array}{l}p+{\pi }^{0}\to p+2\gamma \\ n+{\pi }^{+}\to n+{\mu }^{+}+{\nu }_{\mu }\\ \,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\unicode{x021B3}\,{{e}}^{+}+{\nu }_{{\rm{e}}}+{\bar{\nu }}_{\mu }\end{array}\right.\end{eqnarray} \tag{ 2 }$$

(see Berezinsky & Gazizov 1993). The decays produce charged and neutral pions, and their branching ratios (BRs) are roughly 2:3 and 1:3, respectively. The pions will further decay, producing cascades of lower energy particles. Multi-pion production occurs at higher energies (see Mücke et al. 1999), and when the invariant interaction energy exceeds _r ≈ 3500, charged and neutral pion production can arise with high multiplicities. For the energy range of interest in this work, each type of pion is produced in roughly equal proportion when all processes and the efficiencies of their production channels are taken into account (see Dermer & Menon 2009).

Calculations for interaction lengths and their corresponding path lengths in blackbody intergalactic radiation fields were presented in Owen et al. (2018) and Owen et al. (2019). The same treatment is also adopted in this work. We note that free neutrons are produced in photo-hadronic interactions, as the one described in Equation (2), and these neutrons will undergo a β⁻-decay:

$$\begin{eqnarray}&&n\to p+{e}^{-}+{\bar{\nu }}_{{\rm{e}}},\end{eqnarray} \tag{ 3 }$$

with a mean life of 880 s (Nakamura 2010). A neutron will decay into a proton, an electron, and an electron anti-neutrino if it does not collide and interact with other particles. In astrophysical environments, before β⁻-decay occurs, a neutron could interact with the radiation field, leading to pion production:

$$\begin{eqnarray}&&n+\gamma \to {{\rm{\Delta }}}^{0}\to \left\{\begin{array}{l}p+{\pi }^{-}\\ n+{\pi }^{0}\end{array}\right..\end{eqnarray} \tag{ 4 }$$

The BRs for Δ⁰ → π⁻ and Δ⁰ → π⁰ are 1:3 and 2:3, respectively, near the energy threshold of the interaction (Hümmer et al. 2010).

The rates of photopion interactions for neutrons and protons do not differ much. Thus, the interaction cross-section and path length for high-energy neutrons can be approximated by the corresponding values for protons. (Hümmer et al. 2010; Romero & Gutiérrez 2020). ¹² Despite the short mean life for neutron decay, neutrons can still traverse a long distance to collide with another particle or a photon because of time dilation resulting from their relativistic motion. However, in this work, we seek to assess the propagation distances and spectral form of UHE CRs based on their nucleon number A, with their charge only being consequential to specify the injection spectrum of each species. We assume that the injection of A = 1 UHE CRs is dominated by protons, while CR neutrons only emerge as the interaction products of primary CRs. As such, we model the interactions and propagation of UHE CR neutrons and protons in the same way and do not distinguish between the two species. ¹³

There are several consequences of the photopion processes. An isotropic cosmogenic neutrino component has been attributed to decaying pions arising from photopion interactions of CRs with intergalactic radiation fields (Allard et al. 2006; Rodrigues et al. 2021). Moreover, electrons produced in a photopion cascade, together with those electrons formed by BH pair production, could up-scatter the ambient low-energy photons to γ-rays, producing a diffuse γ-ray flux. In comparison with the γ-rays produced directly from the decay of neutral pions (a product of the photopion process), this Compton-scattered component is expected to be substantially weaker. ¹⁴

2.1.2. Nuclei

CR nuclei interact with radiation fields via BH pair production, (nuclear) photopion production, and nuclear photospallation. BH pair production proceeds as

$$\begin{eqnarray}&&{N}_{1}^{\prime} +\gamma \to {N}_{1}+{l}^{+}+{l}^{-}\end{eqnarray} \tag{ 5 }$$

(Blumenthal 1970), where ${N}_{1}^{\prime}$ and N₁ are the nuclei before and after the interaction, respectively. The BH process is Z² times faster for nuclei with charge number Z than for protons. Photopion production channels are the same for nuclei as they are for protons (see Section 2.1.1), with corresponding interaction products: for example, nuclear photopion processes can also contribute to cosmogenic neutrino production. However, during their interactions, higher kinetic energies are available for nuclei compared to that for protons. The contribution of protons is typically weighted by a factor of $A={m}_{{{\rm{N}}}_{1}}/{m}_{p}$ (e.g., see Stecker 1979; Dermer & Menon 2009), which implies that heavier nuclei would lose more energy during photopion interactions than lighter nuclei.

The main nuclear photospallation channels are

$$\begin{eqnarray}&&{N}_{1}\to \left\{\begin{array}{l}{N}_{2}+p\\ {N}_{3}+2p\end{array}\right.,\end{eqnarray} \tag{ 6 }$$

where N₁ is the primary nucleus of an atomic mass A, and N₂ and N₃ are the secondary nuclei of atomic masses A − 1 and A − 2, respectively (Stecker 1969; Puget et al. 1976; Stecker & Salamon 1999; Ahlers & Salvado 2011). Note that neutrons can also be emitted instead of protons in this process. Other branches, which remove more than one or two nucleons, may occur in a single interaction. These processes, however, occur with much lower probability and are not included in our calculations.

Photospallation (disintegration) is dominated by the formation of the Giant Dipole Resonance (GDR) intermediary. The process may be considered as a two-step mechanism (Puget et al. 1976), where the first stage involves the photoabsorption by a nucleus and the subsequent stage is the emission of a single or multiple nucleons from the nucleus via a (statistical) decay (Hayward 1960). Supposing that each of these emitted nuclei retain a fraction of 1/A of the primary nucleus in a photospallation event, an approximate relation can be obtained from the energy partition between nucleons as

$$\begin{eqnarray}&&{\epsilon }_{A-1}\approx \left(\displaystyle \frac{A-1}{A}\right)\,{\epsilon }_{A},\end{eqnarray} \tag{ 7 }$$

such that the energy of the produced secondary nuclei can be estimated in an interaction (see Dermer & Menon 2009). Here, _A and _A−1 are, respectively, the dimensionless energies of the nuclei before and after a photospallation event, with _A = E_A/m_e c² = γ_A(m_A/m_e) and a similar expression for _A−1. A consequence of the photospallation of heavy nuclei is the production of lighter nuclei in large numbers. Each of these secondary nuclei interacts with ambient radiation, hence initiating a new chain of photospallation. UHE CR nucleons thus litter traces of particles along their paths when they propagate across the Universe.

The propagation of energetic particles in intergalactic space is described by the transport equation. In the regime appropriate for this study, it may take the following form:

$$\begin{eqnarray}&&\displaystyle \frac{\partial {{ \mathcal N }}_{A}}{\partial \,t}=-{\rm{\nabla }}\cdot {{ \mathcal S }}_{A}+\displaystyle \frac{\partial }{\partial \,{\epsilon }_{A}}\left[\,b\,{{ \mathcal N }}_{A}\,\right]+{{ \mathcal Q }}_{A}-{{\rm{\Lambda }}}_{A}{{ \mathcal N }}_{A}.\end{eqnarray} \tag{ 8 }$$

Here, _A is the energy (of the particles) and b is the rate of energy loss (presumably caused by cooling). The index A denotes the nuclear species. In this work, we relate nuclear charge to mass following the correspondence presented in Puget et al. (1976). The charge of a species is not strictly set unambiguously by its mass. However, in most cases, there is a clear relation between the mass number of a UHE CR and its species/charge, with alternative isotopes being unstable and not being preferentially formed. ¹⁵ We further define ${{ \mathcal N }}_{A}$ as the number density of the particles with mass number A, while ${{ \mathcal S }}_{A}$ is the CR streaming term, which is given by $c{{ \mathcal N }}_{A}$ in the observer's frame if diffusion is negligible (e.g., Webb & Gleeson 1979). ${{ \mathcal Q }}_{A}$ is the source term specifying the rate of particle injection. ${{\rm{\Lambda }}}_{A}{{ \mathcal N }}_{A}$ is the sink term specifying the rate of particle loss, with an efficiency set by the parameter Λ_A . In cosmological settings, we assign ${{ \mathcal N }}_{A}$ as the co-moving number density of a particle species. Therefore, ${{ \mathcal N }}_{A}={n}_{A}{(1+z)}^{-3}$, where n_A is the physical number density of a particle species, and z is the cosmological redshift.

The streaming of UHE CRs is effectively the speed of light, c. Thus, Equation (8) becomes

$$\begin{eqnarray}&&\displaystyle \frac{\partial \,{n}_{A}}{\partial \,t}=-c\displaystyle \frac{\partial \,{n}_{A}}{\partial \,s}+\displaystyle \frac{\partial }{\partial \,{\epsilon }_{A}}\left[\,b\,{n}_{A}\right]+{\left(1+z\right)}^{3}{{ \mathcal Q }}_{A}-{{\rm{\Lambda }}}_{A}\,{n}_{A},\end{eqnarray} \tag{ 9 }$$

with s as the co-moving distance. In terms of the cosmological redshift, this may be expressed as

$$\begin{eqnarray}\begin{array}{rcl}\displaystyle \frac{\partial }{\partial \,t}\ {n}_{A}({\epsilon }_{A},z) & = & -c\,{\left(\displaystyle \frac{{ds}}{{dz}}\right)}^{-1}\,\displaystyle \frac{\partial }{\partial z}\,{n}_{A}({\epsilon }_{A},z)\\ & & +\displaystyle \frac{\partial }{\partial \,{\epsilon }_{A}}\left[\,{b}_{A}({\epsilon }_{A},z)\,{n}_{A}({\epsilon }_{A},z)\,\right]\\ & & +{Q}_{A}({\epsilon }_{A},z)-{{\rm{\Lambda }}}_{A}\,{n}_{A}({\epsilon }_{A},z),\end{array}\end{eqnarray} \tag{ 10 }$$

where ${Q}_{A}={{ \mathcal Q }}_{A}{\left(1+z\right)}^{3}$. In a Friedmann–Lemaître–Robertson–Walker (FLRW) universe,

$$\begin{eqnarray}&&\displaystyle \frac{{ds}}{{dz}}=\displaystyle \frac{c\,{ \mathcal E }(z)}{{H}_{0}\,(1+z)},\end{eqnarray} \tag{ 11 }$$

and

$$\begin{eqnarray}&&{ \mathcal E }(z)={\left[{{\rm{\Omega }}}_{{\rm{r}},\ 0}{\left(1+z\right)}^{4}+{{\rm{\Omega }}}_{{\rm{m}},\ 0}{\left(1+z\right)}^{3}+{{\rm{\Omega }}}_{{\rm{\Lambda }},0}\right]}^{-1/2}\end{eqnarray} \tag{ 12 }$$

(see, e.g., Peacock 1999), where Ω_m,0 = 0.315 ± 0.007, Ω_r,0 ≈ 0, and Ω_Λ,0 = 0.685 ± 0.007 are the normalized density parameters for matter, radiation, and dark energy, respectively. The present value of the Hubble parameter H₀ =100 h km s⁻¹Mpc⁻¹, where h = 0.673 ± 0.006 (Planck Collaboration et al. 2020).

A quasi-steady condition is adopted in our calculations, and this greatly simplifies the procedures involved in solving the particle-transport equations. The quasi-steady condition is justified, as the time of propagation of the particles across length-scales of the order of their horizons is much shorter than the Hubble time. Setting ∂_t = 0, Equation (10) becomes

$$\begin{eqnarray}&&\displaystyle \frac{\partial \,{n}_{A}}{\partial z}=\displaystyle \frac{1}{c}\left[\displaystyle \frac{\partial }{\partial {\epsilon }_{A}}\left({b}_{A}{n}_{A}\right)+{Q}_{A}-{{\rm{\Lambda }}}_{A}{n}_{A}\right]\displaystyle \frac{{ds}}{{dz}},\end{eqnarray} \tag{ 13 }$$

which will be used in Section 2.2.1 for the transport calculations of CR protons and nuclei.

2.2.1. Interaction Rates

The source term in Equation (13) has two components,

$$\begin{eqnarray}&&{Q}_{A}={Q}^{{\rm{a}}}+{Q}^{\mathrm{sp}}.\end{eqnarray} \tag{ 14 }$$

Here, Q^a specifies the rate primary nuclei are supplied by a source population of cosmic accelerators, and is defined by the source population, its evolution (specified by redshift, z), and the CR energy _A and mass number A (our adopted models are discussed in Section 3). Q^sp is the rate of production of secondary nuclei from photospallation of heavy nuclei (i.e., the particle cascade contribution). We consider that the injected particles are represented by ¹H, ⁴He, ¹⁴N, ²⁸Si, and ⁵⁶Fe (with their relative abundance fractions detailed in Section 3). The spallation (disintegration) of heavy CR nuclei produces secondary nuclei with mass numbers ranging from 55 to 1, and the injection rate of secondary nuclei (for mass number less than 56) by photospallation particle cascades is given by

$$\begin{eqnarray}&&{{Q}^{\mathrm{sp}}}_{A}(z,{\epsilon }_{A})=\sum _{A^{\prime} }{\int }_{{\epsilon }_{\min }}^{{\epsilon }_{\max }}d\,{\epsilon }_{{\rm{A}}^{\prime} }\ \displaystyle \frac{{R}_{A}^{A^{\prime} }\,{n}_{A^{\prime} }}{{\left(1+z\right)}^{3}}.\end{eqnarray} \tag{ 15 }$$

Here, ${\epsilon }_{\min }$ and ${\epsilon }_{\max }$ are the minimum and maximum energies considered in our calculations. We set the lowest energy to be ${\epsilon }_{\min }{m}_{{\rm{e}}}{c}^{2}=3.98\times {10}^{18}\,\mathrm{eV}$, as it is not clear if the cosmic-ray spectrum is dominated by extragalactic particles at lower energies than this (Giacinti et al. 2012; Aloisio et al. 2014). The maximum energy is set to be ${\epsilon }_{\max }{m}_{{\rm{e}}}{c}^{2}=3.16\times {10}^{20}\,\mathrm{eV}$, corresponding to the most energetic UHE CRs detected at Earth; see Bird et al. (1995). The summation in ${{Q}^{\mathrm{sp}}}_{A}(z,{\epsilon }_{A})$ accounts for the primary nuclei from $A^{\prime}$ up to A = 56, where the photospallation injection rate for the most massive species under consideration (i.e., ${{Q}^{\mathrm{sp}}}_{56}$) is zero. ${{Q}^{\mathrm{sp}}}_{A}$ is dependent on ${n}_{A^{\prime} }$ (where $A^{\prime} \gt A$). As such, we solve Equation (13) from the CR nuclear species with largest A(=56) to those with the smallest A(=1 for protons/neutrons). ${R}_{A}^{A^{\prime} }({\epsilon }_{A^{\prime} },{\epsilon }_{A},z)$ is the differential rate of the photospallation of species with mass number $A^{\prime}$ into species with mass number A, in the lab frame, in the presence of a soft photon field n_ph(_ph). If the radiation field is considered to be isotropic, this can be expressed as

$$\begin{eqnarray}\begin{array}{rcl}{R}_{A}^{A^{\prime} } & = & \displaystyle \frac{c}{2{\gamma }_{A}^{2}}{\int }_{0}^{\infty }d{\epsilon }_{\mathrm{ph}}\,\left\{\,\displaystyle \frac{{n}_{\mathrm{ph}}({\epsilon }_{\mathrm{ph}})}{{{\epsilon }_{\mathrm{ph}}}^{2}}\right.\\ & & \left.\times {\int }_{0}^{2{\gamma }_{A^{\prime} }{\epsilon }_{\mathrm{ph}}}d{\epsilon }_{{\rm{r}}}\,{\epsilon }_{{\rm{r}}}{\sigma }_{A^{\prime} \to A}({\epsilon }_{{\rm{r}}})\,\right\}\end{array}\end{eqnarray} \tag{ 16 }$$

(see Dermer & Menon 2009). Here, ${\gamma }_{A^{\prime} }$ is the Lorentz factor of the primary nucleus, ${\sigma }_{A^{\prime} \to A}({\epsilon }_{{\rm{r}}})$ is the photospallation cross-section of the primary nucleus of $A^{\prime}$ to produce a secondary nucleus of A, and _r is the total invariant interaction energy (in a dimensionless form). For a blackbody radiation field ¹⁶ , the integrals in Equation (16) can be evaluated analytically. This simplifies the calculations for the region where the CMB dominates the radiation field, or the EBL if it is approximated as a summation of blackbody spectra (see Appendix).

In the two-step mechanism for photospallation, proper accounting for the physical processes in both steps are required in order to obtain an accurate value for the cross-section (Puget et al. 1976). For the case of single nucleon emission, the cross-section can, however, be adequately approximated by a Lorentzian function:

$$\begin{eqnarray}&&{\sigma }_{1}({\epsilon }_{{\rm{r}}})=\displaystyle \frac{{\sigma }_{0,\ {\rm{A}}}{{\epsilon }_{{\rm{r}}}}^{2}{{{\rm{\Delta }}}_{n}}^{2}}{{\left({{\epsilon }_{0}}^{2}-{{\epsilon }_{{\rm{r}}}}^{2}\right)}^{2}+{{\epsilon }_{{\rm{r}}}}^{2}{{{\rm{\Delta }}}_{n}}^{2}}\end{eqnarray} \tag{ 17 }$$

(see Karakula & Tkaczyk 1993; Wang et al. 2008). Photospallation may result in multiple nucleon emission when _r is above 30 MeV or below the pion-production threshold (≃150 MeV), albeit with a much smaller cross-section and lower rate. A more complicated expression is available to describe this variation in the cross-section (see Puget et al. 1976), which is appropriate when the ambient radiation field has a very hard spectrum (see Wang et al. 2008). However, for GZK attenuation calculations, the ambient radiation fields which dominates the interactions are the CMB and EBL, for which the simplified expression for the cross-section is sufficient.

For analytical tractability, this formulation can be simplified further by replacing the Lorentzian with a delta function:

$$\begin{eqnarray}&&{\sigma }_{1}({\epsilon }_{{\rm{r}}})={\sigma }_{0,{\rm{A}}}\,{{\rm{\Delta }}}_{n}\delta ({\epsilon }_{{\rm{r}}}-{\epsilon }_{0})\ \end{eqnarray} \tag{ 18 }$$

(Wang et al. 2008; Dermer & Menon 2009). This still yields a photospallation rate correct to within an order of magnitude (Wang et al. 2008), which is sufficient for the purposes of our demonstrative model. In Equation (18), the width of the cross-section in energy Δ_n = 8 MeV, is given by the energy bandwidth of the giant resonance, and the (dimensionless) energy at which the giant resonance cross-section peaks is given by ₀ = 83.46 A^−0.21 for A > 4, and ₀ = 1.81A^2.433 for A ≤ 4 (Karakula & Tkaczyk 1993). The maximum value of the cross-section is σ_0,A = 1.45 × 10⁻²⁷ A cm² (Wang et al. 2008). Substitution of this cross-section, given in Equation (18), into Equation (16) yields

$$\begin{eqnarray}\begin{array}{rcl}{R}_{1} & = & -\displaystyle \frac{4\pi \,f(z)\,c{\sigma }_{0,{\rm{A}}}{{\rm{\Delta }}}_{n}{\epsilon }_{0}{\rm{\Theta }}(z)}{{{\gamma }_{A}}^{2}{{\lambda }_{c}}^{3}}\\ & \times & \mathrm{ln}\left[\,1-\exp \left(\displaystyle \frac{-{\epsilon }_{0}}{2{\gamma }_{A}{\rm{\Theta }}(z)}\right)\right]\ \end{array}\end{eqnarray} \tag{ 19 }$$

(see also Dermer & Menon 2009, for an expression derived under similar circumstances). Here, λ_c is the Compton wavelength of electron, Θ = k_B T/m_e c² is the dimensionless temperature of the radiation field, and f = 1 for a blackbody radiation field, e.g., the CMB, but f(z) < 1 in the case of diluted blackbody radiation (e.g., redshift-dependent EBL components; see Appendix). This rate can also be used to compute the specific absorption rate of UHE CR nuclei due to photospallation. For this, we adopt the single-nucleon emission approximation of Wang et al. (2008), giving

$$\begin{eqnarray}&&{{{\rm{\Lambda }}}^{\mathrm{sp}}}_{A}\approx \displaystyle \frac{{R}_{1}}{A}.\end{eqnarray} \tag{ 20 }$$

This forms a component of the sink (or absorption) term in the transport equation, i.e., Equation (13).

Adopting the single-nucleon approximation simplifies Equation (15) to

$$\begin{eqnarray}&&{{Q}^{\mathrm{sp}}}_{A}(z,{\epsilon }_{A})\approx \displaystyle \frac{1}{{\left(1+z\right)}^{3}}{\int }_{{\epsilon }_{\min }}^{{\epsilon }_{\max }}d{\epsilon }_{A^{\prime} }{R}_{1}{n}_{A^{\prime} }\end{eqnarray} \tag{ 21 }$$

for nuclei with mass number A > 1, and

$$\begin{eqnarray}\begin{array}{rcl}{{Q}^{\mathrm{sp}}}_{1}(z,{\epsilon }_{1}) & \approx & \displaystyle \frac{1}{{\left(1+z\right)}^{3}}\ \left\{2{\int }_{{\epsilon }_{\min }}^{{\epsilon }_{\max }}d{\epsilon }_{2}{R}_{1}{n}_{2}\right.\\ & & \left.+\displaystyle \sum _{A^{\prime} }{\int }_{{\epsilon }_{\min }}^{{\epsilon }_{\max }}d{\epsilon }_{A^{\prime} }{R}_{1}{n}_{A^{\prime} }\right\}\ \end{array}\end{eqnarray} \tag{ 22 }$$

for A = 1 nuclei. Here we do not distinguish between protons and neutrons (see Section 2.1.1), and the first term in Equation (22) accounts for the photosplitting of A = 2 nuclei into two individual nuclei.

Nuclei can also be absorbed by photopion interactions. The expression of the interaction reaction rate in this case takes the same form as that given by Equation (16), but with the substitution of the pion-production cross-section, σ^{A π} . In this work, we also approximate this with a delta function

$$\begin{eqnarray}&&{\sigma }^{A\pi }=\displaystyle \frac{\pi {\sigma }_{0}{{\rm{\Gamma }}}_{\mathrm{res}}A}{2}\ \delta (\epsilon -{\epsilon }_{\mathrm{res}}),\end{eqnarray} \tag{ 23 }$$

where ${{\rm{\Gamma }}}_{\mathrm{res}}$ is introduced as the resonance width. Adopting the same parameterization of ${\epsilon }_{\mathrm{res}}$ and σ₀ as Unger et al. (2015) ¹⁷ (see also, e.g., Moncada et al. 2017), we obtain ${\epsilon }_{\mathrm{res}}\approx 665$, σ₀ = 5.0 × 10⁻²⁸ cm², and ${{\rm{\Gamma }}}_{\mathrm{res}}=150\,\ \mathrm{MeV}$. By substituting Equation (23) into Equation (16), we obtain an analytical expression for the interaction rate:

$$\begin{eqnarray}\begin{array}{rcl}{{R}^{A\pi }}_{A} & = & \displaystyle \frac{4\pi \,f\,c}{{{\gamma }_{A}}^{2}{{\lambda }_{c}}^{3}}{\int }_{0}^{\infty }d{\epsilon }_{\mathrm{ph}}\,\left\{\displaystyle \frac{1}{\exp \left[{\epsilon }_{\mathrm{ph}}/{\rm{\Theta }}\right]-1}\right.\\ & & \left.\times {\int }_{0}^{2{\gamma }_{A}{\epsilon }_{\mathrm{ph}}}d{\epsilon }_{{\rm{r}}}\,{\epsilon }_{{\rm{r}}}\pi {\sigma }_{0}\displaystyle \frac{{{\rm{\Gamma }}}_{\mathrm{res}}}{2}A\,\delta ({\epsilon }_{{\rm{r}}}-{\epsilon }_{\mathrm{res}})\right\}\\ & = & -\displaystyle \frac{2{\pi }^{2}\,f\,c{\sigma }_{0}{{\rm{\Gamma }}}_{\mathrm{res}}A{\epsilon }_{\mathrm{res}}{\rm{\Theta }}}{{{\gamma }_{A}}^{2}{{\lambda }_{c}}^{3}}\\ & & \times \mathrm{ln}\left[1-\exp \left(\displaystyle \frac{-{\epsilon }_{\mathrm{res}}}{2{\gamma }_{A}{\rm{\Theta }}}\right)\right],\end{array}\end{eqnarray} \tag{ 24 }$$

where the last line expresses the rate in a similar form to that for photospallation, which is given in Equation (19). The corresponding specific absorption rate of CRs due to photopion interactions then follows simply as

$$\begin{eqnarray}&&{{{\rm{\Lambda }}}^{A\pi }}_{A}=\displaystyle \frac{{{R}^{A\pi }}_{A}}{A},\end{eqnarray} \tag{ 25 }$$

such that the overall sink term for nuclei, for use in Equation (13), may be written as

$$\begin{eqnarray}&&{{\rm{\Lambda }}}_{A}={{{\rm{\Lambda }}}^{\mathrm{sp}}}_{A}+{{{\rm{\Lambda }}}^{A\pi }}_{A}.\end{eqnarray} \tag{ 26 }$$

The other CR interactions described in Section 2.1 typically operate as continuous loss processes, or such that only a small fraction of the particle's energy is transferred in a single interaction. These can be modeled as cooling processes and the total cooling rate experienced by a UHE CR, in terms of the dimensionless energy of the nucleon (again, we retain our earlier convention where energies are in units of electron rest mass), therefore consists of three components:

$$\begin{eqnarray}&&b=-d\epsilon /{dt}={{b}^{\gamma A}}_{A}+{b}_{\mathrm{rad}}+{b}_{\mathrm{ad}},\end{eqnarray} \tag{ 27 }$$

where ${{b}^{\gamma A}}_{A}$ is the cooling rate due to photo-pair production, b_rad is the radiative cooling rate, and b_ad is the adiabatic cooling rate of UHE CR nuclei as they propagate through an expanding cosmology.

In our calculations for photo-pair production, the invariant dimensionless interaction energy _r ≳ 60, so a fitted approximation for the inelastic cross-section of photo-pair losses is available (see Jost et al. 1950; Bethe & Maximon 1954; Blumenthal 1970; Stepney & Guilbert 1983):

$$\begin{eqnarray}&&{\sigma }^{\gamma A}({\epsilon }_{{\rm{r}}})\approx \left\{\displaystyle \frac{7}{6\pi }\,{\alpha }_{{\rm{f}}}\mathrm{ln}\left[\displaystyle \frac{{\epsilon }_{{\rm{r}}}}{{k}_{\gamma e}}\right]\right\}{\sigma }_{{\rm{T}}},\end{eqnarray} \tag{ 28 }$$

where α_f is the fine-structure constant, σ_T is the Thomson cross-section, and k_{γ e} is a fitting constant, which we set as 2.0 (see Owen et al. 2018). The CR cooling rate due to photo-pair losses can then be expressed as

$$\begin{eqnarray}\begin{array}{rcl}{{b}^{\gamma A}}_{A} & = & -\displaystyle \frac{d\,{\epsilon }_{A}}{d\,t}\\ & = & -\displaystyle \frac{f(z)\,{Z}_{A}^{2}\,c}{2{\gamma }_{A}A}{\int }_{{\gamma }_{A}^{-1}}^{\infty }d{\epsilon }_{\mathrm{ph}}\ \left\{\displaystyle \frac{1}{\exp \left[{\epsilon }_{\mathrm{ph}}/{\rm{\Theta }}(z)\right]-1}\right.\\ & & \left.\times {\int }_{2}^{2{\gamma }_{A}{\epsilon }_{\mathrm{ph}}}\displaystyle \frac{d{\epsilon }_{{\rm{r}}}\,{\epsilon }_{{\rm{r}}}{\sigma }^{\gamma A}({\epsilon }_{{\rm{r}}})}{\sqrt{1\,+\,2{\epsilon }_{{\rm{r}}}}}\right\}\end{array}\end{eqnarray} \tag{ 29 }$$

(Dermer & Menon 2009; Aloisio et al. 2013a), where the lower boundary value of the inner integral follows from at least one pair of electrons being formed by the process. As in Owen et al. (2019), we assume that the electron–positron pairs are produced in the zero-momentum frame of the interaction, and the interaction energy is completely dominated by the contribution from the interacting nucleus (see also Protheroe & Johnson 1996; Dermer & Menon 2009). This yields

$$\begin{eqnarray}&&{{b}^{\gamma A}}_{A}({\epsilon }_{A})\approx -\displaystyle \frac{f\,{Z}_{A}^{2}\,c{\alpha }_{{\rm{f}}}{\sigma }_{{\rm{T}}}{{ \mathcal F }}_{\gamma e}(u)}{{{\lambda }_{{\rm{c}}}}^{3}{{\gamma }_{A}}^{2}A\,{u}^{5/2}},\end{eqnarray} \tag{ 30 }$$

where the dimensionless variable $u={({\gamma }_{A}{\rm{\Theta }})}^{-1}$, and the function ${{ \mathcal F }}_{\gamma e}(u)$ takes the same form as that in Owen et al. (2019).

Radiative cooling is dominated by Compton scattering with the ambient radiation field. It arises at a rate of

$$\begin{eqnarray}&&{b}_{\mathrm{rad}}({\epsilon }_{A},z)=-\displaystyle \frac{4}{3}\displaystyle \frac{{\sigma }_{{\rm{T}}}\,U(z)}{{m}_{{\rm{e}}}c}{\left[\displaystyle \frac{{{\epsilon }_{A}}^{1/2}{Z}_{A}}{A}\left(\displaystyle \frac{{m}_{{\rm{e}}}}{{m}_{{\rm{p}}}}\right)\right]}^{4}\ \end{eqnarray} \tag{ 31 }$$

(see the expression given in Puget et al. 1976), where U(z) is the energy density of the radiation field (at redshift z). Note that, here, b_rad is dependent on the charge of the nuclei Z_A .

Adiabatic losses due to cosmological expansion occur at a rate

$$\begin{eqnarray}&&{b}_{\mathrm{ad}}({\epsilon }_{A},z)=-{\epsilon }_{A}\,{H}_{0}\ { \mathcal E }(z)\end{eqnarray} \tag{ 32 }$$

(see Gould 1975; Berezinsky & Grigor'eva 1988; Berezinsky et al. 2006), where ${ \mathcal E }(z)$ is given by Equation (12).

We compare the relative importance of each of these cooling and interaction processes in Figure 1, where the effective path length of a UHE CR nucleus (ℓ_A ) is defined as the characteristic distance over which the CR loses its energy. This can be estimated from the cooling rate b using ${{\ell }}_{A}=c\,{[| b| /{\epsilon }_{A}]}^{-1}$. In Figure 1, we show the effective loss lengths for UHE CR nuclei (A = 1 for protons, A = 4 for ⁴He nuclei, and A = 56 for ⁵⁶Fe nuclei) in intergalactic conditions at z = 0, with continuous losses due to pair production and adiabatic losses shown by the solid lines. These are compared to conditions at z = 3 in Figure 2, where the background radiation fields are more intense (leading to faster photo-pair losses and higher photopion and photospallation rates, but lower adiabatic loss rates due to the relatively lower rate of cosmic expansion at z = 3 compared to the present Universe). Radiative losses are inconsequential in all cases. For the purposes of comparison, we also estimate corresponding effective path lengths for UHE CRs due to processes that we regard as absorption in our calculations, where a substantial amount of energy is lost by the CR in a single interaction. These are photopion production and photospallation; note that photospallation only affects nuclei with A > 1, and so protons (A = 1) are only affected by absorption through photopion interactions. These effective path lengths are distinguished as dashed lines in Figure 1, and are estimated by approximating these stochastic absorption process as a continuous loss process. This allows an equivalent cooling rate to be estimated as ${{b}^{A\pi }}_{A}\approx -{\epsilon }_{A}{{{\rm{\Lambda }}}^{A\pi }}_{A}$ for pion production, or as ${{b}^{\mathrm{sp}}}_{A}\approx -{\epsilon }_{A}{{{\rm{\Lambda }}}^{\mathrm{sp}}}_{A}$ for photospallation. ¹⁸ Although approximated as continuous losses for the comparative purpose of Figures 1 and 2, our later calculations fully account for photopion and photospallation processes as absorption.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Four source populations are considered, and defined to a maximum redshift ${z}_{\max }=3$. This maximum redshift is chosen because it would be extremely unlikely that UHE CRs could reach us from greater distances, given their suppression by photo-hadronic interactions with the CMB and radiation emitted from astrophysical objects during their propagation in intergalactic space. Each of the source population models we consider follows a different redshift evolution. The first model (referred to as the star formation rate, or SFR, model) is based on the redshift evolution of cosmic star formation (also adopted by Muzio et al. 2019), and takes the form:

$$\begin{eqnarray}&&{\psi }_{\mathrm{SFR}}(z)={\psi }_{\mathrm{SFR}}^{0}\,\displaystyle \frac{{\left(1+z\right)}^{{k}_{1}}}{1+{[(1+z)/{k}_{2}]}^{{k}_{3}}},\end{eqnarray} \tag{ 34 }$$

where k₁ = 2.7, k₂ = 2.9, and k₃ = 5.6 are empirical fit parameters inherited from the best-fit function of Madau & Dickinson (2014) (see also Robertson et al. 2015) ¹⁹

Our second model choice is a parameterized redshift distribution of gamma-ray bursts (and is referred to as the GRB model hereafter). It follows the form presented in Wang et al. (2011):

$$\begin{eqnarray}&&{\psi }_{\mathrm{GRB}}(z)={\psi }_{\mathrm{GRB}}^{0}\,{\left(1+z\right)}^{{k}_{4}}\,{\psi }_{\mathrm{SFR}}(z),\end{eqnarray} \tag{ 35 }$$

where k₄ = 1.4. This is based on indications from Swift observations that the GRB rate is enhanced relative to the SFR at high-redshift (Le & Dermer 2007; see also Yüksel & Kistler 2007) ²⁰

We introduce a third model to describe an AGN source population (referred to as the AGN source model hereafter). For this, we adopt the parameterization of Hasinger et al. (2005) (see also Ahlers et al. 2009; Wang et al. 2011; Alves Batista et al. 2019a),

$$\begin{eqnarray}{\psi }_{\mathrm{AGN}}(z)={\psi }_{\mathrm{AGN}}^{0}\left\{\begin{array}{ll}{\left(1+z\right)}^{{k}_{5}} & (z\lt {z}_{1})\ \\ {{z}_{2}}^{{k}_{5}} & ({z}_{1}\leqslant z\lt {z}_{2})\ \\ {{z}_{2}}^{{k}_{5}}\cdot {{z}_{2}}^{{z}_{2}-z} & (z\geqslant {z}_{2})\ \end{array}\right.,\end{eqnarray} \tag{ 36 }$$

where k₅ = 5.0, z₁ = 1.7, z₂ = 2.7, which is largely consistent with the results in the later survey studies on AGN populations, e.g., Silverman et al. (2008), Ajello et al. (2012) (which were both adopted by Jacobsen et al. 2015, in the study of spectral and evolutionary properties of UHE neutrinos), and Ueda et al. (2014). There are various uncertainties in modeling other AGN contributions to the UHE CR flux. While it is unlikely that the results obtained from our parameterized model would need substantial revision, these uncertainties, e.g., those concerning the AGN evolution, could have detectable effects on the spectral properties of UHE CRs at Earth.

As the origins of UHE CRs are still to be determined, some studies adopted a simple power-law parameterization of the redshift evolution of the source model, in the form ${(1+z)}^{{k}_{\mathrm{PLW}}}$, where k_PLW is a fit parameter to be determined (e.g., Taylor et al. 2015; Aab et al. 2017; Jiang et al. 2021). Although our other three source models are more physically motivated, for comparison, we also consider such a power-law parameterization (referred to as the PLW model), defined as

$$\begin{eqnarray}&&{\psi }_{\mathrm{PLW}}(z)={\psi }_{\mathrm{PLW}}^{0}\,{\left(1+z\right)}^{{k}_{\mathrm{PLW}}},\end{eqnarray} \tag{ 37 }$$

with k_PLW = − 1.6, the best-fit value obtained by Alves Batista et al. (2019a). In Equations (34)–(37), ${\psi }_{\mathrm{SFR}}^{0}$, ${\psi }_{\mathrm{GRB}}^{0}$, ${\psi }_{\mathrm{AGN}}^{0}$, and ${\psi }_{\mathrm{PLW}}^{0}$ are the normalization constants computed by integrating the respective source distribution over redshift up to a cutoff (${z}_{\max }=3$). Their resulting values are presented in Table 1, and the redshift evolutionary histories for our four source models are shown in Figure 3.

Zoom In Zoom Out Reset image size

Table 1. Summary of Parameter Choices Adopted for Each of the Four Redshift Source Distribution Models

Model	Normalization ${\psi }_{{\rm{x}}}^{0}$	Spectral index α	$\mathrm{log}({R}_{\max }/V)$	ρCR/1048 erg Mpc−3 yr−1
(1) SFR	${\psi }_{\mathrm{SFR}}^{0}=0.054$	−1.3	18.2	0.5
(2) GRB	${\psi }_{\mathrm{GRB}}^{0}=0.013$	−1.5	18.2	2.0
(3) AGN	${\psi }_{\mathrm{AGN}}^{0}=0.0041$	−1.0	18.2	0.04
(4) PLW	${\psi }_{\mathrm{PLW}}^{0}=1.1$	+ 1.0	18.7	15.0

Note. See text for the choices of the model parameters.

Download table as:  ASCIITypeset image

The normalization term ${{ \mathcal C }}_{0}$ in Equation (33) can be determined from the UHE CR luminosity density, ρ_CR. This is the rate of UHE CR energy density generated by the source population and, at z = 0, it may be expressed as

$$\begin{eqnarray}&&{\rho }_{\mathrm{CR}}={{ \mathcal C }}_{0}\sum _{A}{\int }_{{\epsilon }_{\min }}^{{\epsilon }_{\max }}d\,{\epsilon }_{A}\,{\epsilon }_{A}\,{\psi }_{{\rm{x}}}(z)\,\displaystyle \frac{d}{d\,{\epsilon }_{A}}\left(\displaystyle \frac{d\,{n}^{\star }}{d\,t}\right){\Space{0ex}{2.65ex}{0ex}| }_{z=0}\end{eqnarray} \tag{ 38 }$$

(see Jiang et al. 2021). The exact value of ρ_CR is poorly constrained, however it may be roughly estimated by dividing the measured UHE CR energy density at z = 0, around 10⁵⁴ erg Mpc⁻³ dex⁻¹ (see Aab et al. 2020), by a typical CR energy-loss time, ℓ/c ∼ 1 Gpc/c ∼ 1 Gyr (see Figure 1), giving ∼10⁴⁵ erg Mpc⁻³ yr⁻¹. This is comparable with more rigorous calculations (e.g., Berezinsky et al. 2006; Aloisio et al. 2007; Wang et al. 2011; Jiang et al. 2021), although the exact number is dependent on assumed model parameters, in particular, the lower energy cutoff of the CR spectrum, and the redshift distribution of the sources (Wang et al. 2011; Jiang et al. 2021). This variation is evident, as shown in the studies of Berezinsky et al. (2006) and Aloisio et al. (2007), where substantially larger values than our estimate, of 3.5 × 10⁴⁶ erg Mpc⁻³ yr⁻¹ and 7 × 10⁴⁶ erg Mpc⁻³ yr⁻¹ (at energies above 10¹⁸ eV), respectively, are derived by assuming that there is no redshift dependence in the UHE CR source populations. By contrast, much lower values, in the range 0.2–2 × 10⁴⁴ erg Mpc⁻³ yr⁻¹ (at 3.16 × 10¹⁹ eV), were obtained by Jiang et al. (2021), who adopted a negative redshift evolution for the UHE CR sources population. In Wang et al. (2011), a range of values depending on the assumed source model and energy cutoff are presented: for energies above 10¹⁸ eV, they found the source luminosity density to be ∼0.3 − 1 × 10⁴⁵ erg Mpc⁻³ yr⁻¹, for a source redshift distribution following that in Yüksel et al. (2008) (corresponding to the SFR scenario in this work); ∼6 × 10⁴⁴ erg Mpc⁻³ yr⁻¹ assuming a redshift distribution following the AGN scenario as in this work, cf Equation (36); and ∼0.8 − 3 × 10⁴⁵ erg Mpc⁻³ yr⁻¹ assuming a source distribution following that of Yüksel & Kistler (2007) (cf. the GRB source scenario considered here). While we do not perform detailed model fitting, we set the UHE CR luminosity density for each of the source models to give a reasonable flux at z = 0 when compared to Pierre Auger Observatory data (Verzi 2019); see Figures 5 and 6. The adopted values (at energies above 5 × 10¹⁸ eV) are shown in Table 1 for each of the source models, and these were used to set the normalization ${{ \mathcal C }}_{0}$ in Equation (38) above.

We model the energy spectrum of the injected CRs using the spectral form:

$$\begin{eqnarray}\displaystyle \frac{d\,{n}^{\star }}{d\,{\epsilon }_{A}}\propto {f}_{A}{\left(\displaystyle \frac{{\epsilon }_{A}}{{\epsilon }_{\min }}\right)}^{-\alpha }\left\{\begin{array}{ll}1 & (x\lt 1)\\ \exp \left(1-x\right) & (x\geqslant 1)\ \end{array}\right.\end{eqnarray} \tag{ 39 }$$

(Taylor et al. 2015; Alves Batista et al. 2019a; Jiang et al. 2021), where $x={\epsilon }_{{\rm{A}}}{m}_{{\rm{e}}}{c}^{2}/{Z}_{A}{R}_{\max }$, ${R}_{\max }(V)$ is the rigidity and the energy limits ${\epsilon }_{\min }{m}_{{\rm{e}}}{c}^{2}$ and ${\epsilon }_{\max }{m}_{{\rm{e}}}{c}^{2}$ retain their earlier values. Our parameter choices for the spectral index α and rigidity are shown in Table 1. Our adopted values are comparable to those derived from fits to Pierre Auger Observatory data in the analysis of Alves Batista et al. (2019a). We note that α is strongly dependent on the adopted source population's redshift evolution (see Taylor et al. 2015), leading to substantial differences in our choice of parameter value between the four models.

We adopt a simplified CR composition source model, where the full range of injected species are represented by the abundances of ¹H, ⁴He, ¹⁴N, ²⁸Si, and ⁵⁶Fe. The individual component spectral forms are indicated in Figure 4, and their total contribution is normalized such that f₁ + f₄ + f₁₄ + f₂₈ + f₅₆ = 1. The fitted composition fraction values obtained by Alves Batista et al. (2019a) are adopted for species abundances for each of the four source populations (SFR, GRB, AGN, and PLW). The CR spectral models used in Alves Batista et al. (2019a) and this work are slightly different, but the fitted values given in Alves Batista et al. (2019a) are still a reasonable choice for our model. Note that we calculate the production of all secondary nuclear CR species between ¹H and ⁵⁶Fe formed as the primary UHE CRs propagate and interact with cosmological radiation fields. The abundance fractions of the species used in our calculations are shown in Table 2.

Zoom In Zoom Out Reset image size

The CR spectrum, expected to be observed at z = 0, was computed by solving the transport equation (Equation (13) in Section 2.2) numerically using a Runge–Kutta integrator. ²¹ Cases for the four source populations, SFR, GRB, AGN, and PLW, were calculated, with the numerical integration proceeding from ${z}_{\max }=3$ to ${z}_{\min }=0$, with an explicit treatment of the injection of the CR compositions and of the evolution of the source populations along z. The computed CR flux and composition at z = 0 took account of CR propagation and interaction effects, including cooling, absorption, and disintegration/spallation (outlined in Section 2). The results are are presented in Figure 5.

Zoom In Zoom Out Reset image size

The production of secondary nuclei by photospallation when CRs propagate in intergalactic space drives the evolution of the composition of CRs, with heavy nuclei gradually being eroded to spawn lighter nuclear species, i.e., those of smaller atomic numbers. When the CRs are injected in our model, they are comprised of five nuclear species (i.e. ¹H, ⁴He, ¹⁴N, ²⁸Si, and ⁵⁶Fe) as specified by the source models (see Section 3.4). Other species are produced via spallation as the CRs propagate from their sources to z = 0. This evolution is captured in Figure 5, which shows the compositions in four broad bands to represent the binned distribution of nuclear species when the CRs complete their journey to reach z = 0.

When setting the abundance of ⁵⁶Fe to be zero, the heaviest species injected by the source model (see Section 3) is ²⁸Si (A = 28), so no CR nuclei would have A > 28. As the CRs propagate, nuclei initially having A = 28 would become nuclei with A ≤ 28. The four broadband composition spectra in Figure 5 and the corresponding four injection species spectra in Figure 4 have strikingly similar morphology and relative strengths. The flux in the band 1 ≤ A ≤ 2 is marginally higher, relatively, than the injected H fraction by the source model, while that in the 3 ≤ A ≤ 6, 7 ≤ A ≤ 19, and 20 ≤ A ≤ 39 is marginally lower than the injected fractions of ⁴He, ¹⁴N, and ²⁸Si, respectively. At z = 0, the flux in each band is dominated by the injected source species, indicating that the evolution of species in an ensemble of CR nuclei propagating over vast distances through intergalactic space is only marginal (subject to the effects of cosmological expansion, cooling, absorption, and photospallation). We may therefore conclude that, in the absence of other additional factors, CR species arriving at z = 0 are overwhelmingly primary CRs rather than secondaries reprocessed by photospallation. This result is not particularly surprising given that the injection process persists to z = 0, and the rate of photospallation interactions is not comparable to the injection of primary CRs by the sources for the four populations considered. It also explains the limited evolution seen between the injected source spectrum and the post-propagation spectrum, with cooling rates being easily overcome by new injections at all energies (see Figure 1, which also clarifies that absorption processes dominate over any cooling process for CRs with energy above 10²⁰ eV).

Our calculations for all considered models yield similar spectra at z = 0 to that obtained by observations, and this agreement is evident (see Figure 5) when comparing the computed total flux spectra with data obtained by the Pierre Auger Observatory (PAO) (Verzi 2019). Figure 6 shows a comparison between the total flux spectra obtained for the four models, SFR, GRB, AGN, and PLW. The AGN source population model was found to produce (summed) spectra best matching the observation, while the flux levels in all four model cases are similar to the PAO observations and other studies (e.g., Alves Batista et al. 2019a). We note that an improved fit cannot be achieved with the parameter range considered in our models. Model refinement to obtain a better fit with available data is left to future studies.

Zoom In Zoom Out Reset image size

Figure 7 further compares the average UHE CR spectral mass composition, $\langle \mathrm{ln}A\rangle$, for the four source classes with data (Pierre Auger Collaboration 2013; Yushkov 2019). This shows a clearer distinction between the different source classes, and also clearly rules out the PLW model as parameterized with Equation (37). The AGN model is also disfavored in this comparison (particularly at around E_A ∼ 10¹⁹ eV), however possible variations in the adopted hadronic model, ²² or the possibility of combined source classes with a refined AGN-type component may offer scope for an improved fit.

Zoom In Zoom Out Reset image size

Anisotropy in the arrival directions of UHE CRs above 8 × 10¹⁸ eV, at a significance > 5.2σ was recently reported (see Abreu et al. 2010; Aab et al. 2014, 2015, 2018; Pierre Auger Collaboration et al. 2017; Abbasi et al. 2018). There are also indications of correlations with certain candidate high-energy cosmic accelerators (see Pierre Auger Collaboration et al. 2007; Abreu et al. 2010), which include AGN detected in γ-rays and X-rays (Nemmen et al. 2010; Terrano et al. 2012), ²³ as well as source distributions tracing large-scale structures, with a source population being biased with respect to the distribution of galaxies (Kashti & Waxman 2008). A CR "hotspot" was recently reported by Abbasi et al. (2014). Studies by Fang et al. (2014) and He et al. (2016) showed that this hotspot is consistent with a single source rather than a chance signal. Clusters of excess arrivals of UHE CRs above an isotropic background from the direction of Centaurus A ²⁴ have been found in some analyses (Fargion & D'Armiento 2011; Kim 2013a, 2013b). ²⁵ However, there is no clear detection of excesses around other nearby radio galaxies, e.g., Centaurus B and Virgo A (see Fraija et al. 2019; Kobzar et al. 2019).

If UHE CRs are produced by sources located nearby to our Galaxy with a fraction originating from within the conventional GZK horizon, as well as sources at large (cosmological) distances, then their arrival flux could have two components: an isotropic component corresponding to unresolved distant sources and an anisotropic component associated with relatively nearby sources. The reported anisotropies detected in UHE CRs are consistent with this scenario, and consideration of a foreground and background component in the UHE CR flux was shown in Kim & Kim (2011) to improve anisotropy correlations with the positions of potential nearby sources. This is complimented by the blind test performed by Rubtsov et al. (2012) for verifying the hypothesis that AGN were the origin of UHE CRs, using event sets from Yakutsk, AGASA, and HiRes. Their findings were consistent with a random background (see also High Resolution Fly'S Eye Collaboration et al. 2008, which also found 13 UHE CR events detected by HiRes in the northern hemisphere to be consistent with isotropy).

Here, we assess the fraction of UHE CR flux that could originate at large distances, from sources located beyond a GZK "horizon" distance of a few tens of megaparsecs (see, e.g., Kachelrieß & Semikoz 2019), and determine their contribution to a UHE CR background. Distant CR sources are homogeneously distributed across the celestial sphere over a range of redshifts. Regardless of the fraction of their flux that can survive to reach Earth, CRs originating from these sources should not, at least in principle, show resolvable anisotropy. To assess the relative contributions to the CRs observed at z = 0 from the CRs initially produced by distant sources at different redshifts, we computed the flux spectra (for each of the source models) with several assigned ${z}_{\max }$, taking values from 3.0, 2.0, 1.0, and 0.5. The spectra computed at z = 0 for four source models are presented in Figure 8 for four values of ${z}_{\max }$. This shows that changing either the source model or ${z}_{\max }$ leads to discernible changes in the spectra at z = 0. The contribution of distant source populations is clearly substantial and naturally accounts for the emergence of an isotropic (background) component in the observed UHE CR flux on Earth.

Zoom In Zoom Out Reset image size

The significance of the distant sources' contribution varies among the different source evolution models, which can be seen from comparing the set of spectra of different ${z}_{\max }$ shown in different panels of Figure 8. The amount of reduction between the spectrum from sources below z = 3 to that from sources below z = 2 accounts for a large fraction of the total CR flux arriving at z = 0, which is at a level between 3% and 50% in the SFR, GRB, and AGB models. The PLW model however, has a relatively smaller reduction in the flux for the same redshift interval, implying that CR flux is dominated by sources at lower redshifts in that case.

Figure 9 shows the fractional contribution from sources above a given redshift, f_CR(>z), for the four source models. Most noticeably, the curves for the GRB and AGN models are practically indistinguishable. The morphology of the f_CR(>z) curve is a manifestation of the redshift locations where a specific population of sources has the most significant contribution. The dominant redshift locations of the GRB and AGN CR sources are very similar to one another (cf. the curves of these two source populations in Figure 3, particularly above (z + 1) ≈ 2.5). This, together with the lacking of substantial numbers of GRB and AGN below (z + 1) ∼ 2.5 leads to their almost identical f_CR(>z) curves. Figure 9 indicates that, for all the source population models considered, a large faction of UHE CR flux observed at Earth would originate from regions far beyond the conventional GZK horizon of a few tens of megaparsecs, and would not be attributed to local sources. If the UHE CRs are predominantly produced as a consequence of star formation, which is represented by the SFR model here, then ∼50% of the CRs could originate from as far as z ∼ 2 or above.

Zoom In Zoom Out Reset image size

These findings are robust. The results obtained from our calculations are not sensitive to the assumed composition of the injected UHE CR flux. The outcome is similar even if the CR injections are restricted to a single species (of any of the five nuclei considered). This is in line with the chemical evolution due to photospallation processes being insignificant, otherwise, the composition evolution of CRs would be strongly manifested in the CR energy spectra of the species observed at Earth.