Curious case of the maximum rigidity distribution of cosmic-ray accelerators

In many models, the sources of ultrahigh-energy cosmic rays (UHECRs) are assumed to accelerate particles to the same maximum energy. Motivated by the fact that candidate astrophysical accelerators exhibit a vast diversity in terms of their relevant properties such as luminosity, Lorentz factor, and magnetic field strength, we study the compatibility of a population of sources with nonidentical maximum cosmic-ray energies with the observed energy spectrum and composition of UHECRs at Earth. For this purpose, we compute the UHECR spectrum emerging from a population of sources with a power-law, or broken-power-law, distribution of maximum energies, applicable to a broad range of astrophysical scenarios. We find that the allowed source-to-source variance of the maximum energy must be small to describe the data if a power-law distribution is considered. Even in the most extreme scenario, with a very sharp cutoff of individual source spectra and negative redshift evolution of the accelerators, the maximum energies of 90% of sources must be identical within a factor of 3—in contrast to the variance expected for astrophysical sources. Substantial variance of the maximum energy in the source population is only possible if the maximum energies follow a broken-power-law distribution with a very steep spectrum above the break. However, in this scenario, the individual source energy spectra are required to be unusually hard with increasing energy output as a function of energy.

Figure 1

Left: cosmic-ray source spectra for the different rigidity cutoff functions. $R_{\max }$ denotes the maximum rigidity and the $y$ axis is scaled to show the ratio to an unmodified power law with spectral index ${\gamma }_{\mathrm{src}}$. Right: population spectra resulting from the convolution of a power-law distribution of maximum rigidities above rigidity $R_0$ and the source spectra displayed in the left panel (${\beta }_{\mathrm{pop}}=4$).

Figure 2

Schematic description of how the parameters of the considered astrophysical scenarios ($\eta$, $\alpha$, $\xi$, $y_1$, $y_2$, ${\beta }_1$, ${\beta }_2$, and ${\gamma }_{\mathrm{src}}$) are connected to the parameters (${\beta }_{\mathrm{pop}}$ and ${\gamma }_{\mathrm{pop}}$) that describe the effective population spectrum ${\varphi }_{\mathrm{pop}}$. In the first step, the single-PL or BPL distributions of physical parameters are converted to a PL or BPL distribution of maximum rigidities. This conversion is exact, assuming that Eqs. (19) and (24) are valid. If $R_{\max }$ follows a PL distribution, with slope $R_{\max }^{-{\beta }_{\mathrm{pop}}}$ above the threshold at $R_0$, we obtain our default case for the population spectrum. If $R_{\max }$ is distributed according to a broken power law, with break at $R_T$ and slope $R_{\max }^{-{\beta }_1}/R_{\max }^{-{\beta }_2}$ before or after the break, the same parametrization of ${\varphi }_{\mathrm{pop}}$ is possible but only after a reinterpretation of the parameters. This approach is exact except around the break at $R_0$.

Figure 3

Predicted spectrum and composition at Earth for the best-fit scenario of the fiducial model [sibyll2.3c, $⟨X_{\max }⟩-{\sigma }_{\mathrm{syst}}$, $\sigma (X_{\max })+{\sigma }_{\mathrm{syst}}$]. The colored bands indicate the contributions of the separate mass groups with $[A_{\min },A_{\max }]$, including the 68% uncertainties (1 d.o.f.). Hatched areas indicate systematic uncertainties of the data. Data points at $E<{10}^{18.4}\mathrm{eV}$ (crosses) are taken from [83] and are only shown for visual guidance. Only points above ${10}^{18.8}\mathrm{eV}$ are used in the fit.

Figure 4

Results of the source parameter scan for the fiducial model marginalized along all but two axes respectively. The surface plot shows the agreement between prediction and Auger observations in terms of the ${\chi }^2$ estimator and the contour lines indicate the one (green) and three (red) sigma confidence interval for two degrees of freedom. The best fit is marked with a white cross.

Figure 5

Fit quality for the population model with broken-power-law distribution of maximum rigidities, marginalized onto ${\beta }_1\times {\beta }_2$ space (the spectral index of the population spectrum below and beyond the break). Contours indicate the one (green) and three (red) sigma confidence intervals (2 d.o.f.). The white-shaded region denotes the parameter space where $R_{\max }^{0.90}/R_0<10$, i.e. where the spread in maximum rigidity is less than a decade for 90% of sources. The best fit is marked with a white cross. Parameters of potential source classes predicted based on their luminosity functions [Eq. (24)] are shown as black points. The allowed values of ${\beta }_2$ for blazars in the $R_{\max }(\mathrm{\Gamma })$ scenario [Eq. (19)] for the Hillas-constrained case and under the assumption that ${\gamma }_{\mathrm{src}}\ge 2$ are indicated on the left side.

Figure 6

Best-fit source spectral index ${\gamma }_{\mathrm{src}}$ (top) and break rigidity $R_0$ (bottom) as a function of ${\beta }_1$ and ${\beta }_2$ for the population model with a broken-power-law distribution of maximum rigidities. Best-fit confidence contours and source candidates are indicated as in Fig. 5. The structures in the upper left of the plots are artifacts from the limited resolution of our sampling grid and some degeneracy in the choice of $(R_0,{\gamma }_{\mathrm{src}})$ for a particular combination of $({\beta }_1,{\beta }_2)$.

Figure 7

Results of the source parameter scan for the population model with redshift evolution of the distribution of maximum rigidities marginalized onto ${\beta }_{\mathrm{pop}}\times {\gamma }_{\mathrm{src}}$ space. The agreement between prediction and Auger observations in terms of the ${\chi }^2$ estimator is displayed with different color levels, and the contour lines indicate the one (green) and three (red) sigma confidence intervals for two degrees of freedom. The best fit is marked with a white cross.

Figure 8

Predicted cosmogenic neutrino flux associated with the best fit (solid line) of the population model with redshift evolution of $p(R_{\max })$. The shaded one and three sigma uncertainty bands correspond to the contours in Fig. 7. We show the observed IceCube HESE flux [67], upper limits from IceCube [68] and Auger [96], and predicted sensitivities of planned detectors [93, 94, 97, 98] as a reference.

Figure 9

Results of the source parameter scan for the population model with superexponential source cutoff function, with exponent ${\lambda }_{\mathrm{cut}}$, marginalized onto ${\beta }_{\mathrm{pop}}\times {\lambda }_{\mathrm{cut}}$ space. The surface plot shows the agreement between prediction and Auger observations in terms of the ${\chi }^2$ estimator and the contour lines indicate the one (green) and three (red) sigma confidence interval for two degrees of freedom. The best fit is marked with a white cross.

Figure 10

Power-law (solid blue) and broken-power-law (solid red) distribution of maximum rigidities $p(R_{\max })$ with break or starting point at $R_0$, and ${\gamma }_{\mathrm{src}}=2$, ${\beta }_1=-0.8$, ${\beta }_2=4.5$. The one-sided 90% quantile of the PL (blue dashed line) and two-sided 90% quantile of the BPL (red dotted lines) are also displayed. In this scenario, the PL approach underestimates the population variance by a factor of $\sim 3.8$.