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

D. Ehlert  0000-0002-4322-6400<EMAIL_ADDRESS>Institutt for fysikk,
Norwegian University of Science and Technology (NTNU), 7491 Trondheim, Norway
F. Oikonomou  0000-0002-0525-3758 Institutt for fysikk, Norwegian University
of Science and Technology (NTNU), 7491 Trondheim, Norway  M. Unger
0000-0002-7651-0272 Institute for Astroparticle Physics, Karlsruhe Institute
of Technology (KIT), 3640 Karlsruhe, Germany  Institutt for fysikk, Norwegian
University of Science and Technology (NTNU), 7491 Trondheim, Norway

###### Abstract

In many models, the sources of ultra-high-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 non-identical 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 three – 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.

##  I Introduction

Ultra-high-energy cosmic rays (UHECRs) are charged particles that reach Earth
with energies of up to several $10^{20}$ $\mathrm{eV}$. The identification of
the astrophysical sources capable of accelerating particles to these energies
is one of the unsolved mysteries of high energy astrophysics (see e.g. [1, 2]
for recent reviews). A correlation between astrophysical objects and the
measured arrival directions of cosmic rays has not yet been established at
high significance [3], but the source properties are constrained by
measurements of the diffuse particle flux and composition at Earth; see e.g.
[4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20].

Most of these studies assume an acceleration mechanism that is universal in
rigidity 111The rigidity of a particle with charge $Ze$ and momentum $p$ is
$R=p\,c/(Ze)\simeq E/Z$ (using natural units and with energy $E$). up to a
maximum rigidity of $R_{\text{max}}$, leading to consecutive flux suppressions
of the elemental spectra at energies of $E_{\text{max}}=ZR_{\text{max}}$,
where $Z$ denotes the cosmic-ray charge. Assuming such a “Peters cycle” [21,
22] at the sources gives a good description of the flux and composition
measured at Earth at ultra-high energies, see e.g. [11]. Depending on the
source environment, the maximum energy can follow a different functional form
[23]. Here, we consider primarily the canonical scenario where the maximum
energy scales in proportion to the nuclear charge and briefly discuss the
effect of alternate scalings on our results.

A major caveat of many studies is that the sources are typically assumed to be
identical – a description that is unlikely to hold for realistic sources. The
most probable astrophysical candidates for the sources of UHECRs, e.g. active
galactic nuclei (AGN) and gamma-ray bursts (GRBs), are generally not very
similar – even within a single source class – but exhibit an enormous
diversity in terms of key parameters like luminosity, size, magnetic field and
jet power.

Only a few studies have relaxed the assumption of identical sources in the
past, by focusing on a low number of discrete local sources [24, 25, 26] or by
considering the superposition of a few ($\leq$ 3) source classes, e.g. [15,
27, 18, 28, 29, 20]. The time variation of $R_{\text{max}}$ in AGN jets was
studied in Ref. [30], and the effective spectrum produced by sources with non-
identical spectral shapes and spectral indices has been discussed in the
context of gamma-ray spectra [31] and Galactic cosmic-ray sources [32].
Populations of sources with non-identical maximum cosmic-ray rigidities were
considered previously, for a pure proton UHECR composition [33], for Galactic
cosmic rays [34], and for gamma-ray bursts [16].

Here we present, for the first time, a rigorous exploration of the population
variance of $R_{\text{max}}$ compatible with current observations of the
spectrum and composition of UHECRs. This is achieved by convolving the
distribution of source properties, parameterised by the maximum rigidity, with
the individual source spectra to obtain an analytical description of the total
population spectrum, as detailed in Sec. II. We simulate the propagation of
UHECRs to Earth through the extragalactic photon fields and find the best
source parameters by comparing the model predictions to UHECR data in Sec.
III. From this, we derive lower limits on the source variance allowed by the
data, as described in Sec. IV. We conclude in Sec. V that only a limited
amount of population variance is permitted, and UHECR sources are required to
be nearly identical in terms of maximum rigidity under realistic choices of
the model parameters if a power-law distribution of maximum rigidities is
assumed. The variance can be large for broken power-law distributions,
provided the source spectra are sufficiently hard. However, we find the
distributions predicted for most candidate source classes – blazars, gamma-ray
bursts, and tidal disruption events – to be incompatible with the limits
obtained from the UHECR fit. Only sources with luminosity distribution similar
to the one of Seyfert galaxies can potentially satisfy the constraints.

##  II Population Spectrum of non-identical Sources

In this study, we assume that the rigidity spectra of individual sources of
UHECRs are well described by the aforementioned Peters Cycle, and thus we
assume a power law with a high-rigidity cutoff,

$\phi_{\text{src}}=\frac{\operatorname{d\\!}^{\,2}N}{\operatorname{d\\!}R\operatorname{d\\!}t}=\sum_{i}\phi_{0}(Z_{i})\,R^{-\gamma_{\text{src}}}\,f(R,R_{\text{max}})\,,$
(1)

where the sum runs over all accelerated chemical elements with charge $Z_{i}$
and the spectral index is assumed to be
universal.222$\gamma_{\text{src}}\approx 2$ for diffusive shock acceleration
[35, 36], but in this study the value of $\gamma_{\text{src}}$ is a free
parameter. The term $f(R,R_{\text{max}})$ describes the high-rigidity cutoff
at maximum rigidity $R_{\text{max}}$. We refer to the sum of the spectra of
all sources within a certain volume as the population spectrum
$\phi_{\text{pop}}$.

In the limit of identical sources, the population spectrum will necessarily
have the same shape as the spectra of individual sources. A source-by-source
variation of the normalisation factors $\phi_{0}(Z_{i})$ does not lead to a
qualitatively different population spectrum since it is equivalent to
identical sources with the source-averaged normalisations $\phi_{0}$, and thus
we will not consider it in the following. It should, however, be kept in mind
that the elemental fractions obtained from the fits presented later in this
paper should be understood as source-averaged fractions. A phenomenologically
more interesting source property is the maximum rigidity. If the probability
for an individual accelerator to reach a certain maximum rigidity is
distributed as
$p(R_{\text{max}})\equiv\operatorname{d\\!}p/\operatorname{d\\!}R_{\text{max}}$
and the source spectra follow $\phi_{\text{src}}(R,R_{\text{max}})$, then the
combined spectrum of the entire population is given by the convolution

$\phi_{\text{pop}}(R)=\int_{0}^{\infty}\phi_{\text{src}}(R,R_{\text{max}})\,p(R_{\text{max}})\,\operatorname{d\\!}R_{\text{max}}.$
(2)

Here and in subsequent occurrences of $\phi_{\text{src}}$, the sum over all
chemical elements of charge $Z_{i}$ (cf. Eq. 1) is assumed implicitly. In the
following, we specify the functional forms of individual source spectra and
probability distribution of $R_{\text{max}}$ that will be studied in this
work.

###  II.1 Source Spectra

In general, the rigidity cutoff of astrophysical accelerators depends on the
acceleration mechanism and the source environment, in particular on the
dominant energy-loss process, see e.g. [37, 38]. The simplest description of
the shape is given by a sharp termination when particles exceed the maximum
rigidity

$\phi_{\text{src}}^{\text{hs}}=\phi_{0}\,R^{-\gamma_{\text{src}}}\,\theta(R_{\text{max}}-R),$
(3)

where $\theta(x)$ denotes the Heaviside step function. The population spectrum
corresponding to this function has been studied previously in Ref. [33]. Since
$\phi_{\text{src}}^{\text{hs}}$ describes the sharpest-possible rigidity
cutoff, it provides a useful extreme case that will allow for a maximum
variation of maximum rigidity among the sources in the population.

A more commonly used choice, expected under certain astrophysical conditions,
see e.g. [38, 39], is given by an exponential cutoff

$\phi_{\text{src}}^{\text{exp}}=\phi_{0}\,R^{-\gamma_{\text{src}}}\exp{\left(-\frac{R}{R_{\text{max}}}\right)}\,.$
(4)

However, this function has the disadvantage that the effect of the cutoff
already starts to become noticeable well below the maximum rigidity and the
interpretation of the spectral index $\gamma_{\text{src}}$ is complicated. For
this reason, some phenomenological studies assume a broken-exponential source
spectrum, e.g. [11]

$\phi_{\text{src}}^{\text{b-exp}}=\phi_{0}\,R^{-\gamma_{\text{src}}}\,\begin{cases}1&R<R_{\text{max}}\\\
\exp\left(1-\frac{R}{R_{\text{max}}}\right)&\text{otherwise},\end{cases}$ (5)

which alleviates the issue but lacks physical motivation.

Finally, we consider spectra with an exponential cutoff raised to the power of
${\lambda_{\text{cut}}}$,

$\phi_{\text{src}}^{\text{s-exp}}=\phi_{0}\,R^{-\gamma_{\text{src}}}\,\exp\left(-\frac{R}{R_{\text{max}}}\right)^{\lambda_{\text{cut}}},\quad{\lambda_{\text{cut}}}>0\,.$
(6)

We refer to this description as a “super-exponential” cutoff. The function can
be used to interpolate continuously between classical, exponential cutoffs
(${\lambda_{\text{cut}}}=1$) and sharp, Heaviside-like terminations
(${\lambda_{\text{cut}}}=\infty$). Additionally, for
${\lambda_{\text{cut}}}<1$ the cutoff shape becomes sub-exponential up to no
cutoff when ${\lambda_{\text{cut}}}\rightarrow 0$. Super-exponential cutoff
profiles were obtained e.g. in Ref. [39] with ${\lambda_{\text{cut}}}=2$ for
synchrotron losses during acceleration.

An illustration of the source spectrum for different choices of the cutoff is
shown in the left panel of Fig. 1.

Figure 1: Left: Cosmic-ray source spectra for the different rigidity cutoff
functions. $R_{\text{max}}$ denotes the maximum rigidity and the $y$-axis is
scaled to show the ratio to an unmodified power law with spectral index
$\gamma_{\text{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_{\text{pop}}}=4$).

###  II.2 Distribution of Maximum Rigidities

##### Power Law

In this paper, we mainly consider a population of sources with a distribution
of maximum rigidities that follows a power law (PL) with spectral index
${\beta_{\text{pop}}}$ above a minimum allowed maximum rigidity $R_{0}$

$p(R_{\text{max}})=\begin{cases}0&R_{\text{max}}<R_{0}\\\
\frac{{\beta_{\text{pop}}}-1}{R_{0}}\,\left(\frac{R_{\text{max}}}{R_{0}}\right)^{-{\beta_{\text{pop}}}}&\text{otherwise,}\end{cases}$
(7)

which is also known as a Pareto distribution. This $R_{\text{max}}$
distribution was previously considered in Ref. [33]. Because of the asymmetric
nature of the power-law distribution, the standard deviation is of limited use
to characterize the source variance, and instead we will report the one-sided
$90\%$ quantile $R_{\text{max}}^{0.90}$ defined as

$\int_{R_{0}}^{R_{\text{max}}^{0.90}}\operatorname{d\\!}R_{\text{max}}\,p(R_{\text{max}})=0.90\,.$
(8)

For a power-law distribution of maximum rigidities, the quantile $q$ is given
by the relation

$R_{\text{max}}^{q}/R_{0}=[1-q]^{1/(1-{\beta_{\text{pop}}})}.$ (9)

For illustration, population diversity of more than a decade, i.e.
$R_{\text{max}}^{0.90}/R_{0}\geq 10$, is obtained if ${\beta_{\text{pop}}}<2$.

##### Broken Power Law

Alternatively, the distribution of maximum rigidities can be modelled as a
broken power law (BPL) which can be written as

$p(R_{\text{max}})=\frac{R_{0}^{-1}}{C}\times\begin{cases}\left(\frac{R_{\text{max}}}{R_{0}}\right)^{-\beta_{1}}&R_{\text{max}}\leq
R_{0}\\\
\left(\frac{R_{\text{max}}}{R_{0}}\right)^{-\beta_{2}}&R_{\text{max}}>R_{0},\end{cases}$
(10)

with $C$ a normalisation constant. A detailed discussion of the BPL scenario
is provided in Appendix B.

###  II.3 Population Spectrum

##### Power Law

Assuming a power-law distribution of $R_{\text{max}}$, it is possible to
derive an analytical description of the population spectrum for all source
spectra presented in Sec. II.1. In the case of sources with a Heaviside
termination in rigidity, the population spectrum is given by,

$\phi_{\text{pop}}^{\text{hs}}=\phi_{0}\,R^{-\gamma_{\text{src}}}\begin{cases}1&R<R_{0}\\\
\left(\frac{R}{R_{0}}\right)^{-{\beta_{\text{pop}}}+1}&\text{otherwise}.\end{cases}$
(11)

For source spectra with a broken exponential cutoff, the population spectrum
becomes,

$\displaystyle\phi_{\text{pop}}^{\text{b-exp}}=\phi_{0}\,R^{-\gamma_{\text{src}}}\begin{cases}1&R<R_{0}\\\
\left(\frac{R}{R_{0}}\right)^{-{\beta_{\text{pop}}}+1}f(\frac{R}{R_{0}},{\beta_{\text{pop}}})&\text{otherwise}\end{cases}$
(12)

with

$\displaystyle
f\left(\textstyle\frac{R}{R_{0}},{\beta_{\text{pop}}}\right)=\,$
$\displaystyle 1+e\,({\beta_{\text{pop}}}-1)$
$\displaystyle\times\left[\gamma\left({\beta_{\text{pop}}}-1,\frac{R}{R_{0}}\right)-\gamma\big{(}{\beta_{\text{pop}}}-1,1\big{)}\right].$

Finally, for sources with a (super-)exponential cutoff, the population
spectrum is given by,

$\displaystyle\phi_{\text{pop}}^{\text{s-exp}}=\phi_{0}\,R^{-\gamma_{\text{src}}}\,$
$\displaystyle\left(\textstyle\frac{R}{R_{0}}\right)^{-{\beta_{\text{pop}}}+1}\,\textstyle\frac{{\beta_{\text{pop}}}-1}{{\lambda_{\text{cut}}}}$
$\displaystyle\times\,\gamma\left(\textstyle\frac{{\beta_{\text{pop}}}-1}{{\lambda_{\text{cut}}}},\left(\frac{R}{R_{0}}\right)^{{\lambda_{\text{cut}}}}\right).$
(13)

For a standard exponential distribution, i.e. ${\lambda_{\text{cut}}}=1$ this
simplifies to

$\displaystyle\phi_{\text{pop}}^{\text{exp}}=\phi_{0}\,R^{-\gamma_{\text{src}}}$
$\displaystyle\left(\textstyle\frac{R}{R_{0}}\right)^{-{\beta_{\text{pop}}}+1}({\beta_{\text{pop}}}-1)$
$\displaystyle\times\gamma\left(\textstyle{\beta_{\text{pop}}}-1,\frac{R}{R_{0}}\right),$
(14)

whereas for ${\lambda_{\text{cut}}}\to\infty$, Eq. 11 is recovered. Here
$\gamma$ denotes the lower incomplete gamma function, not to be confused with
the source spectral index $\gamma_{\text{src}}$.

Population spectra for a particular choice of parameters are shown in Fig. 1,
right. The main impact of the cutoff parameter ${\lambda_{\text{cut}}}$ is a
later onset and faster turnover of the spectral break around $R_{0}$. The
limiting behaviour of the population spectrum is independent of the source
cutoff function. The asymptotic rigidity dependencies of the population
spectra are

$\lim\limits_{R\to 0}\phi_{\text{pop}}(R)\propto R^{-\gamma_{\text{src}}}$
(15)

and

$\lim\limits_{R\to\infty}\phi_{\text{pop}}(R)\propto
R^{-\gamma_{\text{src}}-{\beta_{\text{pop}}}+1}.$ (16)

##### Broken Power-Law

In the broken power-law scenario, the resulting population spectrum for
sources with super-exponential cutoff is

$\displaystyle\phi_{\text{pop}}^{\text{s-exp}}$
$\displaystyle=\frac{\phi_{0}}{C\cdot{\lambda_{\text{cut}}}}\,R^{-\gamma_{\text{src}}}\times\Big{[}L+H\Big{]},$
(17)

with

$\displaystyle
L=\left(\frac{R}{R_{0}}\right)^{-\beta_{1}+1}\cdot\Gamma\left(\frac{\beta_{1}-1}{{\lambda_{\text{cut}}}},\left(\frac{R}{R_{0}}\right)^{{\lambda_{\text{cut}}}}\right)$
$\displaystyle
H=\left(\frac{R}{R_{0}}\right)^{-\beta_{2}+1}\cdot\gamma\left(\frac{\beta_{2}-1}{{\lambda_{\text{cut}}}},\left(\frac{R}{R_{0}}\right)^{{\lambda_{\text{cut}}}}\right).$
(18)

Here “$C$” is a normalisation factor (see Appendix B), “$L$” reflects the
contribution from sources with $R_{\text{max}}\leq R_{0}$ and “$H$” represents
sources above the break. $\Gamma$ and $\gamma$ are the upper and lower
incomplete gamma function respectively. The single power-law scenario is
retrieved as a limiting case of the more general broken power-law case for
$\beta_{1}\to-\infty$, i.e. when the contribution “$L$” of
low-$R_{\text{max}}$ sources vanishes. We discuss the broken power-law
scenario in more detail in Appendix B. Here we only mention that the same
expression for the population spectrum as in the power-law picture holds
approximately after defining the effective spectral index
$\gamma_{\text{pop}}=\gamma_{\text{src}}+\beta_{1}-1$ and
$R_{\text{max}}$-distribution index
${\beta_{\text{pop}}}=\beta_{2}-\beta_{1}+1$.

###  II.4 Relation to Astrophysical Quantities

The population spectra derived in the previous section provide simple analytic
expressions that are well suited for fits to UHECR observations at Earth, with
which the key parameters ${\beta_{\text{pop}}}$ and $\gamma_{\text{src}}$ can
be derived. The connection of these parameters to the properties of UHECR
sources is discussed in the following for a few examples. We will show that
the assumed power-law distribution in maximum rigidity can be attributed to
different acceleration scenarios. The relevant parameters of these scenarios
are summarised in Table 1 and the re-interpretation of the fitted parameters
${\beta_{\text{pop}}},\,\gamma_{\text{src}}$ in terms of proposed underlying
physical properties are listed in Table 2. Through the re-interpretation of
the fitted parameters, all considered physical scenarios can be reduced to the
same population spectra that are obtained for a power-law distribution of
$R_{\text{max}}.$ The process is illustrated in Fig. 2.

Table 1:  Summary of parameters used for the three illustrative scenarios. See sections given in the first column for further details. Scenario | Parameter | Description | Equation
---|---|---|---
Power law | $\gamma_{\text{src}}$ | True spectral index of the sources | $\phi_{\text{src}}\propto R^{-\gamma_{\text{src}}}$
(Sec. II.2) | $\gamma_{\text{pop}}$ | Effective spectral index of the sources (same as the true spectral |
|  | index for a power-law distribution of $p(R_{\text{max}})$) | $\phi_{\text{pop}}\propto R^{-\gamma_{\text{pop}}}$
| ${\beta_{\text{pop}}}$ | (Effective) spectral index of the $p(R_{\text{max}})$ distribution | $p(R_{\text{max}})\propto R_{\text{max}}^{-{\beta_{\text{pop}}}}$
| $\beta_{1},\,\beta_{2}$ | Spectral index of the distribution of maximum rigidities below (above) |
|  | the break for a broken power law in $p(R_{\text{max}})$ | $p(R_{\text{max}})\propto R_{\text{max}}^{-\beta_{1}}\,(R_{\text{max}}^{-\beta_{2}})$
Lorentz | $\eta$ | Spectral index of the power-law distribution of Lorentz-factors | $\operatorname{d\\!}p/\operatorname{d\\!}\Gamma\propto\Gamma^{-\eta}$
factor | $\alpha$ | Energy boosting by the relativistic motion of the jet | $E=E^{\prime}\,\Gamma^{\alpha}$
(Sec. II.4.1) | $\xi$ | Time dilation caused by the relativistic motion of the source region | $t=t^{\prime}\,\Gamma^{\xi}$
Luminosity | $y_{1},\,y_{2}$ | Spectral index of the broken power-law luminosity function of sources |
(Sec. II.4.2) |  | below (above) the break | $\operatorname{d\\!}p/\operatorname{d\\!}L\propto L^{-y_{1}}\,(L^{-y_{2}})$
ptPL$\frac{\operatorname{d\\!}p}{\operatorname{d\\!}\Gamma}(\eta)$ptPL$\frac{\operatorname{d\\!}p}{\operatorname{d\\!}L}(y_{2})$ptBPL$\frac{\operatorname{d\\!}p}{\operatorname{d\\!}L}(y_{1},y_{2})$ptLorentz
Factor$\Gamma$Sec. II.4.1 Luminosity $L$ Sec. II.4.2 PL
$p(R_{\text{max}}|{\beta_{\text{pop}}})$
ptBPL$p(R_{\text{max}}|\beta_{1},\beta_{2})$pt${\beta_{\text{pop}}}=$$\frac{\eta-1}{\alpha}+2$$-\gamma_{\text{src}}+\frac{\xi}{\alpha}$pt${\beta_{\text{pop}}}=$$2y_{2}-3$pt$\beta_{1}=2y_{1}-3$$\beta_{2}=2y_{2}-3$
Population Spectrum
$\phi_{\text{pop}}({\beta_{\text{pop}}},\gamma_{\text{pop}})$
pt${\beta_{\text{pop}}}={\beta_{\text{pop}}}$$\gamma_{\text{pop}}=\gamma_{\text{src}}$pt$\bf\underline{\beta_{1}<1}$${\beta_{\text{pop}}}\approx\beta_{2}$$\gamma_{\text{pop}}\approx\gamma_{\text{src}}$$\bf\underline{\beta_{1}>1}$${\beta_{\text{pop}}}=\beta_{2}-\beta_{1}+1$$\gamma_{\text{pop}}=\gamma_{\text{src}}+\beta_{1}-1$
Figure 2: Schematic description of how the parameters of the considered
astrophysical scenarios ($\eta$, $\alpha$, $\ \xi$, $y_{1}$, $y_{2}$,
$\beta_{1}$, $\beta_{2}$, $\gamma_{\text{src}}$) are connected to the
parameters $({\beta_{\text{pop}}},\gamma_{\text{pop}})$ that describe the
effective population spectrum $\phi_{\text{pop}}$. In the first step, the
single power-law (PL) or broken power-law (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_{\text{max}}$ follows a PL distribution, with slope
$R_{\text{max}}^{-{\beta_{\text{pop}}}}$ above the threshold at $R_{0}$, we
obtain our default case for the population spectrum. If $R_{\text{max}}$ is
distributed according to a broken power law, with break at $R_{T}$ and slope
$R_{\text{max}}^{-\beta_{1}}/R_{\text{max}}^{-\beta_{2}}$ before/after the
break, the same parameterisation of $\phi_{\text{pop}}$ is possible but only
after a re-interpretation of the parameters. This approach is exact except
around the break at $R_{0}$.

####  II.4.1 Jet Lorentz Factor

In some scenarios, for sources with relativistic jets, the maximum rigidity is
directly related to the bulk Lorentz factor of the motion, $\Gamma_{\rm jet}$.
For instance, the Hillas criterion [40] for relativistic sources gives
$R_{\text{max}}=R_{0}\Gamma_{\rm jet}$, with $R_{0}\propto lB$ where $l$ is
the size of the source and $B$ the magnetic field.

It is also possible that UHECRs are galactic cosmic rays that receive a “one-
shot” boost of a factor of $\sim\Gamma_{\rm jet}^{2}$ in the jet of their host
galaxies, in which case $R_{\text{max}}\sim R_{0}\,\Gamma_{\rm jet}^{2}$,
where $R_{0}$ is the maximum energy of the cosmic rays before re-acceleration.
This is referred to as the _Espresso_ mechanism [41, 42, 43]. Cosmic rays that
do not enter the most relativistic parts of the jet are only partially boosted
with $R_{\text{max}}\propto\Gamma_{\rm jet}^{\alpha},\,\alpha<2$. Thus, here
we investigate the general case of

$R_{\text{max}}=R_{0}\,\Gamma_{\rm jet}^{\alpha}$ (19)

where the aforementioned cases are described by $\alpha=1$ (Hillas) and
$\alpha\leq 2$ (Espresso).

Assuming, for example, a power-law distribution of the Lorentz factors, as
found consistent with observations of jetted AGN in Refs. [44, 45],

$\frac{\operatorname{d\\!}p}{\operatorname{d\\!}\Gamma_{\rm
jet}}=(\eta-1)\,\Gamma_{\rm jet}^{-\eta},$ (20)

the distribution of maximum rigidities can be calculated as

$\displaystyle p(R_{\text{max}})$
$\displaystyle=\frac{\operatorname{d\\!}p}{\operatorname{d\\!}\Gamma_{\rm
jet}}\left|\mathinner{\dfrac{\operatorname{d\\!}{}\Gamma_{\rm
jet}}{\operatorname{d\\!}{R_{\text{max}}}}}\right|$ (21)
$\displaystyle=\frac{\eta-1}{\alpha}\,R_{0}^{-1}\,\left(\frac{R_{\text{max}}}{R_{0}}\right)^{\frac{1-\eta}{\alpha}-1}\theta(R_{\text{max}}-R_{0}).$

Of course, $R_{0}$ is also expected to vary from source to source, and
therefore the distribution of $R_{\text{max}}$ should be broader and the above
equation can be understood as a lower limit on the source-to-source variation
of $R_{\text{max}}$.

The boosting of particle energies also affects the expected flux emitted by
individual sources. This introduces additional terms into Eq. 21 but we show
in Appendix A that the convolution of the source spectra and the
$R_{\text{max}}$ distribution lead to the same functional forms as derived in
the last section. However, the parameter ${\beta_{\text{pop}}}$ can now be
related to physical properties of the source population, namely the spectral
index $\gamma_{\text{src}}$, the Lorentz-boosting factor
$R_{\text{max}}\propto\Gamma_{\text{jet}}^{\alpha}$, and the distribution of
Lorentz factors $p(\Gamma_{\text{jet}})\propto\Gamma_{\text{jet}}^{-\eta}$ via

${\beta_{\text{pop}}}=\frac{\eta-1}{\alpha}+2-\gamma_{\text{src}}+\xi/\alpha\,,$
(22)

where the time dilation factor $\xi=1$ for an acceleration region co-moving
with the jet and $\xi=0$ for espresso-type re-acceleration.

####  II.4.2 Luminosity

Another plausible distribution of the maximum rigidity can be derived from the
minimum luminosity requirement for particle acceleration in expanding flows.
Here, the minimum luminosity, $L_{\rm 0}$, needed to accelerate CRs to maximum
rigidity $R_{0}$ is given by [46, 47, 48, 49, 50, 51, 52],333The normalisation
value varies slightly between different papers.

$L_{\rm 0}\approx
10^{45.5}~{}\frac{1}{\beta}\left(\frac{R_{0}}{10^{20}\,$\mathrm{V}$}\right)^{2}~{}{\rm
erg~{}s^{-1}},$ (23)

where $\beta$ is the speed of the flow in units of $c$. In this scenario, we
can relate $R_{\text{max}}$ and observed luminosity $L$ of a source via

$R_{\text{max}}\sim R_{0}\,\beta^{1/2}\,\left(\frac{L}{L_{0}}\right)^{1/2}.$
(24)

The impact of the variance of $\beta$ in a population can be approximately
neglected for highly-, as well as mildly- and non-relativistic sources.444If
the acceleration region is highly relativistic, then $\beta\sim 1$ and there
is no additional variance introduced by the non-identical outflow speeds. Even
in non- or mildly-relativistic source environments, the impact is expected to
be small. For example, the authors of [53] found that the relation between
luminosity and outflow speed in a sample of AGN outflows is $L\sim
v_{\text{out}}^{4.6}$. The additional contribution in Eq. 24 is thus
$\beta^{1/2}\sim L^{1/9}$, which constitutes only a subdominant effect
compared to the original $R_{\text{max}}\propto L^{1/2}$.

In this scenario, where we relate the maximum rigidity to the source
luminosity,

$p(R_{\text{max}}|z)=\dfrac{\operatorname{d\\!}{}p}{\operatorname{d\\!}{L}}(z)\left|\mathinner{\dfrac{\operatorname{d\\!}{}L}{\operatorname{d\\!}{R_{\text{max}}}}}\right|$
(25)

where $\operatorname{d\\!}p(z)/\operatorname{d\\!}L$ is the luminosity
function of the sources. For a single power-law distribution of luminosities,
which can adequately describe many proposed source classes, and without taking
into account the redshift evolution,555For our purposes, luminosity and
density redshift evolution of sources are indistinguishable in terms of the
total contribution to the observed UHECR energy flux. However, as the maximum
rigidity is related to the source luminosity (Eq. 24), an evolution of maximum
rigidities, $L(z)\to R_{\text{max}}(z)$, would be introduced. The impact of
such an evolution is studied in Sec. IV.4, where we find that the cosmic-ray
fit is not very sensitive to this behaviour but negative evolutions, i.e.
$R_{\text{max}}\propto z^{\zeta}\,,\zeta<0$, are preferred at moderate
significance. we can write

$\frac{\operatorname{d\\!}p}{\operatorname{d\\!}L}=\frac{y_{2}-1}{L_{0}}\,\left(\frac{L}{L_{0}}\right)^{-y_{2}}\,.$
(26)

We assume that the emitted flux of a single source scales with the luminosity
as $\phi_{\text{src}}\propto L/L_{0}$. Noting Eq. 24, this introduces an
additional dependency of the source flux on the maximum rigidity, which can be
absorbed into the $p(R_{\text{max}})$ distribution by adjusting the definition
of the effective slope ${\beta_{\text{pop}}}$. The distribution of maximum
rigidities is then

$p(R_{\text{max}})=\frac{2\,(y_{2}-1)}{R_{0}}\,\left(\frac{R_{\text{max}}}{R_{0}}\right)^{-2y_{2}+3},$
(27)

which, except for an additional normalisation constant $\kappa$, reduces to
the PL expression in Eq. 7 after defining

${\beta_{\text{pop}}}=2y_{2}-3\quad\text{and}\quad\kappa=\frac{{\beta_{\text{pop}}}+1}{{\beta_{\text{pop}}}-1}\,.$
(28)

The situation is more complex if sources follow a broken power-law luminosity
function $\operatorname{d\\!}p/\operatorname{d\\!}L(y_{1},y_{2})$. By defining
$\beta_{1}=2y_{1}-3,\,\text{and}\,\beta_{2}=2y_{2}-3$, it is possible to
express the maximum rigidity distribution as a broken power law.

The expected values of ${\beta_{\text{pop}}}$ (labelled $\beta_{\rm pop,max}$)
are listed in the last column of Table 2 for a selection of possible source
candidates. These can be compared directly to the fitted values of
${\beta_{\text{pop}}}$ discussed below. For the investigated sources
$\beta_{\rm pop,max}$ is in general low, meaning that we would expect to
observe the effect of the variance of the population in the UHECR data. It
should be kept in mind that the estimates given in this section are only a
lower limit on the source variance (upper limit on ${\beta_{\text{pop}}}$) as
we focused only on the variation of a few key parameters and treated others as
a constant (e.g. $R_{0}$) and therefore the real source variance will be
larger. In addition, for sources observing a broken power-law distribution the
variance predicted from Eq. 8 with the tabled ${\beta_{\text{pop}}}$ will
underestimate the true population variance since only sources above the break
are considered if the approximation as a single power law is made (case III.2
of Table 2). Depending on the distribution of sources below the break, the
true variance can be much larger if $\beta_{1}$ is not small. The discrepancy
is smaller if the sub-break distribution is inverted ($\beta_{1}<0$) and
approaches zero for $\beta_{1}\to-\infty$, in which case the expression
reduces to a single power law (Eq. 7).

Table 2:  Effective fit parameters, ${\beta_{\text{pop}}}$, the spectral index of the maximum rigidity distribution of the UHECR source population, $\gamma_{\text{src}}$, the assumed spectral index of the UHECR spectrum of individual sources, and their interpretation in terms of source properties for various scenarios considered in this work. The scenarios are: (I) a distribution of $p(R_{\text{max}})$ that follows an ad hoc single power law (PL) or broken power law (BPL); (II) maximum rigidity that scales as $R_{\text{max}}\propto\Gamma^{\alpha}$ with $\Gamma$ the PL-distributed bulk Lorentz factor of the acceleration region (see Sec. II.4.1) with $\operatorname{d\\!}p/\operatorname{d\\!}\Gamma\propto\Gamma^{-\eta}$; and (III) $R_{\text{max}}$ as a function of source luminosity, $R_{\text{max}}\propto\sqrt{L}$, with PL or BPL distribution of $\operatorname{d\\!}p/\operatorname{d\\!}L$ (see Sec. II.4.2). For power-law distributions the parameter in brackets denotes the slope, e.g. $\operatorname{d\\!}p/\operatorname{d\\!}\Gamma(\eta)\propto\Gamma^{-\eta}$, while for broken power-law distributions the parameters give the slope before and after the break, respectively. Scenario I represents our baseline model that we use to compute the population spectra for different source spectral cutoff functions. Cases II and III can be reduced to the former after re-interpretation of the source and population parameters $(\eta,\alpha,y_{1},y_{2},\beta_{1},\beta_{2},$ and $\gamma_{\text{src}})$ in terms of the parameters $\beta_{\rm pop,max}$ and $\gamma_{\text{pop}}$. For scenario III we quote the slopes for $\operatorname{d\\!}p/\operatorname{d\\!}L$, as opposed to $\operatorname{d\\!}p/\operatorname{d\\!}\log(L)$ which is often used in the literature. This introduces a factor of $L^{-1}$, i.e. $y_{1/2}\to y_{1/2}+1$. ID | Param. | Distribution | ${\beta_{\text{pop}}}$ | $\gamma_{\text{pop}}$ | Sources | $\beta_{\rm pop,max}$
---|---|---|---|---|---|---
I.1 | $R_{\text{max}}$ | PL, $p(R_{\text{max}}|{\beta_{\text{pop}}})$ | ${\beta_{\text{pop}}}$ | $\gamma_{\text{src}}$ |  |
I.2 | $R_{\text{max}}$ | BPL, $p(R_{\text{max}}|\beta_{1},\beta_{2})$ |  |  |  |
|  | $\beta_{1}<1$ | $\approx\beta_{2}$ | $\approx\gamma_{\text{src}}$ |  |
|  | $\beta_{1}>1$ | $\beta_{2}-\beta_{1}+1$ | $\gamma_{\text{src}}+\beta_{1}-1$ |  |
II | $R_{\text{max}}\propto\Gamma^{\alpha}$ | PL, $\operatorname{d\\!}p/\operatorname{d\\!}\Gamma(\eta)$ | $(\eta-1)/\alpha+2$ | $\gamma_{\text{src}}$ | Blazars [45] 111A steeper distribution of $\eta=2.1\pm 0.4$ was found in [60] when fitting only blazars with $\Gamma=1-40$, resulting in reduced population variance of ${\beta_{\text{pop}}}=4.1\pm 0.4$ for the Hillas and ${\beta_{\text{pop}}}=2.6\pm 0.2$ for the Espresso scenario.: $\eta=1.4$$\pm{0.2}$ |
|  |  | $-\gamma_{\text{src}}+\xi/\alpha$ |  | \+ Hillas: $\quad\,\alpha=1,~{}\xi=1$ |  $3.4$$\pm 0.2$$-\gamma_{\text{src}}$
|  |  |  |  | \+ Espresso: $\alpha=2,~{}\xi=0$ |  $2.2$$\pm 0.1$$-\gamma_{\text{src}}$
III.1 | $R_{\text{max}}\propto\sqrt{L}$ | PL, $\operatorname{d\\!}p/\operatorname{d\\!}L(y_{2})$ | $2\,y_{2}-3$ | $\gamma_{\text{src}}$ | BL Lacs [54] 222Assuming the pure-luminosity-evolution (s/m)PLE model.: $\quad\,y_{2}=2.61$$\pm 0.37$ |  $2.22$$\pm 0.74$
|  |  |  |  | FSRQs [55] 222Assuming the pure-luminosity-evolution (s/m)PLE model.: $y_{2}=2.36$$\pm 0.10$ |  $1.72$$\pm 0.20$
|  |  |  |  | Blazars [55] 222Assuming the pure-luminosity-evolution (s/m)PLE model.: $\quad\,~{}y_{2}=2.32$$\pm 0.08$ |  $1.64$$\pm 0.16$
|  |  |  |  | TDEs [56, 57]: $\quad~{}y_{2}=2.30$$\pm 0.20$ |  $1.60$$\pm 0.40$
III.2 | $R_{\text{max}}\propto\sqrt{L}$ | BPL, $\operatorname{d\\!}p/\operatorname{d\\!}L(y_{1},y_{2})$ |  |  |  |
|  | $y_{1}<2$ | $\approx 2\,y_{2}-3$ | $\approx\gamma_{\text{src}}$ | GRBs [58]: $y_{1}=1.2$${}^{+0.2}_{-0.1}$$,\,y_{2}=2.4$${}^{+0.3}_{-0.6}$ |  $1.8$${}^{+0.6}_{-1.2}$
|  |  |  |  | FSRQs [55] 222Assuming the pure-luminosity-evolution (s/m)PLE model.: $y_{1}=0$$\pm 2.07$$,\,y_{2}=2.67$$\pm 0.17$ |  $2.34$$\pm 0.34$
|  |  |  |  | Blazars [55] 222Assuming the pure-luminosity-evolution (s/m)PLE model.: $y_{1}=0.49$$\pm 1.15$$,\,y_{2}=2.79$$\pm 0.19$ |  $2.58$$\pm 0.38$
|  |  |  |  | Seyferts [59]: $y_{1}=1.96$$\pm 0.04$$,\,y_{2}=3.71$$\pm 0.09$ |  $4.42$$\pm 0.18$

##  III Methods

###  III.1 UHECR Data

We use the latest publicly available data from the Pierre Auger Observatory
for comparison with our numerical simulations. These are the energy spectrum
of UHECRs from Ref. [61], and the mean and standard deviation of the maximum
depth of air showers [62, 63] that are sensitive to the composition of cosmic
rays, see e.g. [64].

###  III.2 UHECR Propagation

UHECR injection and propagation are simulated with the numerical Monte Carlo
framework CRPropa3 [65], including the production of cosmogenic neutrinos and
gamma rays. Upper limits and measurements of the latter are qualitatively
taken into account in what follows. These are the Fermi-LAT isotropic diffuse
gamma-ray background [66], the observed IceCube high-energy starting-event
neutrino flux [67], and the IceCube $90\%$ upper limits above
$5\text{\times}{10}^{6}\text{\,}\mathrm{GeV}$ [68]. The UHECR sources are
simulated in the continuous-source approximation out to maximum redshift
$z_{\text{max}}=4$. For UHECRs in the energy range that we fit, the effective
horizon is much closer at no more than $z\lesssim 1$, e.g. [69], but sources
at larger distances can have a strong impact on the predicted flux of
cosmogenic neutrinos.

All relevant interactions are taken into account during propagation [65, 11,
70]; these are (i) redshift energy loss, (ii) photo-pion production and (iii)
electron-positron pair production on the cosmic microwave background (CMB) and
the infrared background (IRB [71]), and for heavier cosmic rays also (iv)
photo-disintegration on CMB & IRB and (v) nuclear decay.

We assume that UHECRs propagate in the ballistic regime and neglect the
effects of extragalactic magnetic fields on the trajectories of UHECRs. Based
on the results of previous studies [72, 73, 74, 29, 75, 76], these propagation
effects are mainly important at low rigidities, where they can lead e.g. to an
apparent hardening of the UHECR flux at Earth, but are not expected to alter
our conclusions about the source variance of maximum rigidities.

###  III.3 Model Fit

We compare the model-predicted UHECR spectrum and composition after
propagation to observations by the Pierre Auger Observatory. For that purpose
we convert the model composition into the air shower observables – the mean
depth of the shower maximum $\langle X_{\text{max}}\rangle$ and its standard
deviation $\sigma(X_{\text{max}})$ – following Ref. [77].

The agreement between simulation and observations is evaluated with a standard
$\chi^{2}$ estimator plus additional penalty terms

$\displaystyle\chi^{2}=\sum_{E_{i}\geq
E_{\text{min}}}\left(\frac{d_{i}-m(E_{i},{\bf
p})}{\sigma_{\text{stat}}(d_{i})}\right)^{2}+\chi^{2}_{\text{UL}}+\chi^{2}_{\text{zero}}+\chi^{2}_{\text{shifts}}\,,$
(29)

and minimized adjusting the model parameters ${\bf p}$. The sum runs over all
Auger data points at energies $E_{i}$ above the threshold energy. Here,
$d_{i}$ denotes the three measured quantities, i.e. the energy spectrum,
average $X_{\text{max}}$ and standard deviation of $X_{\text{max}}$. We select
a high value of $E_{\text{min}}=10^{18.8}\,$\mathrm{eV}$$ as the minimum
fitted UHECR energy to reduce the impact of a possible low-energy cosmic-ray
component, different from the one responsible for the highest energies
(Hillas’ “component B” [78]). We have verified that the results are consistent
within uncertainties under small changes of $E_{\text{min}}$.

The smallest set of free fit parameters $\bf{p}$ are the minimum rigidity
$R_{0}$ and slope ${\beta_{\text{pop}}}$ of the single power-law distribution
of maximum rigidities, source spectral index $\gamma_{\text{src}}$, total
population emissivity $L_{0}$, and elemental injection fractions $f_{A}^{R}$
which are defined as relative flux ratios at the same rigidity 666Fractions
are sometimes defined at the same energy instead of rigidity. A transformation
between the two parametrisations can be achieved via $f_{A}^{R}=f_{A}^{E}\cdot
Z(A)^{-\gamma_{\text{src}}+1}$.. A combination of five injection elements -
1H, 4He, 14N, 28Si and 56Fe - is used as an effective approximation of mass
groups in the cosmic-ray composition.

For alternative “interaction limited” scenarios, e.g. [23], a scaling of the
maximum energy with the CR mass is typically predicted. Yet, for most stable
injection elements, except hydrogen, there is an approximately constant
relation between CR mass and charge, with $A/Z\approx 2$, and the transition
from $E_{\text{max}}\propto Z$ to $E_{\text{max}}\propto A$ would thus only
result in a re-scaling of the reference rigidity $R_{0}$ by a constant factor.
For proton injection, where $A/Z\neq 2$, the maximum energy will not follow
the same linear scaling but will be offset by a relative factor of
$(A/Z)_{\text{p}}/(A/Z)_{\text{elem}\neq\text{p}}$. However, since in our fits
this component consistently sits below the ankle, the impact on the results is
expected to be negligible.

We treat spectral data points below the threshold as one-sided $\chi^{2}$
penalty terms that only contribute to the overall goodness-of-fit if the
predicted flux exceeds the observations. This component is denoted as
$\chi^{2}_{\text{UL}}$ in Eq. 29 and is analogous to the first term but only
evaluated if the model exceeds the data.

No cosmic rays were observed in the two highest-energy bins at
$E\geq\,10^{20.2}\,$\mathrm{eV}$$ and only $90\%$ upper limits are given in
[61]. The $\chi^{2}$ penalty derived for this type of zero-event data points
follows from the asymptotic $\chi^{2}$ term assuming a Poissonian distribution
of events [79], and is estimated as

$\chi^{2}_{\text{zero}}=\sum_{i=1}^{\rm ULs}2n_{i}^{\text{model}},$ (30)

where $n_{i}^{\text{model}}$ is the number of particles predicted by the
simulation at the energy $E(\text{UL}_{i})$ after taking into account the
detector exposure [61].

Finally, we consider the systematic uncertainties in the absolute scale of the
energy, $\langle X_{\text{max}}\rangle$ and $\sigma(X_{\text{max}})$. They are
included in the fit as nuisance parameters with

$\chi^{2}_{\text{shifts}}=\sum_{k\in\\{E,\langle
X_{\text{max}}\rangle,\sigma(X_{\text{max}})\\}}\left(\frac{\delta_{k}}{\sigma_{k}}\right)^{2},$
(31)

where the energy uncertainty is assumed as $\sigma_{E}\leq 14\%$ [80], and the
shower depth uncertainties are taken directly from the data set [63, 62]. The
scale shifts can be fitted freely in the range
$\delta_{k}\in[-\sigma_{k},\sigma_{k}]$, but unless indicated otherwise, we
fix the values of the systematic shifts $\delta_{k}$ to fiducial values, as
detailed below.

##  IV Results and Discussion

###  IV.1 Fiducial Model

The observed variance of $X_{\text{max}}$ consists of two separate
contributions: (i) shower-to-shower variations, and (ii) intrinsic shower
variability. To allow for stronger source diversity, it is necessary that
$\sigma(X_{\text{max}})$ is not already dominated by the latter. Air showers
originating from light primary cosmic rays have a larger variability [77], and
intrinsic shower variations are minimised by shifting the observational data
to the heaviest composition that is allowed within systematic uncertainties –
corresponding to an adjustment of $\langle X_{\text{max}}\rangle$ by about
$-8.5\text{\,}\mathrm{g}\text{\,}{\mathrm{cm}}^{-2}$ at the ankle and
$-7\text{\,}\mathrm{g}\text{\,}{\mathrm{cm}}^{-2}$ at the highest energies
[63, 62]. In addition, we shift the observed variance of the shower maximum up
by the systematic uncertainty to allow for the largest reasonable source
variance. We adopt these shifts as our fiducial model to allow for maximum
intrinsic source diversity. An adjustment of the energy scale is also
possible, but we found our conclusions to be invariant under such a shift and
neglect it for simplicity.

We furthermore select Sibyll2.3c [81] as our default hadronic interaction
model to convert the predicted UHECR composition to mean shower depth and
variance. This is motivated by the fact that interpreting the $X_{\text{max}}$
data of Auger with Sibyll yields a heavier composition than when using Epos-
LHC [82]. As default, we assume a flat redshift evolution of the source
density and an exponential source cutoff. Alternative redshift evolutions are
explored in Sec. IV.3 and different cutoff functions in Sec. IV.5.

In line with our outlined considerations, we find that the fiducial scale
shifts allow for increased population variance compared to the un-shifted
observations, as shown in Table 3. This is true independent of the choice of
hadronic interaction model.

Figure 3: Predicted spectrum and composition at Earth for the best-fit
scenario of the fiducial model (Sibyll2.3c, $\langle
X_{\text{max}}\rangle-\sigma_{\text{syst}}$,
$\sigma(X_{\text{max}})+\sigma_{\text{syst}}$). The coloured bands indicate
the contributions of the separate mass groups with
$[A_{\text{min}},A_{\text{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 marginalised 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.

The predicted spectrum and composition at Earth for the best-fit source
parameters of the fiducial model are shown in Fig. 3. The viable range of
${\beta_{\text{pop}}}$ is sharply bounded from below (Fig. 4), approximately
as

${\beta_{\text{pop}}}\gtrsim-\gamma_{\text{src}}+4\,$ (32)

and appreciable source diversity is only possible for soft source spectra
$\gamma_{\text{src}}\gtrsim 1$. Yet, when assuming acceleration following a
Peters cycle, soft source spectra imply a significant amount of mixing between
the different mass groups, which leads to an increase in shower variance, a
prediction that is in tension to the low values of $\sigma(X_{\text{max}})$
measured by the Pierre Auger Collaboration, and the resulting agreement
between simulation and observations is poor. This problem is exacerbated for
the soft spectra expected from diffusive shock acceleration, and values above
$\gamma_{\text{src}}\approx 2$ are excluded.

Realisations of the source model where Eq. 32 is violated are characterised by
a population spectrum with an extreme UHE tail; a consequence of
$\lim_{R\to\infty}(\phi_{\text{pop}})\propto
R^{-{\beta_{\text{pop}}}-\gamma_{\text{src}}+1}$. Such extremely-UHE cosmic
rays experience strong interactions during propagation, producing a large flux
of light secondary cosmic rays with energies up to the Greisen-Zatsepin-Kuzmin
(GZK) [84, 85] limit. In combination with the remaining non-disintegrated
component of heavy primaries, the predicted flux at Earth exhibits a large
amount of mixing between the mass groups, which leads to strong intrinsic
shower variance in excess of observations and an overall bad fit to the
observed spectral shape.

The parameter space where a good fit to the measured UHECR spectrum and
composition is achieved can be divided into two different regimes; one that
runs approximately parallel to the boundary with
${\beta_{\text{pop}}}+\gamma_{\text{src}}\approx 4-6$ in the range
$\gamma_{\text{src}}\in[-1,0.5]$, and a second that is effectively degenerate
in the population variance ${\beta_{\text{pop}}}\gtrsim 5$ with
$\gamma_{\text{src}}\in[0,1]$. The former is associated with a sub-EV maximum
rigidity threshold $R_{0}$ and a heavy composition dominated by nitrogen-like
nuclei with little contribution from lighter elements. The second regime
allows for a lighter composition of up to $\sim 50\%$ protons/helium with
$R_{0}\sim$2\text{\,}\mathrm{EV}$$. Only the second region is present in the
scenario without fiducial scale shifts applied. In both regimes, sources are
effectively identical and population variance of half a decade or more, i.e.
$R_{\text{max}}^{0.90}/R_{0}>5$, and ${\beta_{\text{pop}}}\lesssim 2.4$, can
be excluded at a confidence level above $6\sigma$ 777A penalty factor $S$ that
corrects for the quality of the global best-fit point $\chi^{2}_{\text{min}}$
is included in this estimate. We adopt the form
$S^{-1}=\sqrt{\chi^{2}_{\text{min}}/N_{\text{d.o.f.}}}$ proposed in [86]. In
essence, the penalty factor reduces the rejection strength of sub-optimal fit
points if the overall best fit is poor.. With Epos-LHC, we obtain increased
source variance at the best fit $\Delta{\beta_{\text{pop}}}\approx-0.8$ but at
the cost of a significantly reduced fit quality.

A statistically significant lower limit on the population variance cannot be
established, and identical sources cannot be rejected.

Table 3:  Best-fit parameters for several variations of the source model. From left to right: the base scenario with Sibyll2.3c as air shower model and no shifts of the energy- and $X_{\text{max}}$-scales of the data; our fiducial model, Sibyll2.3c and the $\langle X_{\text{max}}\rangle$ / $\sigma(X_{\text{max}})$ data points shifted by $-1/+1\,\sigma_{\text{syst}}$ respectively; the same scale shifts but with Epos-LHC as air shower model. The injection fractions are given in descending order for p, He, N, Si and Fe. An asterisk indicates that the confidence interval extends to the edge of the scan range and the parameter is not properly constrained in that direction. Model |  Sibyll2.3c |  Sibyll2.3c | Epos-LHC
---|---|---|---
| (no shifts) | (fid. shifts) | (fid. shifts)
$R_{0}$ [EV] | $1.73^{+0.20}_{-0.18}$ | $0.57^{+1.88}_{-0.11}$ | $1.6^{+0.6}_{-0.4}$
${\beta_{\text{pop}}}$ | $29.9^{+1.7*}_{-18.1}$ | $5.2^{+26.4*}_{-0.5}$ | $4.4^{+0.5}_{-0.5}$
$\gamma_{\text{src}}$ | $-0.23^{+0.18}_{-0.26}$ | $-0.8^{+1.4}_{-0.5}$ | $0.1^{+0.4}_{-0.5}$
$f_{A}^{R}[\%]$ | $0^{+0}_{-0}$ | $0^{+36.4}_{-0}$ | $0^{+0}_{-0}$
| $58.1^{+0.4}_{-1.9}$ | $0^{+51.3}_{-0}$ | $36.9^{+7.4}_{-22.8}$
| $35.0^{+1.6}_{-0.2}$ | $93.7^{+0.5}_{-53.5}$ | $50.3^{+16.3}_{-5.4}$
| $5.7^{+0.5}_{-0.6}$ | $0.3^{+7.7}_{-0.3}$ | $11.3^{+6.6}_{-3.8}$
| $1.16^{+0.12}_{-0.11}$ | $6.0^{+0.2}_{-3.8}$ | $1.41^{+0.27}_{-0.04}$
$R_{\text{max}}^{0.90}$ $[R_{0}]$ | $1.083^{+0.155}_{-0.005}$ | $1.72^{+0.13}_{-0.64}$ | $1.97^{+0.22}_{-0.17}$
$\chi^{2}/\text{d.o.f.}$ | $45.0/26$ | $40.4/26$ | $56.3/26$

###  IV.2 Broken Power-Law Distribution of the Maximum Rigidity

We find that a broken power-law distribution in $R_{\text{max}}$ can fit the
observed UHECR flux in the regime of highly non-identical sources
$(\beta_{1}\gtrsim 1)$ if

  1. I.a

The $R_{\text{max}}$ distribution satisfies

$\beta_{2}\gtrsim\beta_{1}+3$

as can be seen in Fig. 5 and,

  2. I.b

The effective index of the rigidity spectra of the sources satisfies

$\gamma_{\text{pop}}=1.22^{+0}_{-0.04},$

as illustrated in Fig. 6.

This suggests that the effective population spectra $R^{-\gamma_{\text{pop}}}$
must be moderately soft and the break must steepen the $R_{\text{max}}$
distribution by at least $\times R_{\text{max}}^{-3}$. Assuming reasonable
source spectra, i.e. $\gamma_{\text{src}}>0\,(>1)\,[>2$], we can place an
upper limit on the steepness of the sub-break $R_{\text{max}}$ distribution of
$\beta_{1}\lesssim 2.2\,(1.2)\,[0.2]$.

If sources are in the regime of limited population variance
$(\beta_{1}\lesssim 1)$ we find that they must satisfy

  1. II.a

$\beta_{2}\gtrsim 4.5$,
and

  2. II.b

$\gamma_{\text{src}}=-0.6^{+1.2}_{-1.0},$

consistent with the results for standard-candle sources. In summary, the key
constraints that must be observed by potential source classes for $68\%$ of
realisations are

$\displaystyle\gamma_{\text{src}}$
$\displaystyle=\begin{cases}2.22^{+0}_{-0.04}&\text{for }\beta_{1}\geq 1,\\\
-0.6^{+1.2}_{-1.0}&\text{otherwise,}\end{cases}$ (33) $\displaystyle\beta_{2}$
$\displaystyle\geq\begin{cases}\beta_{1}+3&\text{for }\beta_{1}\geq 1,\\\
4.5&\text{otherwise.}\end{cases}$ (34)

Our results clearly reflect the correlation of $\gamma_{\text{src}}$ and
$\beta_{1}$ expected from Eq. 50. Large population variance (large
$\beta_{1}$) is possible if individual source spectra are hard
($\gamma_{\text{src}}\lesssim 1$). Pareto distributions ($\beta_{2}\ll 1$)
with sufficient flatness to produce substantial source variance
($\beta_{2}\lesssim 3$) are disfavoured at $4.7\sigma$. All studied
astrophysical source classes (Table 2) are located in a region of the
parameter space where “normal” (non-inverted) spectral indices are preferred
(Fig. 6) based on the values of $\beta_{1}$ and $\beta_{2}$ that we can infer
from the studied luminosity functions. However, we found the predicted maximum
rigidity distributions to be generally incompatible with the constraints of
the UHECR fit. Only for Seyfert-like galaxies is the predicted
$R_{\text{max}}$ distribution above the break approximately compatible with
the fit to the UHECR data.

The horizontal band at $\beta_{2}\approx 6$ in Fig. 6 corresponds to a
secondary family of viable solutions separate from the global minimum, with
sub-EV break rigidity and harder injection spectra.

The best-fit parameters of the BPL and fiducial PL fit are compared in Table
5. The more-general, broken power-law approach provides an improvement of
$\Delta\chi^{2}\approx 5.7$ – which corresponds to a weak preference with
respect to a single power-law description at $2\sigma$.

Figure 5: Fit quality for the population model with broken power-law
distribution of maximum rigidities, marginalised 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_{\text{max}}^{0.9}/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_{\text{max}}(\Gamma)$ scenario (Eq. 19) for the Hillas-constrained case and
under the assumption that $\gamma_{\text{src}}\geq 2$ are indicated on the
left side.

Figure 6: Best-fit source spectral index $\gamma_{\text{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 artefacts from the limited
resolution of our sampling grid and some degeneracy in the choice of
$(R_{0},\gamma_{\text{src}})$ for a particular combination of
$(\beta_{1},\beta_{2})$.

###  IV.3 Redshift Evolution of the Source Density

For simplicity, we have so far assumed the distribution of sources to be flat
in redshift. However, the most probable source classes of UHECRs do not
exhibit this behaviour but have densities evolving as a function of redshift.
A common parametrisation of the evolution is given by

$n(z)=\begin{cases}(1+z)^{m}&\text{for }m\leq 0,\\\ (1+z)^{m}&\text{for
}m>0\text{ and }z<z_{0},\\\ z_{0}^{m}&\text{for }m>0\text{ and
}z_{0}<z<z_{\text{max}},\\\ 0&\text{otherwise}\end{cases}$ (35)

which captures the general trends for the expected source classes: a power-law
increase or decrease in density up to some break point and an approximate
flattening above that. Positive redshift evolution ($m>0$) is observed e.g.
for active galactic nuclei [87] and gamma-ray bursts [88, 58], and negative
evolution for some BL Lac sub classes [89] and tidal disruption events [88,
90, 91]. Source densities following the star formation rate are approximately
reproduced for $m=3.4$ [92].

We study the effect of source density evolution on the allowed level of
population variance by evaluating the fiducial model also for redshift
evolutions of $m=-3,3,6$ with $z_{0}=1.5$ and $z_{\text{max}}=4$. Results are
shown in Table 4.

Best agreement with observations is found for predominantly local sources
($m=-3$), and a continuous decrease in fit quality is identified for stronger
density evolutions. The improved fit for small $m$ is driven primarily by a
better match of the observed composition, in particular $\langle
X_{\text{max}}\rangle$, but the difference in $\chi^{2}$ only becomes large
once $m>3$. Redshift evolutions stronger than $m=6$ could be excluded at more
than $3\sigma$ based on their cosmogenic neutrino signature by future neutrino
detectors such as GRAND-200k [93] or IceCube Gen2 [94].

We find a clear anti-correlation between source density evolution $m$ and
spectral index $\gamma_{\text{src}}$, in agreement with previous studies, e.g.
[10, 11, 69, 95]. This is caused by the, on average, larger source distance
for stronger density evolutions and consequently increased interactions during
propagation. Since interactions soften the spectrum, a harder injection
spectrum is required at the sources. The same argument applies to the
progressively heavier source composition at the best fit. For strong
evolution, the viability of the ${\beta_{\text{pop}}}$-degenerate regime of
the fiducial model is reduced and the
${\beta_{\text{pop}}}+\gamma_{\text{src}}\approx 4$ regime is preferred more
strongly but shifted to harder source spectra $\gamma_{\text{src}}<0$.
Extremely identical sources are disfavoured in this case because they would
lead to a worse description of the observed spectral shape and an
underestimation of the shower variance.

As established previously for the fiducial model, there exists an approximate
boundary of ${\beta_{\text{pop}}}>-\gamma_{\text{src}}+4$, which dictates that
larger population variance requires softer source spectra. Local source
distributions allow for softer spectra, and the source density redshift
evolution and population variance are therefore positively correlated in the
sense that smaller values of $m$ allow for smaller values of
${\beta_{\text{pop}}}$.

To summarise, negative redshift evolutions of the source density provide many
attractive benefits: (i) a quantitatively better fit to the observed UHECR
spectrum and composition, (ii) lighter required injection composition, (iii)
more natural spectral indices $\gamma_{\text{src}}>0$, and (iv) a potentially
higher, but still not very large, population variance. This makes classes of
astrophysical objects with a negative redshift evolution, such as tidal
disruption events [88, 91] and high-spectral-peak BL Lacs [89], appealing as
sources of ultra-high-energy cosmic rays.

Table 4:  Best-fit parameters for the fiducial model but with different source density redshift evolutions, where $n(z)\sim(1+z)^{m},\,z<1.5$ (see Eq. 35). The second column ($m=0$) is identical to the fiducial scenario presented in Table 3. The injection fractions are given in descending order for p, He, N, Si and Fe. An asterisk indicates that the confidence interval extends to the edge of the scan range and the parameter is not properly constrained in that direction. Redshift |  |  |  |
---|---|---|---|---
evolution $m$ | -3 | 0 | 3 | 6
$R_{0}$ [EV] | $0.80^{+1.88}_{-0.16}$ | $0.57^{+1.88}_{-0.11}$ | $0.46^{+0.05}_{-0.09}$ | $0.52^{+0.06}_{-0.05}$
${\beta_{\text{pop}}}$ | $4.4^{+23.9}_{-0.5}$ | $5.2^{+26.4*}_{-0.5}$ | $6.46^{+0.36}_{-0.34}$ | $6.46^{+0.36}_{-0.34}$
$\gamma_{\text{src}}$ | $0.2^{+0.8}_{-0.4}$ | $-0.8^{+1.4}_{-0.5}$ | $-2.0^{+0.4}_{-0.5}$ | $-2.24^{+0.35}_{-0.18}$
$f_{A}^{R}[\%]$ | $3.5^{+46.8}_{-3.5}$ | $0^{+36.4}_{-0}$ | $0^{+0.01}_{-0}$ | $0^{+0}_{-0}$
| $8.7^{+49.8}_{-8.7}$ | $0^{+51.3}_{-0}$ | $2.6^{+17.0}_{-2.6}$ | $0^{+0}_{-0}$
| $81.3^{+11.5}_{-46.7}$ | $93.7^{+0.5}_{-53.5}$ | $90.5^{+2.0}_{-16.2}$ | $38.5^{+1.8}_{-15.7}$
| $1.7^{+3.7}_{-0.8}$ | $0.3^{+7.7}_{-0.3}$ | $0^{+0.9}_{-0}$ | $53.0^{+16.2}_{-3.9}$
| $4.8^{+0.8}_{-2.8}$ | $6.0^{+0.2}_{-3.8}$ | $6.8^{+0.5}_{-1.3}$ | $8.5^{+2.1}_{-0.5}$
$R_{\text{max}}^{0.90}$ $[R_{0}]$ | $1.97^{+0.22}_{-0.88}$ | $1.72^{+0.13}_{-0.64}$ | $1.53^{+0.04}_{-0.04}$ | $1.53^{+0.04}_{-0.04}$
$\chi^{2}/\text{d.o.f.}$ | $37.3/26$ | $40.6/26$ | $42.5/26$ | $68.9/26$

###  IV.4 Redshift Evolution of the Maximum Rigidity

In addition to the interactions with ambient photon fields, cosmic rays lose
energy due to the adiabatic expansion of the Universe, with
$E_{\text{obs}}=E_{\text{inj}}/(1+z)$. For a population of sources, this will
lead to different effective maximum rigidities for sources at different
distances and result in a naturally broadened population spectrum at Earth,
even in the limit of identical sources.

We have previously assumed that the distribution of maximum rigidities
$p(R_{\text{max}})$ does not depend on distance. However, most classes of
astrophysical objects exhibit larger luminosities at higher redshifts [87, 88,
58]. If the maximum rigidity and luminosity of a cosmic-ray source are
connected, as outlined in Sec. II.4.2, then $R_{\text{max}}$ should also
evolve as a function of redshift. We study this scenario by evolving the
starting point of the $R_{\text{max}}$ distribution with redshift,

$R_{0}(z)=R_{0}\,(1+z)^{q},~{}q\in{\rm I\\!R}\,.$ (36)

In the limit of $q=0$, we obtain the default no-redshift-scaling case while
for $q=1$ adiabatic losses are exactly compensated, and sources would have the
same effective maximum rigidity at all redshifts. Overcompensation ($q>1$) and
even enhancement of local sources ($q<0$) are also possible.

Figure 7: Results of the source parameter scan for the population model with
redshift evolution of the distribution of maximum rigidities marginalised onto
${\beta_{\text{pop}}}\times\gamma_{\text{src}}$ space. The agreement between
prediction and Auger observations in terms of the $\chi^{2}$ estimator is
displayed with different colour 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.

We find that the cosmic-ray fit has only moderate sensitivity to the value of
$q$, and no appreciable correlation with $R_{0}$, ${\beta_{\text{pop}}}$ or
$\gamma_{\text{src}}$ is observed. Nevertheless, negative evolutions are
preferred, with the best fit at $q=-4.3^{+1.0}_{-0.8*}$, and positive values
of $q\geq 1$ excluded at $3\sigma$ confidence level. The difference is
explained by a marginally better fit of the mean shower depth and a better
description of the observed spectral shape, which is related to the stacking
of contributions from different redshift shells.

Intuitively, the largest possible population rigidity variance should be
allowed for a redshift scaling of
$R_{\text{max}}(z)=R_{\text{max}}^{z=0}\cdot(1+z)$ as this would compensate
the intrinsic broadening of the maximum rigidity termination via adiabatic
energy losses. This expectation is not reflected in the results (Fig. 7), and
we find a lower limit of ${\beta_{\text{pop}}}^{\text{LL}}\approx 4$ almost
independent of $q$. A weak positive correlation of
${\beta_{\text{pop}}}^{\text{LL}}$ and $q$ can be observed, especially for the
$3\sigma$ contour. This trend is related to the reduced fragmentation of heavy
cosmic rays during propagation if the highest energies are reached only at the
most local sources. As a consequence, negative evolutions of
$R_{\text{max}}(z)$ allow for more extrinsic mixing of the mass groups due to
non-identical sources.

Perhaps the most tantalising result of this scan is not the precise best fit
but rather the realisation that strongly positive rigidity-redshift scalings
lead to an important multimessenger signature in the form of a large flux of
cosmogenic high-energy and ultra-high-energy neutrinos (Fig. 8). If the
evolution is strong, then the more distant, high-$R_{\text{max}}$ sources are
screened from our view in cosmic rays because of the interactions experienced
during propagation. At the same time, these interactions produce a large flux
of cosmogenic neutrinos that can reach us even from high redshifts. Around the
peak at $E_{\nu}\approx 10^{17.4}\,$\mathrm{eV}$$ more than $95\%$ of the
predicted neutrino flux is produced by UHECRs from sources at redshift $z>1$.
Based on this prediction, existing UHE limits of IceCube [68] and Auger [96]
are able to constrain the redshift evolution of maximum rigidities to
$q\lesssim 2$. However, we stress that this high neutrino flux is only
obtained as the most extreme scenario within $3\sigma$ and only for a low
number of specific realisations. With the increased sensitivity of future
detectors, such as IceCube Gen2 and GRAND-200k, this upper limit can be
reduced to $q\lesssim 1$ assuming a non-detection of UHE neutrinos.

The expected neutrino spectra differ substantially from predictions that are
obtained without redshift scaling of the maximum rigidity both in shape and
magnitude, especially for non-identical sources and for strongly positive
evolutions of $R_{\text{max}}$. The $99.7\%$ upper limit (2 d.o.f.) in Fig. 8,
with a peak in the neutrino flux between
$0.1\text{\,}\mathrm{EeV}1\text{\,}\mathrm{EeV}$ of up to
$E^{2}\,\operatorname{d\\!}N/\operatorname{d\\!}E\gtrsim
10^{-8}\,$\mathrm{GeV}\text{\,}{\mathrm{cm}}^{-2}\text{\,}{\mathrm{sr}}^{-1}\text{\,}{\mathrm{s}}^{-1}$$,
is obtained for approximately
$R_{0}=$1.8\text{\,}\mathrm{EV}$,\,{\beta_{\text{pop}}}=14.5,\,\gamma_{\text{src}}=0.3,\,q=2$,
and injection fractions $f_{A}^{R}\mathrel{\mathop{\mathchar
58\relax}}0,0.020,0.875,0.055,0.050$ (1H,4He,14N,28Si,56Fe). The large value
of ${\beta_{\text{pop}}}$ means that population variance is not an important
ingredient for the large predicted neutrino flux in this scenario. The
positive redshift evolution of $R_{\text{max}}$ decreases the quality of the
UHECR fit. For the same source parameters but no evolution ($q=0$), and
matching injection fractions, the fit improves by $\Delta\chi^{2}=-7.9$.

We have verified that the predicted neutrino flux is in qualitative agreement
with Refs. [69, 95] if no redshift evolution of $R_{\text{max}}$ and the same
maximum source distance are assumed. Our upper limit of the cosmogenic
neutrino flux is similar in magnitude to the flux predicted from interactions
in the source environment in Ref. [17] but offset in energy by a factor of
$\sim 10$.

Also visible in Fig. 8 is a new “shoulder” feature at ultra-high neutrino
energies beyond the classical CMB-induced peak. This feature is produced by
cosmic rays with trans-GZK energies from the highest-energy tail of the
population spectrum and is also present in the other scenarios we have
studied. The extent and magnitude are linked to the strength of the UHECR tail
and consequently to the amount of population variance.

The neutrino spectra in Fig. 8 are derived under the assumption of a source
density that does not evolve as a function of distance. For sources that are
more abundant in the local Universe the predicted flux is reduced while it is
enhanced for most other realistic source classes (AGN [87], GRB [88, 58]).

Figure 8: Predicted cosmogenic neutrino flux associated with the best fit
(solid line) of the population model with redshift evolution of
$p(R_{\text{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 [94, 93, 94, 97, 98] as a reference.

###  IV.5 Other Variations of the Source Model

We have considered additional variations of the source model to study the
impact on the allowed level of population diversity in maximum rigidity (Table
5). They are described briefly in the following. In all scenarios, a small but
generally non-zero level of population variance on the order of
${\beta_{\text{pop}}}\sim 4-6$ is recovered at the best fit to Auger
observations. The largest amount of source diversity – the most extreme case –
is obtained for negative redshift evolution of the source density and
Heaviside cutoff of the source spectra, which yields ${\beta_{\text{pop}}}\sim
3-3.5$ depending on the choice of hadronic interaction model.

##### Minimum Source Distance

Cosmic rays are attenuated during propagation depending on their source
distance. UHECRs with energies around the cutoff and with the observed heavy
composition are expected to reach us only from relatively local sources. To
avoid artifacts in the simulations we have chosen $z_{\text{min}}=10^{-3}$, or
about $4.3\text{\,}\mathrm{Mpc}$, as minimum source distance. Setting a larger
minimum distance, i.e. $43\text{\,}\mathrm{Mpc}$, we observe a decrease in fit
quality for our fiducial model. This is because nearby sources contribute
primarily at the highest energies. If these are removed, sources from larger
distances must compensate for the loss in flux. However, for a Peters cycle
progression of maximum rigidities with a preference for low maximum rigidities
$\mathcal{O}($1\text{\,}\mathrm{EeV}$)$ this compensation must come
predominantly from heavier elements since only they can reach the required
energies. Because of increased interaction due to the larger source distance,
this also leads to the production of a substantial flux of lighter secondary
cosmic rays and a stronger mixing of the mass groups – in contrast with
observations. The tension can be partially mitigated when sources are
essentially identical and the high-rigidity tail is very small. This shift to
larger values of ${\beta_{\text{pop}}}$ for the best fit is observed in our
simulations; however, the viable range is not affected strongly, and the lower
limit remains approximately the same.

##### Source Cutoff Function

Sources with an exponential UHE cutoff already include an intrinsic dispersion
in the maximum rigidity of the produced cosmic rays even for a single source.
This contribution is reduced for sharper-than-exponential cutoffs and becomes
zero for a sudden, Heaviside-like limit. An increased level of population
variance should therefore be expected for sources with a steeper cutoff
function.

As proposed in Eqs. 6 and II.3, a super-exponential cutoff can be assumed with
adjustable exponent ${\lambda_{\text{cut}}}$. Exactly Heaviside-like sources
are obtained only for ${\lambda_{\text{cut}}}\to\infty$ but effectively the
population spectra become very similar already for
${\lambda_{\text{cut}}}\gtrsim 2$. Beyond that point, the difference manifests
mainly in the increasing sharpness of the break at $R_{0}$.

Re-simulating the fiducial model for a range of exponents,
${\lambda_{\text{cut}}}\in[1,50]$, we find the global best fit at
${\lambda_{\text{cut}}}=5.4^{+1.7}_{-2.3}$, which indicates effectively
Heaviside-like sources. Overall the fit exhibits only moderate sensitivity to
the precise shape of the cutoff except when it is close to an ordinary
exponential (Fig. 9). The latter is disfavoured at a level of $2.3\sigma$. As
expected from previous results, the population variance is poorly constrained
for sources with exponential cutoff function, but even for close-to-Heaviside
source cutoffs an upper limit cannot be placed – only at
${\lambda_{\text{cut}}}\gtrsim 8$ does a significant preference of
intermediate diversity ${\beta_{\text{pop}}}\sim 4$ emerge. This reveals
substantial dependence of the population variance on the precise shape of the
spectrum at the break. We note that for approximately Heaviside-like sources
the best-fit shifts to softer source spectra and hard spectra of
$\gamma_{\text{src}}\lesssim 1$ are disfavoured. This is in better agreement
with the expectation from diffusive shock acceleration. The correlation
between the sharpness of the cutoff and rigidity variance is not as strong as
expected, and even sources with effectively instantaneous limit do not allow
for diversities much greater than ${\beta_{\text{pop}}}\sim 3.5$.

Figure 9: Results of the source parameter scan for the population model with
super-exponential source cutoff function, with exponent
${\lambda_{\text{cut}}}$, marginalised onto
${\beta_{\text{pop}}}\times{\lambda_{\text{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.

##### Fixed Injection Composition

We investigate the scenario with injection fractions reflecting the
composition observed for Galactic cosmic rays (GCRs). It is peculiar to note
that an adequate fit of the UHECR data can be achieved with (near-)GCR-like
composition if fractions are specified at the same energy [10, 27, 99],
however, this is not the case if fractions at the same rigidity are used. In
the latter case, the abundance of heavy elements is insufficient to power the
UHECR flux at the highest energies, or in other words, the predicted
composition at UHE is much lighter than indicated by observations.

We argue that rigidity fractions are the correct choice since the most
relevant interactions for cosmic rays, such as electromagnetic re-
acceleration, magnetic horizon effects, and approximately also disintegration,
are universal in rigidity and thus preserve the rigidity fractions but not the
energy fractions. The fit is possible with fractions at the same energy since
then the effective composition can be modified by changing the spectral index,
due to $f_{A}^{R}=f_{A}^{E}\,Z^{-\gamma_{\text{src}}+1}$.

We select the GCR composition at a fixed rigidity
($f_{\text{gal}}^{\text{R}}\mathrel{\mathop{\mathchar
58\relax}}\,0.794,0.190,0.013,0.002,0.001$) [100] as our baseline scenario but
allow for a free-floating rescaling of the metallicity, i.e. the abundance of
elements heavier than helium, by a factor of $\mu\in[1,100]$. With this model,
the best agreement with observations is obtained for hard source spectra
$\gamma_{\text{src}}\approx 0.7$ and population variance of
${\beta_{\text{pop}}}\approx 4.9$, which is similar to the source variance
allowed in our fiducial scenario. At the best fit, a rescaling of the
injection metallicity by a factor of approx. 60 is required, although values
as low as $\mu\approx 20$ are allowed within the $3\sigma$ confidence
interval.

Table 5:  Best-fit source spectral index $\gamma_{\text{src}}$ and maximum rigidity variance $\beta_{\text{pop}}$ plus corresponding $\chi^{2}$ for different variations of the source model. These are (fd) the fiducial model, (bp) broken power-law distribution of $p(R_{\text{max}})$, (zr) redshift evolution of $p(R_{\text{max}})$, (zn) redshift evolution of the source density, (zm) larger minimum source distance, (sc) super-exponential source cutoff function, (fg) relative injection fractions similar to the composition observed for Galactic cosmic rays, and (ex) the extreme scenario that yields the largest amount of population variance with negative redshift evolution of the source density ($m=-3$) and Heaviside rigidity cutoff at the source. The best-fit values of the additional free parameters are, $\beta_{1}=5.7^{+0.8}_{-0.4}\,\&\,\beta_{2}=23.1^{+8.5*}_{-10.8}$ for (bp), $q=-4.3^{+1.0}_{-0.8}$ for (zr), $\lambda=5.4^{+1.7}_{-2.3}$ for (sc), and $\mu=64.0^{+10.3}_{-8.8}$ for (fg). Confidence intervals that reach a limit of the scan range are marked with an asterisk. Model | Parameter | $\beta_{\text{pop}}$ | $\gamma_{\text{src}}$ | $\chi^{2}$
---|---|---|---|---
fd |  | $5.2^{+26.4*}_{-0.5}$ | $-0.8^{+1.4}_{-0.5}$ | $40.4$
bp | $\beta_{1},\beta_{2}$ | $18.4^{+8.5}_{-11.2}$ | $-3.5^{+0.2}_{-0.8}$ | $34.7$
zr | $q\in[-5,2]$ | $4.8^{+26.9*}_{-0.5}$ | $-0.19^{+0.89}_{-0.18}$ | $33.7$
zn | $m=-3$ | $4.4^{+23.9}_{-0.5}$ | $0.2^{+0.8}_{-0.4}$ | $37.3$
| $m=3$ | $6.46^{+0.36}_{-0.34}$ | $-2.0^{+0.4}_{-0.5*}$ | $42.5$
| $m=6$ | $6.46^{+0.36}_{-0.34}$ | $-2.24^{+0.35}_{-0.18}$ | $68.9$
zm | $z_{\text{min}}=0.01$ | $29.9^{+1.7*}_{-25.5}$ | $0.38^{+0.18}_{-1.22}$ | $46.2$
sc | $\lambda\in[1,50]$ | $4.0^{+3.2}_{-0.4}$ | $1.43^{+0.16}_{-0.16}$ | $33.6$
fg | $f_{A}^{R}$ | $4.9^{+0.5}_{-0.5}$ | $0.73^{+0.16}_{-0.16}$ | $45.5$
ex | Epos-LHC | $3.17^{+0.18}_{-0.17}$ | $1.43^{+0.09}_{-0.09}$ | $40.6$
|  Sibyll2.3c | $3.5^{+0.6}_{-0.5}$ | $1.69^{+0.09}_{-0.09}$ | $34.7$

##  V Conclusion

We have performed the first systematic investigation of the allowed population
variance in maximum UHECR rigidity. To this end, we have derived analytical
expressions for the population spectrum of an ensemble of non-identical UHECR
sources, assuming a (broken-)power-law distribution of maximum rigidities
$p(R_{\text{max}})\propto R_{\text{max}}^{-{\beta_{\text{pop}}}}$ and
different choices of the spectral high-energy cutoff at the sources. For the
first time, we have integrated this approach into a fit of the energy spectrum
and composition data to quantify the constraints on source similarity from
existing observations by the Pierre Auger Observatory.

If maximum rigidities are distributed according to a power law with a sudden
start (Pareto distribution) – which appears as a suitable choice for several
source candidates (e.g. AGN, blazars, TDEs) based on the observed luminosity
functions or under the assumption of power-law distributed Lorentz factors –
our results show that sources are required to be effectively identical if only
Auger data at the nominal energy and composition scale are considered. After
adjusting the measured mean shower depth and variance within systematic
uncertainties to the favourable directions that result in the heaviest
composition interpretation, we find that large yet finite values of
${\beta_{\text{pop}}}\sim 5$ are preferred, corresponding to a dispersion in
maximum rigidity of sources by a factor of approximately two.

Increased levels of population variance up to ${\beta_{\text{pop}}}\sim 3-4$
are possible for sources with sharp UHE cutoff and for source densities
evolving negatively with redshift. Even then, maximum rigidities do not differ
between sources by more than a factor of a few. In contrast, if sources are
more abundant at larger redshifts they are required to be more identical
because the preferred source spectrum becomes harder with redshift due to
increased interactions during propagation. Since the population spectrum
behaves as $\lim\limits_{R\to\infty}(\phi_{\text{pop}})\propto
R^{-\gamma_{\text{src}}-{\beta_{\text{pop}}}+1}$ a smaller source diversity
(larger ${\beta_{\text{pop}}}$) is required to limit the strength of the
spectral UHE tail.

For some source classes (e.g. GRBs, blazars, Seyferts), the luminosity
function motivates a broken power-law distribution of maximum rigidities. In
this scenario, the population variance can be large, driven by sources below
the break rigidity $R_{0}$, provided the break is sharp and the spectral index
of individual sources is sufficiently hard to counteract the variance
introduced by the non-identical sources. This requires hard spectral indices
of $\gamma_{\text{src}}\lesssim 1.2$ if the rigidity distribution below the
break is softer than $\beta_{1}\sim 1$. In addition, for any value of
$\beta_{1}$, the $R_{\text{max}}$ distribution must steepen at the break by at
least $R_{\text{max}}^{-3}$. For $\beta_{1}\to-\infty$ we obtain the power-law
scenario as an asymptotic limit, with $\gamma_{\text{src}}\lesssim 1$ and
$\beta_{2}\gtrsim 4$.

We have derived the UHECR population spectra of plausible astrophysical source
classes by connecting luminosity and maximum rigidity via the
Lovelace–Blandford–Waxman relation [46, 47, 48, 49, 50, 51], Eq. 24. For all
proposed source classes, the preferred spectral index corresponding to their
respective distribution $p(R_{\text{max}}|\beta_{1},\beta_{2})$ is in the
physically plausible range of $\gamma_{\text{src}}=0-2$ if they produce UHECRs
with a maximum rigidity distribution which follows the one we have derived
using their observed luminosity functions and Eq. 24. The predicted UHECR
population spectra produced by blazars, tidal disruption events, and GRBs are
inconsistent with the UHECR data within our formalism. Seyfert-like galaxies
are the only investigated population (see Table 2) with sufficiently steep
post-break slope to explain the required small variance at ultra-high
energies. However, a hard spectral index of $\gamma_{\text{src}}\approx 0-0.5$
is necessary. The variance in maximum rigidity obtained using Eq. 24
represents a lower limit. Additional variance of the maximum rigidity is
expected as a consequence of the distribution of other relevant source
properties that we have not considered here, which will likely reduce this
compatibility.

In summary, we have found that the maximum rigidity distribution of UHECR
sources is remarkably narrow, necessitating nearly identical (“standard-
candle”) sources, or a sharp cutoff in the rigidity distribution of the UHECR
source population. In the latter case, the low-rigidity tail exacerbates the
need for hard injection spectra that has been exposed by prior studies which
performed a combined fit of UHECR observations. Our results place strong
constraints on the most plausible astrophysical source classes of UHECRs.

Alternatively, it is possible that exotic mechanisms limit the maximum
rigidities of accelerators to the same value, e.g. [101], or that the observed
flux of UHECRs is dominated by a single local source. Such a single- or few-
source scenario seems however incompatible with the observed level of
anisotropy of the cosmic-ray arrival directions at UHE unless deflections of
cosmic rays in the Galactic and extragalactic magnetic fields are much larger
than commonly expected. An analysis of the effect of cosmic variance is beyond
the scope of this paper, but we comment that the fitted scenarios result in
typical maximum rigidities that correspond to cosmic-ray energies below the
onset of photo-nuclear interactions with the cosmic microwave background
radiation. Thus the energy-loss lengths of nuclei are large, and the volume of
UHECR sources contributing to the flux at Earth can be
$\mathcal{O}(\text{Gpc}^{3})$.

###### Acknowledgements.

We thank the anonymous referee for the constructive report. We would like to
thank Jonathan Biteau, Damiano Caprioli, Björn Eichmann, Glennys Farrar,
Michael Kachelrieß, Ioana Mariş and Marco Muzio for useful feedback on this
study. This work was made possible by Institut Pascal at Université Paris-
Saclay during the Paris-Saclay Astroparticle Symposium 2021, with the support
of the P2IO Laboratory of Excellence (program “Investissements d’avenir”
ANR-11-IDEX-0003-01 Paris-Saclay and ANR-10-LABX-0038), the P2I research
departments of the Paris-Saclay university, as well as IJCLab, CEA, IPhT,
APPEC, the IN2P3 master project UCMN and EuCAPT. We also wish to thank Sergio
Petrera for providing an update of the parameters of the $X_{\text{max}}$
parameterisation [77] for Sibyll2.3c.

##  Appendix A Transformation of the Emitted Flux

The apparent brightness of a highly relativistic source depends on the angle
$\theta$ between the observer and the direction of motion and is affected by
relativistic beaming (headlight effect) and the relativistic Doppler effect.
However, for charged particles, we can assume that their direction of motion
is isotropised after emission from the source, rendering geometrical beaming
irrelevant. Only the Lorentz boost from the rest frame of the acceleration
region to the observer frame must be considered. For sources with a flux
cutoff $f_{\text{cut}}$, the differential flux within the jet frame is given
by

$\phi^{\prime}(R^{\prime})=\frac{\operatorname{d\\!}^{\,2}N^{\prime}}{\operatorname{d\\!}R^{\prime}\operatorname{d\\!}t^{\prime}}=\phi_{0}^{\prime}\,R^{\prime-\gamma_{\text{src}}}\,f_{\text{cut}}(R^{\prime}/R_{0}).$
(37)

To transform this flux into the observer frame the following transformations
need to be taken into account

$\displaystyle
R=R^{\prime}\,\Gamma^{\alpha}~{}\text{and}~{}t=t^{\prime}\,\Gamma^{\xi}.$ (38)

Conservation of particle number implies $N=N^{\prime}$. The first
transformation is the boost in energy for particles; with $\alpha=1$ for
particles accelerated in the jet frame and emitted isotropically 888For a
general Lorentz boost of $\Gamma\,(1-\beta\cos\theta)$ with isotropic emission
angle $\theta$ the high-energy tail of the rigidity spectrum retains its
spectral shape, and we concentrate on the simpler case of a boost by
$\Gamma$., and $\alpha=2$ for the espresso mechanism. The second
transformation is due to time dilation, which also depends on the acceleration
process. If production of UHECRs within the jet is considered, then the
relativistic motion of the source region will stretch the observed time by a
factor of $\Gamma$ ($\xi=1$). On the other hand, if an espresso-like mechanism
is assumed, where the jet merely re-accelerates a pre-existing flux of cosmic
rays $\operatorname{d\\!}N/\operatorname{d\\!}t$, then no dilation is
expected, assuming that the rate of particles entering and exiting the jet is
the same, i.e.
$\operatorname{d\\!}N/\operatorname{d\\!}t|_{\text{out}}=\operatorname{d\\!}N/\operatorname{d\\!}t|_{\text{in}}$,
and thus $\xi=0$. The observed flux can then be written as

$\displaystyle\phi(R)$
$\displaystyle=\frac{\operatorname{d\\!}^{\,2}N}{\operatorname{d\\!}R\operatorname{d\\!}t}=\frac{\operatorname{d\\!}^{\,2}N(R^{\prime}(R))}{\operatorname{d\\!}R^{\prime}\operatorname{d\\!}t^{\prime}}\,\left|\frac{\operatorname{d\\!}R^{\prime}}{\operatorname{d\\!}R}\right|\,\left|\frac{\operatorname{d\\!}t^{\prime}}{\operatorname{d\\!}t}\right|$
$\displaystyle=\phi_{0}^{\prime}\,R^{-\gamma_{\text{src}}}\,\Gamma^{\alpha(\gamma_{\text{src}}-1)-\xi}\,f_{\text{cut}}{\textstyle\left(-\frac{R}{\Gamma^{\alpha}R_{0}}\right)}$
$\displaystyle=\phi_{0}^{\prime}\,R^{-\gamma_{\text{src}}}\,\left(\textstyle\frac{R_{\text{max}}}{R_{0}}\right)^{\gamma_{\text{src}}-1-\xi/\alpha}\,f_{\text{cut}}{\textstyle\left(-\frac{R}{R_{\text{max}}}\right)},$
(39)

where in the last step Eq. 19 was used. To evaluate the convolution of source
spectra and $R_{\text{max}}$ distributions, the product
$\phi(R,R_{\text{max}})\times p(R_{\text{max}})$ needs to be evaluated using
the $R_{\text{max}}$ distribution from Eq. 21. The resulting product can be
re-written in the “usual” form used in Sec. II.3,

$\displaystyle\phi(R,R_{\text{max}})$ $\displaystyle\times p(R_{\text{max}})=$
$\displaystyle\phi_{0}\,R^{-\gamma_{\text{src}}}\,f_{\text{cut}}{\textstyle\left(-\frac{R}{R_{\text{max}}}\right)}\,\textstyle\frac{{\beta_{\text{pop}}}-1}{R_{0}}\,\left(\textstyle\frac{R_{\text{max}}}{R_{0}}\right)^{-{\beta_{\text{pop}}}},$
(40)

with definitions

$\phi_{0}=\phi_{0}^{\prime}\,\frac{\eta-1}{\eta+\alpha\,(1-\gamma_{\text{src}})-\xi-1}$
(41)

and

${\beta_{\text{pop}}}=\frac{\eta-1}{\alpha}+2-\gamma_{\text{src}}+\xi/\alpha.$
(42)

Therefore, the same analytical forms derived in Sec. II.3 can be used in this
case, but the interpretation of ${\beta_{\text{pop}}}$ is more complex as it
depends on several source properties.

##  Appendix B Broken Power-Law Maximum Rigidity Distribution

Source properties, e.g. luminosity (Sec. II.4.2), often follow broken power-
law (BPL) distributions for some likely UHECR source classes rather than
single power laws (PL). A general broken power-law distribution of maximum
rigidities can be written as

$p(R_{\text{max}})=\frac{R_{0}^{-1}}{C}\cdot\begin{cases}\left(\frac{R_{\text{max}}}{R_{0}}\right)^{-\beta_{1}}&R_{\text{max}}\leq
R_{0}\\\
\left(\frac{R_{\text{max}}}{R_{0}}\right)^{-\beta_{2}}&R_{\text{max}}>R_{0},\end{cases}$
(43)

with a break at $R_{0}$ and slope $R_{\text{max}}^{-\beta_{1}}$
($R_{\text{max}}^{-\beta_{2}}$) below (above). Normalisability imposes
$\beta_{1}<1$ and $\beta_{2}>1$. Under a physically more plausible scenario,
with some minimum and maximum $R_{\text{max}}$
($R_{\text{max}}^{\text{LL}}/R_{\text{max}}^{\text{UL}}$) for the population
of sources, these conditions can be relaxed and the normalisation constant can
be expressed as

$C=\frac{1-\left(\frac{R_{\text{max}}^{\text{LL}}}{R_{0}}\right)^{-\beta_{1}+1}}{1-\beta_{1}}+\frac{1-\left(\frac{R_{\text{max}}^{\text{UL}}}{R_{0}}\right)^{-\beta_{2}+1}}{\beta_{2}-1}\,.$
(44)

If $\beta_{1}<1$ and $\beta_{2}>1$, for $R_{\text{max}}^{\text{LL}}\to 0$ and
$R_{\text{max}}^{\text{UL}}\to\infty$ this simplifies to

$C=\frac{1}{1-\beta_{1}}+\frac{1}{\beta_{2}-1}.$ (45)

Assuming source spectra with super-exponential cutoff (Eq. 6), the associated
population spectrum reads

$\displaystyle\phi_{\text{pop}}$
$\displaystyle=\frac{\phi_{0}}{C\cdot{\lambda_{\text{cut}}}}\,R^{-\gamma_{\text{src}}}\cdot\Big{[}L+H\Big{]}$
(46) $\displaystyle
L=\left(\frac{R}{R_{0}}\right)^{-\beta_{1}+1}\cdot\Gamma\left(\frac{\beta_{1}-1}{{\lambda_{\text{cut}}}},\left(\frac{R}{R_{0}}\right)^{{\lambda_{\text{cut}}}}\right)$
$\displaystyle
H=\left(\frac{R}{R_{0}}\right)^{-\beta_{2}+1}\cdot\gamma\left(\frac{\beta_{2}-1}{{\lambda_{\text{cut}}}},\left(\frac{R}{R_{0}}\right)^{{\lambda_{\text{cut}}}}\right),$

with “$L$” the contribution from sources with $R_{\text{max}}$ below the break
at $R_{0}$ and “$H$” for sources above. For $\beta_{1}\to-\infty$ a single
power law is recovered. Because of the convergence properties of the
incomplete gamma functions ($\Gamma,\gamma$),

$\lim\limits_{R/R_{0}\to\infty}L=0\quad\text{and}\quad\lim\limits_{R/R_{0}\to
0}H=0\,;\\\ $ (47)

i.e. in the limits of very small and very large rigidities, one of the terms
dominates. Only in the vicinity of the break are their contributions of
similar magnitude.

Eq. 46 suggests that the slope of the population spectrum at any rigidity is
completely specified by a set of three parameters, with the most obvious
choice ($\gamma_{\text{src}},\beta_{1},\beta_{2}$). However, an approximate
description with only two parameters is possible, similar to the power-law
scenario. The simplification is exact in the limit of very small and very
large rigidities but not in the transition region where $L$ and $H$ have
competitive levels. Two different cases can be identified depending on the
value of $\beta_{1}$. The behaviour in the high-rigidity limit is unaffected
by the distribution below the break and is always $\phi_{\text{pop}}\propto
R^{-\gamma_{\text{src}}-\beta_{2}+1}$.

  * •

$\beta_{1}\leq 1$ : Because the probability $p(R_{\text{max}})$ increases only
slowly for $R_{\text{max}}\to 0$ it can be shown that $\lim\limits_{R\to
0}\phi_{\text{pop}}\propto R^{-\gamma_{\text{src}}}$, suggesting that after
interpretation ${\beta_{\text{pop}}}=\beta_{2}$ we approximately retrieve the
same population spectrum as for the single power law where
$\beta_{1}=-\infty$. The slope of the population spectrum can then be
adequately described by the set of
$(\gamma_{\text{src}},{\beta_{\text{pop}}}=\beta_{2})$.

  * •

$\beta_{1}>1$ : The contribution from sources below the break does not vanish
and we obtain $\lim\limits_{R\to 0}\phi_{\text{pop}}\propto
R^{-\gamma_{\text{src}}-\beta_{1}+1}$. Initially, this does not appear
compatible with a PL-like description, however, re-interpretation in terms of
effective population spectral index
$\gamma_{\text{pop}}=\gamma_{\text{src}}+\beta_{1}-1$ and effective
$R_{\text{max}}$-distribution ${\beta_{\text{pop}}}=\beta_{2}-\beta_{1}+1$
leads to the usual behaviour

$\displaystyle\lim\limits_{R\to 0}$ $\displaystyle\phi_{\text{pop}}\propto
R^{-\gamma_{\text{pop}}}$ (48) $\displaystyle\lim\limits_{R\to\infty}$
$\displaystyle\phi_{\text{pop}}\propto
R^{-\gamma_{\text{pop}}-\beta_{\text{pop}}+1}\,,$ (49)

with the difference that the distribution is now characterised by the re-
defined parameter set $(\gamma_{\text{pop}},{\beta_{\text{pop}}})$.

The above prescriptions can be used to interpret the fit results for a PL
distribution given in Sec. IV for a population of sources with
BPL($R_{\text{max}}$). Yet, if not $\beta_{1}\to-\infty$, the population
variance predicted from Eq. 8 will underestimate the true variance since
sources below the break are neglected (Fig. 10). All source classes proposed
in Table 2 III.2 belong to the $\beta_{1}<1$ category, enabling the simplified
treatment.

Figure 10: Power-law (solid blue) and broken power-law (solid red)
distribution of maximum rigidities $p(R_{\text{max}})$ with break/starting-
point at $R_{0}$, and
$\gamma_{\text{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$.

Although we do not predict steep sub-break distributions for any evaluated
astrophysical source classes, it is still useful to study the scenario where
$\beta_{1}>1$. The relevant parameters are now
$(\gamma_{\text{pop}},{\beta_{\text{pop}}})$, defined as

$\displaystyle(a)~{}~{}$
$\displaystyle\gamma_{\text{pop}}=\gamma_{\text{src}}+\beta_{1}-1$ (50)
$\displaystyle(b)~{}~{}$
$\displaystyle{\beta_{\text{pop}}}=\beta_{2}-\beta_{1}+1\,.$ (51)

Solutions of this system are degenerate and different choices of
$(\gamma_{\text{src}},\beta_{1},\beta_{2})$ can lead to the same population
spectrum 999Except for the cutoff shape at $R_{0}$. Combinations of steep
source distributions and hard source spectra will lead to sharper cutoffs than
more identical sources with softer spectra.. One of the parameters must be
fixed with additional information to break the degeneracy. This is a
remarkable finding – a similar population spectrum can be achieved for highly
identical sources with softer spectral index and for diverse sources with
harder spectral index. Therefore, the distribution below the cutoff can be
steep and sources very diverse, provided individual sources have ultra-hard
spectra so that Eq. 50 remains fulfilled. If this condition is satisfied and
the break at $R_{0}$ is sufficiently strong, then the UHECR flux is dominated
by sources very close to the break at $R_{0}$, explaining the apparently low
population variance. Sources further below the break still contribute to the
flux at lower rigidities by modifying the population spectrum as
$\phi_{\text{pop}}\propto R^{-\gamma_{\text{src}}}\to
R^{-\gamma_{\text{src}}-\beta_{1}+1}$ if $\beta_{1}>1$. For flatter sub-break
distributions, their contribution is approximately negligible except close to
the break. This suggests that the sources of UHECRs can appear very similar
even if the underlying population is diverse, provided the UHE termination of
the $R_{\text{max}}$-distribution is sufficiently sharp, and the distribution
before the cutoff is not too steep after taking into account the hardness of
individual source-spectra.

## References

  * Anchordoqui [2019] L. A. Anchordoqui, Phys. Rept. 801, 1 (2019), arXiv:1807.09645 [astro-ph.HE] .
  * Alves Batista  _et al._ [2019a] R. Alves Batista _et al._ , Front. Astron. Space Sci. 6, 23 (2019a), arXiv:1903.06714 [astro-ph.HE] .
  * di Matteo  _et al._ [2021] A. di Matteo _et al._ (Telescope Array, Pierre Auger), PoS ICRC2021, 308 (2021), arXiv:2111.12366 [astro-ph.HE] .
  * Allard  _et al._ [2007a] D. Allard, E. Parizot, and A. V. Olinto, Astropart. Phys. 27, 61 (2007a), arXiv:astro-ph/0512345 [astro-ph] .
  * Hooper  _et al._ [2007] D. Hooper, S. Sarkar, and A. M. Taylor, Astropart. Phys. 27, 199 (2007), arXiv:astro-ph/0608085 .
  * Allard  _et al._ [2007b] D. Allard, A. V. Olinto, and E. Parizot, Astron. Astrophys. 473, 59 (2007b), arXiv:astro-ph/0703633 [ASTRO-PH] .
  * Allard  _et al._ [2008] D. Allard, N. G. Busca, G. Decerprit, A. V. Olinto, and E. Parizot, JCAP 0810, 033, arXiv:0805.4779 [astro-ph] .
  * Globus  _et al._ [2015] N. Globus, D. Allard, R. Mochkovitch, and E. Parizot, Mon. Not. Roy. Astron. Soc. 451, 751 (2015), arXiv:1409.1271 [astro-ph.HE] .
  * N. Globus, D. Allard, and E. Parizot [2015] N. Globus, D. Allard, and E. Parizot, Phys. Rev. D 92, 021302 (2015), arXiv:1505.01377 [astro-ph.HE] .
  * Unger  _et al._ [2015] M. Unger, G. R. Farrar, and L. A. Anchordoqui, Phys. Rev. D 92, 123001 (2015), arXiv:1505.02153 [astro-ph.HE] .
  * Aab  _et al._ [2017] A. Aab _et al._ (Pierre Auger), JCAP 04, 038, [Erratum: JCAP 03, E02 (2018)], arXiv:1612.07155 [astro-ph.HE] .
  * K. Fang and K. Murase [2018] K. Fang and K. Murase, Phys. Lett. 14, 396 (2018), [Nature Phys.14,no.4,396(2018)], arXiv:1704.00015 [astro-ph.HE] .
  * Kachelrieß  _et al._ [2017] M. Kachelrieß, O. Kalashev, S. Ostapchenko, and D. V. Semikoz, Phys. Rev. D 96, 083006 (2017), arXiv:1704.06893 [astro-ph.HE] .
  * Boncioli  _et al._ [2019] D. Boncioli, D. Biehl, and W. Winter, Astrophys. J. 872, 110 (2019), arXiv:1808.07481 [astro-ph.HE] .
  * Muzio  _et al._ [2019] M. S. Muzio, M. Unger, and G. R. Farrar, Phys. Rev. D 100, 103008 (2019), arXiv:1906.06233 [astro-ph.HE] .
  * Heinze  _et al._ [2020] J. Heinze, D. Biehl, A. Fedynitch, D. Boncioli, A. Rudolph, and W. Winter, Mon. Not. Roy. Astron. Soc. 498, 5990 (2020), arXiv:2006.14301 [astro-ph.HE] .
  * Muzio  _et al._ [2022] M. S. Muzio, G. R. Farrar, and M. Unger, Phys. Rev. D 105, 023022 (2022), arXiv:2108.05512 [astro-ph.HE] .
  * Abreu  _et al._ [2021] P. Abreu _et al._ (Pierre Auger), PoS ICRC2021, 311 (2021).
  * Bergman [2021] D. Bergman (Telescope Array), PoS ICRC2021, 338 (2021).
  * Halim  _et al._ [2022] A. A. Halim _et al._ (Pierre Auger), (2022), arXiv:2211.02857 [astro-ph.HE] .
  * Peters [1961] B. Peters, Il Nuovo Cimento 22, 800 (1961).
  * Gaisser  _et al._ [2016] T. K. Gaisser, R. Engel, and E. Resconi,  _Cosmic Rays and Particle Physics: 2nd Edition_ (Cambridge University Press, 2016).
  * Ptitsyna and Troitsky [2010] K. V. Ptitsyna and S. V. Troitsky, Phys. Usp. 53, 691 (2010), arXiv:0808.0367 [astro-ph] .
  * Ahlers  _et al._ [2013] M. Ahlers, L. A. Anchordoqui, and A. M. Taylor, Phys. Rev. D 87, 023004 (2013), arXiv:1209.5427 [astro-ph.HE] .
  * Eichmann  _et al._ [2018] B. Eichmann, J. P. Rachen, L. Merten, A. van Vliet, and J. Becker Tjus, JCAP 02, 036, arXiv:1701.06792 [astro-ph.HE] .
  * Eichmann  _et al._ [2022] B. Eichmann, M. Kachelrieß, and F. Oikonomou, JCAP 07 (07), 006, arXiv:2202.11942 [astro-ph.HE] .
  * Rodrigues  _et al._ [2021] X. Rodrigues, J. Heinze, A. Palladino, A. van Vliet, and W. Winter, Phys. Rev. Lett. 126, 191101 (2021), arXiv:2003.08392 [astro-ph.HE] .
  * Das  _et al._ [2021] S. Das, S. Razzaque, and N. Gupta, Eur. Phys. J. C 81, 59 (2021), arXiv:2004.07621 [astro-ph.HE] .
  * Mollerach and Roulet [2020] S. Mollerach and E. Roulet, Phys. Rev. D 101, 103024 (2020), arXiv:2004.04253 [astro-ph.HE] .
  * Matthews and Taylor [2021] J. H. Matthews and A. M. Taylor, Mon. Not. Roy. Astron. Soc. 503, 5948 (2021), arXiv:2103.06900 [astro-ph.HE] .
  * Lipari [2021] P. Lipari, Astropart. Phys. 125, 102507 (2021), arXiv:2001.00982 [astro-ph.HE] .
  * Yuan  _et al._ [2011] Q. Yuan, B. Zhang, and X.-J. Bi, Phys. Rev. D 84, 043002 (2011), arXiv:1104.3357 [astro-ph.HE] .
  * Kachelriess and Semikoz [2006] M. Kachelriess and D. V. Semikoz, Phys. Lett. B 634, 143 (2006), arXiv:astro-ph/0510188 .
  * Shibata  _et al._ [2010] M. Shibata, Y. Katayose, J. Huang, and D. Chen, Astrophys. J. 716, 1076 (2010).
  * Bell [1978a] A. R. Bell, Mon. Not. Roy. Astron. Soc. 182, 147 (1978a).
  * Bell [1978b] A. R. Bell, Mon. Not. Roy. Astron. Soc. 182, 443 (1978b).
  * Protheroe and Stanev [1999] R. J. Protheroe and T. Stanev, Astropart. Phys. 10, 185 (1999), arXiv:astro-ph/9808129 .
  * Protheroe [2004] R. J. Protheroe, Astropart. Phys. 21, 415 (2004), arXiv:astro-ph/0401523 .
  * Zirakashvili and Aharonian [2007] V. N. Zirakashvili and F. Aharonian, Astron. Astrophys. 465, 695 (2007), arXiv:astro-ph/0612717 .
  * Hillas [1984] A. M. Hillas, Ann. Rev. Astron. Astrophys. 22, 425 (1984).
  * Caprioli [2015] D. Caprioli, Astrophys. J. Lett. 811, L38 (2015), arXiv:1505.06739 [astro-ph.HE] .
  * Mbarek and Caprioli [2019] R. Mbarek and D. Caprioli, Astrophys. J. 886, 8 (2019), arXiv:1904.02720 [astro-ph.HE] .
  * Mbarek and Caprioli [2021] R. Mbarek and D. Caprioli, Astrophys. J. 921, 85 (2021), arXiv:2105.05262 [astro-ph.HE] .
  * Lister and Marscher [1997] M. L. Lister and A. P. Marscher, The Astrophysical Journal 476, 572 (1997).
  * Lister  _et al._ [2019] M. L. Lister _et al._ , Astrophys. J. 874, 43 (2019), arXiv:1902.09591 [astro-ph.GA] .
  * Lovelace [1976] R. V. E. Lovelace, Nature 262, 649–652 (1976).
  * Waxman [1995] E. Waxman, Phys. Rev. Lett. 75, 386 (1995), arXiv:astro-ph/9505082 .
  * Waxman [2001] E. Waxman, ICTP Lect. Notes Ser. 4, 309 (2001), arXiv:astro-ph/0103186 .
  * Blandford [2000] R. D. Blandford, Physica Scripta T85, 191 (2000).
  * Lemoine and Waxman [2009] M. Lemoine and E. Waxman, JCAP 11, 009, arXiv:0907.1354 [astro-ph.HE] .
  * Kachelriess [2022] M. Kachelriess, PoS ICRC2021, 018 (2022), arXiv:2201.04535 [astro-ph.HE] .
  * Rieger [2022] F. M. Rieger, Universe 8, 607 (2022), arXiv:2211.12202 [astro-ph.HE] .
  * Santoro  _et al._ [2020] F. Santoro, C. Tadhunter, D. Baron, R. Morganti, and J. Holt, Astronomy & Astrophysics 644, A54 (2020).
  * Ajello  _et al._ [2009] M. Ajello _et al._ , Astrophys. J. 699, 603 (2009), arXiv:0905.0472 [astro-ph.CO] .
  * Marcotulli  _et al._ [2022] L. Marcotulli _et al._ (BASS), Astrophys. J. 940, 77 (2022), arXiv:2209.09929 [astro-ph.HE] .
  * van Velzen [2018] S. van Velzen, Astrophys. J. 852, 72 (2018), arXiv:1707.03458 [astro-ph.HE] .
  * Lin  _et al._ [2022] Z. Lin, N. Jiang, X. Kong, S. Huang, Z. Lin, J. Zhu, and Y. Wang, Astrophys. J. Lett. 939, L33 (2022), arXiv:2210.14950 [astro-ph.HE] .
  * Wanderman and Piran [2010] D. Wanderman and T. Piran, Mon. Not. Roy. Astron. Soc. 406, 1944 (2010), arXiv:0912.0709 [astro-ph.HE] .
  * Ueda  _et al._ [2014] Y. Ueda, M. Akiyama, G. Hasinger, T. Miyaji, and M. G. Watson, Astrophys. J. 786, 104 (2014), arXiv:1402.1836 [astro-ph.CO] .
  * Saikia  _et al._ [2016] P. Saikia, E. Körding, and H. Falcke, Mon. Not. Roy. Astron. Soc. 461, 297 (2016), arXiv:1606.06147 [astro-ph.HE] .
  * Aab  _et al._ [2020] A. Aab _et al._ (Pierre Auger), Phys. Rev. D 102, 062005 (2020), arXiv:2008.06486 [astro-ph.HE] .
  * Aab  _et al._ [2014] A. Aab _et al._ (Pierre Auger), Phys. Rev. D 90, 122005 (2014), arXiv:1409.4809 [astro-ph.HE] .
  * Yushkov [2020] A. Yushkov (Auger), PoS ICRC2019, 482 (2020).
  * Kampert and Unger [2012] K.-H. Kampert and M. Unger, Astropart. Phys. 35, 660 (2012), arXiv:1201.0018 [astro-ph.HE] .
  * Alves Batista  _et al._ [2016] R. Alves Batista _et al._ , JCAP 05, 038, arXiv:1603.07142 [astro-ph.IM] .
  * Ackermann  _et al._ [2015] M. Ackermann _et al._ (Fermi-LAT), Astrophys. J. 799, 86 (2015), arXiv:1410.3696 [astro-ph.HE] .
  * Kopper [2018] C. Kopper (IceCube), PoS ICRC2017, 981 (2018).
  * Aartsen  _et al._ [2018] M. G. Aartsen _et al._ (IceCube), Phys. Rev. D 98, 062003 (2018), arXiv:1807.01820 [astro-ph.HE] .
  * Alves Batista  _et al._ [2019b] R. Alves Batista, R. M. de Almeida, B. Lago, and K. Kotera, JCAP 01, 002, arXiv:1806.10879 [astro-ph.HE] .
  * Alves Batista  _et al._ [2015] R. Alves Batista, D. Boncioli, A. di Matteo, A. van Vliet, and D. Walz, JCAP 10, 063, arXiv:1508.01824 [astro-ph.HE] .
  * Gilmore  _et al._ [2012] R. C. Gilmore, R. S. Somerville, J. R. Primack, and A. Domínguez, Mon. Not. Roy. Astron. Soc. 422, 3189 (2012).
  * Aloisio and Berezinsky [2004] R. Aloisio and V. Berezinsky, Astrophys. J. 612, 900 (2004), arXiv:astro-ph/0403095 .
  * González  _et al._ [2021] J. M. González, S. Mollerach, and E. Roulet, Phys. Rev. D 104, 063005 (2021), arXiv:2105.08138 [astro-ph.HE] .
  * Mollerach and Roulet [2013] S. Mollerach and E. Roulet, JCAP 10, 013, arXiv:1305.6519 [astro-ph.HE] .
  * Globus  _et al._ [2008] N. Globus, D. Allard, and E. Parizot, Astron. Astrophys. 479, 97 (2008), arXiv:0709.1541 [astro-ph] .
  * Wittkowski [2017] D. Wittkowski, PoS ICRC2017, 563 (2017).
  * Abreu  _et al._ [2013] P. Abreu _et al._ (Pierre Auger), JCAP 02, 026, arXiv:1301.6637 [astro-ph.HE] .
  * Hillas [2005] A. M. Hillas, J. Phys. G 31, R95 (2005).
  * Baker and Cousins [1984] S. Baker and R. D. Cousins, Nucl. Instrum. Meth. 221, 437 (1984).
  * Dawson [2020] B. Dawson (Pierre Auger), PoS ICRC2019, 231 (2020).
  * Riehn  _et al._ [2018] F. Riehn, H. P. Dembinski, R. Engel, A. Fedynitch, T. K. Gaisser, and T. Stanev, PoS ICRC2017, 301 (2018), arXiv:1709.07227 [hep-ph] .
  * Pierog  _et al._ [2015] T. Pierog, I. Karpenko, J. M. Katzy, E. Yatsenko, and K. Werner, Phys. Rev. C 92, 034906 (2015), arXiv:1306.0121 [hep-ph] .
  * Verzi [2019] V. Verzi, PoS ICRC2019, 450 (2019).
  * Greisen [1966] K. Greisen, Phys. Rev. Lett. 16, 748 (1966).
  * Zatsepin and Kuzmin [1966] G. T. Zatsepin and V. A. Kuzmin, JETP Lett. 4, 78 (1966).
  * Rosenfeld [1975] A. H. Rosenfeld, Ann. Rev. Nucl. Part. Sci. 25, 555 (1975).
  * Hasinger  _et al._ [2005] G. Hasinger, T. Miyaji, and M. Schmidt, Astron. Astrophys. 441, 417 (2005), arXiv:astro-ph/0506118 .
  * Sun  _et al._ [2015] H. Sun, B. Zhang, and Z. Li, Astrophys. J. 812, 33 (2015), arXiv:1509.01592 [astro-ph.HE] .
  * Ajello  _et al._ [2014] M. Ajello _et al._ , Astrophys. J. 780, 73 (2014), arXiv:1310.0006 [astro-ph.CO] .
  * Guépin  _et al._ [2018] C. Guépin, K. Kotera, E. Barausse, K. Fang, and K. Murase, Astron. Astrophys. 616, A179 (2018), [Erratum: Astron.Astrophys. 636, C3 (2020)], arXiv:1711.11274 [astro-ph.HE] .
  * Kochanek [2016] C. S. Kochanek, Mon. Not. Roy. Astron. Soc. 461, 371 (2016), arXiv:1601.06787 [astro-ph.HE] .
  * Hopkins and Beacom [2006] A. M. Hopkins and J. F. Beacom, Astrophys. J. 651, 142 (2006), arXiv:astro-ph/0601463 .
  * Álvarez-Muñiz  _et al._ [2020] J. Álvarez-Muñiz _et al._ (GRAND), Sci. China Phys. Mech. Astron. 63, 219501 (2020), arXiv:1810.09994 [astro-ph.HE] .
  * Aartsen  _et al._ [2019] M. G. Aartsen _et al._ (IceCube), Neutrino astronomy with the next generation IceCube Neutrino Observatory (2019), arXiv:1911.02561 [astro-ph.HE] .
  * Heinze  _et al._ [2019] J. Heinze, A. Fedynitch, D. Boncioli, and W. Winter, Astrophys. J. 873, 88 (2019), arXiv:1901.03338 [astro-ph.HE] .
  * Pedreira [2021] F. Pedreira (Pierre Auger), PoS ICRC2019, 979 (2021).
  * Allison  _et al._ [2016] P. Allison _et al._ (ARA), Phys. Rev. D 93, 082003 (2016), arXiv:1507.08991 [astro-ph.HE] .
  * Cummings  _et al._ [2021] A. L. Cummings, R. Aloisio, and J. F. Krizmanic, Phys. Rev. D 103, 043017 (2021), arXiv:2011.09869 [astro-ph.HE] .
  * Kimura  _et al._ [2018] S. S. Kimura, K. Murase, and B. T. Zhang, Phys. Rev. D 97, 023026 (2018), arXiv:1705.05027 [astro-ph.HE] .
  * Zyla  _et al._ [2020] P. A. Zyla _et al._ (Particle Data Group), PTEP 2020, 083C01 (2020).
  * Anchordoqui [2022] L. A. Anchordoqui, Phys. Rev. D 106, 116022 (2022), arXiv:2205.13931 [hep-ph] .