"Super GZK" Particles in a Classic Kramers' Diffusion-over-a-barrier Model. I. The Case of Protons

The Fokker–Planck drift (or systematic energy-loss rate, ${{ \mathcal K }}_{1}(E)$) and diffusion (energy-loss dispersion rate, ${{ \mathcal K }}_{2}(E)$) are derived from the measured (Genzel et al. 1973) double-differential cross-sections $d{\sigma }_{i}^{}/d{\rm{\Omega }}$, where i refers to one of three main reaction channels: p + γ_CMB → π⁺ + n, p + γ_CMB → π⁰ + p, and n + γ_CMB → π⁻ + p, and where such measurements are made in the rest frame of the proton. In the laboratory frame, the energy-loss rate is calculated as

$$\begin{eqnarray}&&\displaystyle \frac{\langle {\rm{\Delta }}E\rangle }{{\rm{\Delta }}s}(E)=\int \int \displaystyle \frac{d\sigma }{{dE}^{\prime} }(E^{\prime} ,E,{E}_{\gamma })\,\rho ({E}_{\gamma })(E-E^{\prime} )\,{dE}^{\prime} {{dE}}_{\gamma },\end{eqnarray} \tag{ 1 }$$

where ${\rm{\Delta }}E=E-E^{\prime}$, $d\sigma /{dE}^{\prime}$ is the energy-loss cross-section for each reaction channel, which is a function of $E^{\prime} (\lt E)$, the energy of the proton after its collision with the CMB photon, the proton's initial energy E, and the CMB photon's energy E_γ . ρ(E_γ ) is the CMB photon energy density. (Note that since c, speed of light, is also essentially the relative velocity between the proton and CMB photon in the laboratory frame, the relative velocity term in Equation (1) above cancels out when we write Δs = c Δt, with t being the proper time).

Following Hill & Schramm (1985), we use temperature-scaled, nondimensionalized energy units, where we define, $\epsilon =E\,{k}_{{\rm{B}}}T/{m}_{o}^{2}$, $\epsilon ^{\prime} =E^{\prime} \,{k}_{{\rm{B}}}T/{m}_{o}^{2}$, and _γ = E_γ /k_B T, ${\rm{\Delta }}\epsilon =\epsilon \,-\epsilon ^{\prime}$, with k_B being Boltzmann's constant, m_o being the nucleon's rest mass, and T being the CMB photons' blackbody temperature. Kinematically we can relate $d\sigma /{dE}^{\prime}$ in Equation (1) above to the measured double-differential cross-section d σ^(cm)/dΩ in the center-of-mass frame, in which frame measurements are also expressed. In these scaled energy units, the measured, parameterized double-differential cross is written as (Genzel et al. 1973; Hill & Schramm 1985):

$$\begin{eqnarray}&&\displaystyle \frac{d\sigma }{d\epsilon ^{\prime} }=\sum _{n=0}^{N}\,{a}_{n}\,\left(\displaystyle \frac{2}{\epsilon {{\prime} }^{+}-{\epsilon }^{-}}\right)\,{\left[\displaystyle \frac{2\epsilon -(\epsilon {{\prime} }^{+}+\epsilon {{\prime} }^{-})}{\epsilon {{\prime} }^{+}-\epsilon {{\prime} }^{-}}\right]}^{n},\end{eqnarray} \tag{ 2 }$$

where a_n are the expansion coefficients in the cosine of the pion's angle in the center-of-mass system, and $\epsilon {{\prime} }^{+}$ ($\epsilon {{\prime} }^{-}$) are upper (lower) kinematic limits for the proton (or neutron) in the same frame. These limits were derived by Hill & Schramm 1985) as:

$$\begin{eqnarray}\begin{array}{rcl}\epsilon {{\prime} }^{\pm } & = & \epsilon \times \left\{\displaystyle \frac{2\,\epsilon \,{\epsilon }_{\gamma }+1-{m}_{\pi }^{2}/(2{m}_{o}^{2})}{4\,\epsilon \,{\epsilon }_{\gamma }+1}\right.\\ & & \pm \left.\displaystyle \frac{{\left[{\left(2\epsilon {\epsilon }_{\gamma }-{m}_{\pi }^{2}/(2{m}_{o}^{2})\right)}^{2}-{m}_{\pi }^{2}/{m}_{o}^{2}\right]}^{1/2}}{4\,\epsilon \,{\epsilon }_{\gamma }+1}\right\},\end{array}\end{eqnarray} \tag{ 3 }$$

where m_π is the pion's rest mass. In the same scaled energy units, the CMB photons' energy density we use (Hill & Schramm 1985) is written as

$$\begin{eqnarray}&&\rho ({\epsilon }_{\gamma })=\displaystyle \frac{{n}_{\gamma }}{2{\pi }^{2}}\,{\epsilon }_{\gamma }^{2}{\int }_{1}^{\infty }\displaystyle \frac{\zeta }{\exp (\zeta \,{\epsilon }_{\gamma }/2)-1}d\zeta ,\end{eqnarray} \tag{ 4 }$$

with n_γ being the CMB photon number density.

In such units, the energy-loss rate is now written as

$$\begin{eqnarray}&&\begin{array}{l}\displaystyle \frac{\langle {\rm{\Delta }}\epsilon \rangle }{{\rm{\Delta }}s}(\epsilon )\\ \quad =\ \ \mathrm{const}{\int }_{{\epsilon }_{\gamma ,\min }}^{{\epsilon }_{\gamma ,\max }}{\int }_{\epsilon {{\prime} -}^{}}^{\epsilon {{\prime} +}^{}}\displaystyle \frac{d\sigma }{d\epsilon ^{\prime} }(\epsilon ^{\prime} ,\epsilon ,{\epsilon }_{\gamma })\,\rho ({\epsilon }_{\gamma })\,{\rm{\Delta }}\epsilon \,d\epsilon ^{\prime} d{\epsilon }_{\gamma },\end{array}\end{eqnarray} \tag{ 5 }$$

where the constant is on the order of 10⁻⁴ when Δs is expressed in Mpc and all energies in the temperature-scaled units. ${\epsilon }_{\gamma ,\min }$ is a lower limit on the energy of the photon, which we set arbitrarily to 0.19/ so as to keep the kinematic limits of Equation (3) real. The upper limit is set to 0.64/ for the same reason.

Similar to Equation (5), the energy-loss dispersion rate in scaled energy units is written,

$$\begin{eqnarray}\begin{array}{rcl}{{ \mathcal K }}_{2}(\epsilon ) & = & \displaystyle \frac{\langle {\left({\rm{\Delta }}\epsilon \right)}^{2}\rangle }{{\rm{\Delta }}s}(\epsilon )=\mathrm{const}{\int }_{{\epsilon }_{\gamma ,\min }}^{{\epsilon }_{\gamma ,\max }}{\int }_{\epsilon {{\prime} -}^{}}^{\epsilon {{\prime} +}^{}}\displaystyle \frac{d\sigma }{d\epsilon ^{\prime} }\\ & & \times \ (\epsilon ^{\prime} ,\epsilon ,{\epsilon }_{\gamma })\,\rho ({\epsilon }_{\gamma })\,{\left({\rm{\Delta }}\epsilon \right)}^{2}\,d\epsilon ^{\prime} d{\epsilon }_{\gamma }.\end{array}\end{eqnarray} \tag{ 6 }$$

In scaled energy units, the Fokker–Planck transport equation in energy space can now be expressed as,

$$\begin{eqnarray}&&\displaystyle \frac{\partial f}{\partial s}=-\displaystyle \frac{\partial }{\partial \epsilon }\left[{{ \mathcal K }}_{1}(\epsilon )f(\epsilon ,s)\right]+\displaystyle \frac{{\partial }^{2}}{\partial {\epsilon }^{2}}\left[{{ \mathcal K }}_{2}(\epsilon )f(\epsilon ,s)\right],\end{eqnarray} \tag{ 7 }$$

where ds is the differential pathlength and ${{ \mathcal K }}_{1}(\epsilon )=-\displaystyle \frac{\langle {\rm{\Delta }}\epsilon \rangle }{{\rm{\Delta }}s}(\epsilon )$. Shown in Figure 1 are the differential energy-loss cross-section $d\sigma /d\epsilon ^{\prime}$ as a function of energy at three incident proton energies: = 0.5 (1890 EeV); = 0.2 (757 EeV); and = 0.08 (302 EeV). For each incident energy, shown in red is the sum of the three processes' differential energy-loss cross-section, where green depicts $d\sigma /d\epsilon ^{\prime}$ for the reaction γ_CMB + p → π⁺ + n, blue for γ_CMB + p → π⁰ + p, and brown for γ_CMB + n → π⁻ + p. Each differential cross-section was estimated using the parameterization in Equation (2) where the expansion coefficients a_n are taken from measurements (Genzel et al. 1973) up to and including the third, n = 2, term.

Zoom In Zoom Out Reset image size

Figure 2 shows the calculated (discrete points) energy-loss rate, according to Equation (5). The dashed curve is a best fit to these points;

$$\begin{eqnarray}&&| {{ \mathcal K }}_{1}(\epsilon )| \approx {c}_{1}\arctan ({c}_{2}\,{\epsilon }^{{c}_{3}}),\end{eqnarray} \tag{ 8 }$$

with c₁ = 0.0123 Mpc⁻¹, c₂ = 15.65, and c₃ = 1.325. Similarly, Figure 3 shows the calculated (discrete points) Fokker–Planck diffusion coefficient, ${{ \mathcal K }}_{2}(\epsilon )$, according to Equation (6). The dashed curve is a best fit to these points:

$$\begin{eqnarray}&&{{ \mathcal K }}_{2}(\epsilon )\approx \displaystyle \frac{{d}_{1}{\epsilon }^{2}}{{d}_{2}+{d}_{3}{\epsilon }^{2}},\end{eqnarray} \tag{ 9 }$$

with d₁ = 95.344 Mpc⁻¹, d₂ = 5272.7, and d₃ = 24714. It can be seen from Figure 3 that the diffusive energy losses become important at high energies. However, the propagation effect is a function of proton energy and propagation distance.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Since both the systematic and stochastic energy losses are purely collisional in our model, we can estimate a characteristic dispersive energy-loss distance by noting that the characteristic energy-loss scale is ∼$\sqrt{2}\,{{ \mathcal K }}_{2}(\epsilon )/{{ \mathcal K }}_{1}(\epsilon )$ (Payne 1969; Tsendin 2006). Dividing this energy-loss scale by ${{ \mathcal K }}_{1}(\epsilon )/\sqrt{2}$ should give a corresponding characteristic dispersive energy-loss distance. Using these coefficients we can estimate a characteristic dispersive distance, ℓ(), using:

$$\begin{eqnarray}&&{\ell }(\epsilon )=\displaystyle \frac{2{{ \mathcal K }}_{2}(\epsilon )}{{{ \mathcal K }}_{1}^{2}(\epsilon )}.\end{eqnarray} \tag{ 10 }$$

The calculated ℓ() is shown in Figure 4, which can be considered as a distance scale to measure the importance of diffusive energy losses, as a function of proton energy. If a proton with energy is propagating from a source located at a distance s Mpc, the characteristic dispersive attenuation distance is l(). As an example, for a proton with energy = 0.01 this characteristic distance is ≈19.5 Mpc, while it is ≈5.5 Mpc for a proton with energy = 0.1. If these two protons were to propagate 10 Mpc distance, then the effect of dispersion in their respective energy losses is more important for the proton that has l = 5.5 Mpc. Figure 4 also suggests that over a rather wide energy range, ∈ [0.05, 0.5], ℓ changes only by a factor of 2. Given the definition of ℓ this suggests that dispersion in energy loss remains equally important as the systematic energy loss over this large dynamic energy range. We expound on this feature and implications to the GZK cutoff in Section 2.3. As a result, when studying the transport of UHECR protons, their energy loss and propagation distance cannot be fully separated. In other words, the importance of dispersion in energy losses is determined by the primary proton's energy as well as propagation distance.

Zoom In Zoom Out Reset image size

To derive the Fokker–Planck potential, Equation (7) needs to be re-expressed with a constant diffusion coefficient (Risken 1989) using the transformation in the independent variables, → ξ (see the Appendix):

$$\begin{eqnarray}&&{\left(\displaystyle \frac{d\xi }{d\epsilon }\right)}^{2}\times {{ \mathcal K }}_{2}(\epsilon )={{ \mathcal K }}_{2}^{{\prime} }=\mathrm{constant},\end{eqnarray} \tag{ 11 }$$

which gives, for $d{{ \mathcal K }}_{2}^{{\prime} }/d\epsilon =0$ over a finite domain for ∈ [0.005, 0.5],

$$\begin{eqnarray}&&\xi (\epsilon )={\int }_{.005}^{\epsilon }{\left(\displaystyle \frac{{{ \mathcal K }}_{2}^{{\prime} }}{{{ \mathcal K }}_{2}(\zeta )}\right)}^{1/2}d\zeta .\end{eqnarray} \tag{ 12 }$$

The transformed drift term is given by

$$\begin{eqnarray}&&{{ \mathcal K }}_{1}^{{\prime} }[\xi (\epsilon )]={\left(\displaystyle \frac{{{ \mathcal K }}_{2}^{{\prime} }}{{{ \mathcal K }}_{2}(\epsilon )}\right)}^{1/2}\left({{ \mathcal K }}_{1}(\epsilon )-\displaystyle \frac{1}{2}\displaystyle \frac{d{{ \mathcal K }}_{2}(\epsilon )}{d\epsilon }\right).\end{eqnarray} \tag{ 13 }$$

In terms of the new independent variable, ξ, and the dependent variable, $\tilde{f}(\xi ,s)$, Equation (7) is written (note though that the time-like variable, s, remains the same),

$$\begin{eqnarray}&&\displaystyle \frac{\partial \tilde{f}(\xi ,s)}{\partial s}=-\displaystyle \frac{\partial }{\partial \xi }\left[{{ \mathcal K }}_{1}^{{\prime} }(\xi )\tilde{f}(\xi ,s)\right]+{{ \mathcal K }}_{2}^{{\prime} }\displaystyle \frac{{\partial }^{2}\tilde{f}(\xi ,s)}{\partial {\xi }^{2}}.\end{eqnarray} \tag{ 14 }$$

In formal analogy with the Hamiltonian operator (Englefield 1988; Risken 1989),

$$\begin{eqnarray}&&{ \mathcal H }=-V(\xi )+{{ \mathcal K }}_{2}^{{\prime} }\displaystyle \frac{{\partial }^{2}}{\partial {\xi }^{2}},\end{eqnarray} \tag{ 15 }$$

where ${{ \mathcal K }}_{2}^{{\prime} }$ is an arbitrary constant and V(ξ) is a single-particle Schrodinger potential, we can define a Fokker–Planck "potential" as

$$\begin{eqnarray}&&\phi (\xi )=-{\int }^{\xi }{{ \mathcal K }}_{1}^{{\prime} }[\xi ^{\prime} ]d\xi ^{\prime} ,\end{eqnarray} \tag{ 16 }$$

which is related to the Schrodinger potential by $V{(\xi )=1/(4{ \mathcal K }{{\prime} }_{2}^{2})(d\phi /d\xi )}^{2}-(1/2){d}^{2}\phi /d{\xi }^{2}$. A connection between the Fokker–Planck and Schrodinger equations can help guide solutions, because large number solutions exist for quadratic and cubic potential, as well demonstrated by Englefield (1988). In terms of the Fokker–Planck potential, Equation (14) becomes:

$$\begin{eqnarray}&&\displaystyle \frac{\partial \tilde{f}(\xi ,s)}{\partial s}=\displaystyle \frac{\partial }{\partial \xi }\left[\displaystyle \frac{\partial \phi (\xi )}{\partial \xi }\tilde{f}(\xi ,s)\right]+{{ \mathcal K }}_{2}^{{\prime} }\displaystyle \frac{{\partial }^{2}\tilde{f}(\xi ,s)}{\partial {\xi }^{2}}.\end{eqnarray} \tag{ 17 }$$

Figure 5 shows the derived Fokker–Planck potential using Equations (13) and (16) and setting ${{ \mathcal K }}_{2}^{{\prime} }=\tfrac{1}{2}$. We numerically solved for the energy points ∈ [0.005, 0.5], which correspond to ξ ∈ [0.0, 25.6] and E ∈ [1.89 × 10¹⁹, 1.89 × 10²¹] eV, where data are available, that satisfy Equations (13) and (16) simultaneously. The darker solid curve is a fourth-order polynomial fit to the calculated discrete points. The light solid curve is a third-order polynomial fit to the same points. The fourth-order fit is given by ϕ⁽⁴⁾(ξ) = − 17.309 × 10⁻⁶ ξ⁴ + 0.000399ξ³ + 0.00693ξ² + 0.251ξ + 0.0265. The third-order fit is given by ϕ⁽³⁾(ξ) =− 5.644 × 10⁻⁴ ξ³ + 0.0244ξ² + 0.142ξ + 0.142. The quartic potential has a zero at $\xi ={\xi }_{\max }=41.281$ (6.58 × 10²¹ eV), and a maximum at ξ = ξ_m = 28.67 (2.70 × 10²¹ eV). The value of the quartic potential at maximum is 10.63, in units of ξ² per Mpc, which are the same units ${{ \mathcal K }}_{2}^{{\prime} }$ is expressed as unity in. Similarly, the cubic potential has a zero at $\xi ={\xi }_{\max }=48.54$ (8.93 ×10²¹ eV), and a maximum at ξ = ξ_m = 31.49 (3.51 × 10²¹ eV). The value of the cubic potential at maximum is 11.19 in units of ξ² per Mpc. The two fits are statistically different, the error sum of squares (SSE) for the quartic fit is 0.0385 while it is 0.129 for the cubic fit. Thus, the quartic fit is more significant and is the one we use to carry out the numerical solution of Equation (14) in Section 3. However, both fits are identical in the energy range of interest and we use both of them to estimate the horizon distance and probability to clear the potential barrier.

Zoom In Zoom Out Reset image size

Next we estimate the mean distance, as a function of the proton energy it takes the UHECR proton to overcome this potential barrier, i.e., to go from the energy at which the potential is maximum, which is at ξ = ξ_m , to the equilibrium point of the potential at ξ = 0. This distance, which we call ${ \mathcal S }(\xi )$, is calculated using two different approaches. In the first approach we follow the same formalism (Gardiner 1983) typically used to calculate the mean first passage time,

$$\begin{eqnarray}&&\begin{array}{l}{ \mathcal S }(\xi )=\displaystyle \frac{1}{{{ \mathcal K }}_{2}^{{\prime} }}{\int }_{\xi }^{0}{\int }_{{\xi }_{m}}^{\xi ^{\prime} }\exp (+\phi (\xi ^{\prime} )/{{ \mathcal K }}_{2}^{{\prime} })\\ \,\times \,\exp (-\phi (\xi ^{\prime\prime} )/{{ \mathcal K }}_{2}^{{\prime} })d\xi ^{\prime\prime} d\xi ^{\prime} .\end{array}\end{eqnarray} \tag{ 18 }$$

Note that the expression for ${ \mathcal S }(\xi )$ above is a solution to the second-order ODE (Hanggi et al. 1990):

$$\begin{eqnarray}&&{{ \mathcal K }}_{2}^{\prime} \displaystyle \frac{{d}^{2}{ \mathcal S }}{d{\xi }^{2}}-\displaystyle \frac{d\phi (\xi )}{d\xi }\,\displaystyle \frac{d{ \mathcal S }(\xi )}{d\xi }=-1,\end{eqnarray} \tag{ 19 }$$

where ξ ∈ [0, ξ_m ], and where we assume that ${ \mathcal S }\to 0$ at ξ = 0, i.e., an absorbing boundary, and $d{ \mathcal S }/d\xi \to 0$ at ξ = ξ_m , a reflecting boundary.

The second approach is based on the continuous-slowing-down approximation (CSDA), where ${ \mathcal S }(\xi )$ is calculated using:

$$\begin{eqnarray}&&{{ \mathcal S }}^{\mathrm{CSDA}}(\xi )={\int }_{\xi }^{0}\displaystyle \frac{\,d\xi ^{\prime} }{{{ \mathcal K }}_{1}(\xi ^{\prime} )}.\end{eqnarray} \tag{ 20 }$$

Calculating the horizon distance using Equation (18) involves both the systematic and stochastic energy losses in a compact form (the potential form), while the CSDA calculations involve only the systematic energy losses. Two horizon distances result from Equation (18); ${{ \mathcal S }}^{(3)}(\xi )$ results from using the third-order potential fit, ϕ⁽³⁾, and similarly, ${{ \mathcal S }}^{(4)}(\xi )$ is the horizon distance that results from the fourth-order potential fit, ϕ⁽⁴⁾.

Figure 6 shows the calculated mean horizon distance, ${ \mathcal S }(\xi )$, as a function of proton energy. Several defining characteristics can be seen. First, the two potential fits seem to give very similar distance-versus-energy curves, as their behavior is essentially the same except for energies at and higher than about ξ ≈ 27 (2.2 × 10²¹ eV). Distance inferred from the cubic potential, ${{ \mathcal S }}^{(3)}$, reaches a maximum of around 88 Mpc at an energy of about ξ = 32 (3.67 × 10²¹ eV), while ${{ \mathcal S }}^{(4)}$ reaches its maximum of about 75 Mpc at ξ = 28 (2.51 × 10²¹ eV). The maximum horizon distance corresponds to an energy that maximizes the potential. Noting that the ratio of the two potentials' heights, ${\phi }_{\max }^{(3)}/{\phi }_{\max }^{(4)}$, is only about 1.05, which translates to $\sqrt{1.05}$ along ξ, while the ratio of their widths at equilibrium, ${\xi }_{m}^{(3)}/{\xi }_{m}^{(4)}$ is about 1.2, for this high energy regime the distance inferred from the potential appears more sensitive to the width of the potential than its height.

Zoom In Zoom Out Reset image size

Two energy regimes can be seen from Figure 6. For energies up to about ξ ≈ 25(1.75 × 10²¹ eV), and for both potentials S(ξ) appears to rise exponentially with ξ. For energies between ξ ≈ 25 and ξ ≈ 35 (1.75 × 10²¹ eV–4.58 × 10²¹ eV), the rise reaches a maximum then a drop in the horizon distance is seen.

Additionally, we note that the inferred distance for both potentials increases by only a factor of 3 in about 2 orders of magnitude in proton energy; from 20 Mpc at ξ = 5 (4.89 × 10¹⁹ eV) to 60 Mpc at ξ = 25 (1.75 × 10²¹ eV). This clearly suggests that the GZK cutoff is quite broad. We quantify the degree of this broadness in Section 2.4.

Figure 6 also shows the calculated horizon distance based on the CSDA approximation, which is purely determined by the systematic energy losses. A difference between the calculated horizon distances with and without the diffusive energy-loss effects are seen at high energies. As can be seen from Figures 3 and 4, stochastic energy losses are important at high energies; so at high enough energies the diffusive effect is dominant and has a clear impact on the horizon distance calculations. These are indicators of the importance of dispersion in energy losses in UHECR propagation. Likewise, we expect to see an energy-loss dispersion effect on the protons' energy spectrum, which we address in Section 3.

In addition to the horizon distance, using the derived Fokker–Planck potential we can explicitly calculate the likelihood of a proton at a certain energy (ξ > 0) reaching the lower limit of the ξ-domain, which we set to zero (see Section 2.2). This probability measure is calculated using (Gardiner 1983):

$$\begin{eqnarray}&&\displaystyle { \mathcal M }(\xi )=\frac{{\int }_{0}^{\xi }\exp (+\phi ({\xi }^{{\rm{{\prime} }}})/{{ \mathcal K }}_{2}^{{\rm{{\prime} }}})\,d{\xi }^{{\rm{{\prime} }}}}{{\int }_{0}^{{\xi }_{max}}\exp (+\phi ({\xi }^{{\rm{{\prime} }}{\rm{{\prime} }}})/{{ \mathcal K }}_{2}^{{\rm{{\prime} }}})\,d{\xi }^{{\rm{{\prime} }}{\rm{{\prime} }}}.}\end{eqnarray} \tag{ 21 }$$

Figure 7 shows the calculated probability for a proton with energy ξ to clear the Fokker–Planck potential ϕ(ξ). This figure shows that if a particle with energy less than ξ ≈ 20 (8.05 × 10²⁰ eV) cleared the potential barrier, then the probability of this particle being a proton is almost zero, and hence, it is more likely to be a heavy ion if a charged particle is observed at this energy. Conversely, at much higher energies, ξ > 35 (4.58 × 10²¹ eV), there is almost a 100% chance that these particles are protons. This estimated probability agrees with the fact that heavy ions would have experienced photodisintegration below this energy. On the other hand, if a charged particle with energy ξ ≈ 30 (3.08 × 10²¹ eV) is observed, then there is around 50% probability for this particle to be a proton and 50% probability for it to be a heavy ion. Previous analysis of the Pierre Auger Observatory measured composition indicates heavier elements at energies beyond 2 × 10¹⁸ eV. Additionally, Aab et al. (2014, 2016) concluded that this fraction of the spectrum being proton is zero above 10¹⁹ eV with an indication of higher fractions at higher energies. These findings are qualitatively consistent with Figure 7.

Zoom In Zoom Out Reset image size

This estimated probability is for protons at steady-state. In order to compare our model results with other experiments, similar calculations as a function of distance are needed. Additionally, similar calculations are needed for each heavy ion. This will enable a more precise interpretation of the observations and identify particles based on their energies. Unfortunately, these calculations are not available at present due to the dearth of knowledge of measured photodisintegration differential cross-section of many heavy ions of interest. However, Figure 7 can still be used to distinguish between protons and heavy ions as observed charged particles.