# Prospects for joint cosmic ray and neutrino constraints
on the evolution of trans-GZK proton sources

Marco Stein Muzio  0000-0003-4615-5529<EMAIL_ADDRESS>Department of Physics,
Pennsylvania State University, University Park, PA 16802, USA  Department of
Astronomy and Astrophysics, Pennsylvania State University, University Park, PA
16802, USA  Institute of Gravitation and the Cosmos, Center for Multimessenger
Astrophysics, Pennsylvania State University, University Park, PA 16802, USA
Michael Unger  0000-0002-7651-0272<EMAIL_ADDRESS>Institute for
Astroparticle Physics, Karlsruhe Institute of Technology (KIT), Karlsruhe,
Germany  Institutt for fysikk, Norwegian University of Science and Technology
(NTNU), Trondheim, Norway  Stephanie Wissel  0000-0003-0569-6978
<EMAIL_ADDRESS>Department of Physics, Pennsylvania State University,
University Park, PA 16802, USA  Department of Astronomy and Astrophysics,
Pennsylvania State University, University Park, PA 16802, USA  Institute of
Gravitation and the Cosmos, Center for Multimessenger Astrophysics,
Pennsylvania State University, University Park, PA 16802, USA

###### Abstract

We consider the prospects for future ultrahigh energy cosmic ray and neutrino
observations to constrain the evolution of sources producing a proton flux
above $10$ EeV (1 EeV$=10^{18}$ eV). We find that strong constraints on the
source evolution can be obtained by combining measurements of the cosmic ray
proton fraction above $30$ EeV with measurement of the neutrino flux at $1$
EeV, if neutrinos are predominantly of cosmogenic origin. In the case that
interactions in the source environment produce a significant astrophysical
neutrino flux, constraints on the source evolution may require measurement of
the observed proton fraction, as well as, the neutrino flux at multiple
energies, such as $1$ EeV and $10$ EeV. Finally, we show that fits to current
UHECR data favor models which result in a $>30$ EeV proton fraction and $1$
EeV neutrino flux that could realistically be discovered by the next
generation of experiments.

##  I Introduction

Over the past decade the Pierre Auger Observatory (Auger) has significantly
rewritten our understanding of the cosmic ray (CR) spectrum at ultrahigh
energies (UHEs). In particular, precise measurements of air shower properties
have led to the conclusion that UHECRs are not predominantly protons, but that
the fraction of heavier nuclei increases with energy above $10^{18.3}$ eV [1,
2, 3]. However, there is still observational and phenomenological motivation
for a flux of protons in the spectrum at the highest energies. Analysis of the
distribution of depths of shower maximum, $X_{\mathrm{max}}$, show that the
proton fraction above $10^{19}$ eV could be as high as $10\%$ in some energy
bins [4, 5]. Furthermore, a combined analysis of cosmic-ray composition and
flux results in a non-zero proton fraction above $10^{19.5}$ eV [6].
Phenomenological studies have also shown that a subdominant proton component
peaking above $10^{19}$ eV can significantly improve the fit to UHECR spectrum
and composition data [7].

Previously it was suggested (e.g. [8]) that a measurement of cosmogenic
neutrinos, i.e. neutrinos that are produced during the extragalactic
propagation of protons in interactions with cosmic photon fields, can be used
to determine the cosmic ray proton fraction. But, as pointed out in Ref. [9],
the cosmological evolution of the sources introduces a strong degeneracy that
cannot be resolved by measurements of the neutrino flux alone. On the other
hand, this implies that multimessenger studies of UHE neutrinos and cosmic
rays provide a unique opportunity to determine the evolution of sources.
Moreover, since each candidate source class exhibits a unique redshift
evolution (see Section II), constraints on the source evolution will provide
valuable insights on the thus far elusive source of UHECRs.

In this paper we consider the prospects for using both cosmic ray and neutrino
measurements to constrain the evolution of a population of UHE proton sources,
while taking into account constraints imposed by UHECR spectrum and
composition, neutrino, and gamma-ray data. We show that such constraints are
possible even when significant source interactions are considered.

##  II Model

We adopt the phenomenological Unger-Farrar-Anchordoqui (UFA) CR source model
[10], as elaborated in [7, 11]. The UFA model accounts for UHECR interactions
with photons and gas in the environment surrounding the accelerator to explain
the observed UHECR spectrum and composition, without assuming a particular
astrophysical source type. Instead, this model uses general parameters to
characterize the source’s environment, such as the average number of
interactions before escape and the temperature of the ambient photon field.
For this study we consider the superposition of two UFA-like source
populations: 1) a baseline population which accounts for the majority of the
observed UHECR spectrum and composition; and 2) a population which accelerates
a pure-proton spectrum to energies $\gtrsim 10$ EeV. To minimize the number of
free parameters we assume both populations follow the same source evolution,
but a more detailed study could be done to explore the effect of a
superposition of CR source populations with distinct evolutions. Additionally,
in order to set conservative neutrino constraints, we assume the spectral
shape of the ambient photon field to be well-characterized by a black-body
spectrum for both populations 111Other parametrizations of the source photon
field, including a broken power-law, were explored but did not have a
significant effect on our results, in agreement with the results of [10, 47]..

We consider a two-parameter model of the source evolution $\xi(z)$, the
comoving CR power density at redshift $z$ relative to its value today,
consisting of a simple power law and an exponential cutoff

$\displaystyle\xi_{m,z_{0}}(z)=\begin{cases}(1+z)^{m}&z\leq z_{0}\\\
(1+z_{0})^{m}e^{-(z-z_{0})}&z>z_{0}\end{cases},$ (1)

where $-7\leq m\leq 7$ and $1\leq z_{0}\leq 5$. This simple parametrization
sufficiently captures the qualitative features of many observed source
evolutions considered when modeling UHECRs and the neutrinos they produce.
Additionally, several observationally-informed source evolutions are
considered. These include a star formation rate (SFR) evolution [13],

$\displaystyle\xi_{\mathrm{SFR}}(z)\propto\frac{(1+z)^{3.26}}{1+[(1+z)/2.59]^{5.68}},$
(2)

an active galactic nuclei (AGN) evolution [14],

$\displaystyle\xi_{\mathrm{AGN}}(z)\propto\begin{cases}(1+z)^{5}&z\leq 1.7\\\
(1+1.7)^{5}&1.7<z\leq 2.7\\\ (1+1.7)^{5}e^{-(z-2.7)}&z>2.7\end{cases},$ (3)

and a gamma-ray burst (GRB) evolution [15, 16],

$\displaystyle\xi_{\mathrm{GRB}}(z)\propto\frac{(1+z)^{1.5}}{\left[(1+z)^{-34}+\left(\frac{1+z}{5160}\right)^{3}+\left(\frac{1+z}{9}\right)^{35}\right]^{0.1}}.$
(4)

Throughout this work we fit the observed spectrum and composition data of
Auger [17, 18, 19, 1], adopting a $+20\%$ shift of the Auger energy scale and
a $-10$ g/cm2 average shift of $\langle X_{\mathrm{max}}\rangle$ following
[11]. The goodness-of-fit is determined by calculating a combined $\chi^{2}$
to the UHECR spectrum and the first two moments of the depth of shower maximum
distributions, $\langle X_{\mathrm{max}}\rangle$ and
$\sigma\left(X_{\mathrm{max}}\right)$ [mapped into $\langle\ln{A}\rangle$ and
$\mathrm{V}(\ln{A})$, where $A$ is the CR mass number, according to the
parametrization of [20]].

The baseline model for a given source evolution is determined by fitting the
UHECR spectrum above $10^{17.5}$ eV and composition data above $10^{17.8}$ eV,
assuming a single-mass injection into the source environment. Once the best-
fit model is determined all source parameters of the baseline model are fixed
for the remainder of the analysis. The source parameters for baseline models
assuming an observationally-informed source evolution can be found in Appendix
C.

The pure-proton population is given its own set of source parameters and an
injected spectrum with a maximum rigidity in $10-1000$ EeV range. Importantly,
the average number of interactions before escape is a parameter of the model,
allowing our analysis to capture both the possibility of significant source
interactions and of a “naked” accelerator, free of any significant source
environment. In the latter case, CRs produced by the pure-proton population
only experience interactions during extragalactic propagation, as was explored
in [7, 8].

The relative contribution of the two populations is set by a parameter
$f_{pp}$,

$\displaystyle
f_{pp}=\frac{\int_{E_{\mathrm{ref}}}^{\infty}E\phi_{pp}dE}{\int_{E_{\mathrm{ref}}}^{\infty}E(\phi_{pp}+\phi_{\mathrm{BL}})dE},$
(5)

controlling the fraction of energy escaping both source populations produced
by the pure-proton population, where $E_{\mathrm{ref}}=10^{17}$ eV,
$\phi_{pp}$ is the escaping spectrum produced by the pure-proton population,
and $\phi_{\mathrm{BL}}$ is the escaping spectrum produced by the baseline
population.

Figure 1: Correlation between the observed proton fraction above $30$ EeV,
$f_{p}^{\mathrm{obs}}$, and the $1$ EeV neutrino flux, $\phi_{18}$. Each point
represents a separate model realization. Gray points are model realizations
excluded by multimessenger constraints. A power law source evolution was
assumed for each model realization with the value of $m$ indicated by the
color and $z_{0}=2$. The Sibyll2.3c hadronic interaction model (HIM) was
assumed. Figure 2: The range to which the power-law index of the source
evolution can be constrained, $\Delta m$, for a given measurement of the
observed proton fraction $f_{p}^{\mathrm{obs}}$ above $30$ EeV and $1$ EeV
cosmogenic neutrino flux, marginalizing over the cutoff redshift $z_{0}$ and
assuming Sibyll2.3c. Contours indicate the lower-bound on the power law index
for a given measurement. The decrease in $\Delta m$ below
$m_{\mathrm{min}}=-7$ is due to the finite range of negative values of $m$
explored. White regions indicate combinations of the observed proton fraction
and neutrino flux that are either incompatible with multimessenger data, are
not realizable physically, or require a source evolution with $m<-7$. Magenta
stars indicate the predicted values for the best-fit models assuming SFR, AGN,
and GRB source evolutions. Also indicated are current $90\%$ confidence level
(CL) upper-limits on the neutrino flux for Auger [21] (horizontal solid light-
blue line) and IceCube [22] (horizontal solid orange line), as well as, $90\%$
CL limit forecasts from [23] for a variety of ongoing and future neutrino
experiments [24, 25, 26, 27, 28, 29, 23, 30] (horizontal dotted colored
lines). The Auger measurement [4, 2] of the observed proton fraction above
$30$ EeV is also shown for Sibyll2.3c (dashed vertical light blue line
indicates the $1\sigma$ upper-limit). Also indicated are the $90\%$ CL limit
forecasts for a variety of ongoing and future UHECR experiments [31, 32, 33]
(vertical dotted colored lines; details in Appendix A).

To explore the range of multimessenger signals which can be produced by the
pure-proton population for each source evolution and value of $f_{pp}$, all
model parameters of this population are randomly sampled, all parameters of
the baseline population are held fixed 222This was done for computational
efficiency, as we found that also refitting the source parameters of the
baseline model had a negligible effect on our results., and only those
controlling properties of the Galactic CR spectrum (specifically its
composition, spectral index, cutoff energy, and normalization) are tuned to
obtain the best-fit to the UHECR spectrum and composition data above $10^{18}$
eV. Once all parameters have been set, several criteria are used to determine
whether the resulting multimessenger signals are compatible with
multimessenger data. First, we require that the fit to UHECR spectrum and
composition data result in a $\chi^{2}/\mathrm{ndf}<5$. This cut was chosen to
ensure a standard on the absolute quality of the fit, while also accommodating
the varying quality of fit possible for a baseline population alone assuming
different hadronic interaction models (HIMs) and source evolutions [7].
Second, pure-proton models which degrade the quality of fit by more than
$3\sigma$ compared to the baseline model alone are considered to be in
conflict with UHECR data 333We follow the PDG [48, 49] defining the number of
sigma from the best-fit as
$N_{\sigma}^{\prime}=S^{-1}\sqrt{\chi^{2}_{\mathrm{model}}-\chi^{2}_{\mathrm{min}}}$,
where $S=\sqrt{\chi^{2}_{\mathrm{min}}/N_{\mathrm{dof}}}$ is the scale factor
introduced to enlarge the uncertainties to account for a
$\chi^{2}_{\mathrm{min}}/N_{\mathrm{dof}}>1$ [49], $\chi^{2}_{\mathrm{model}}$
is the $\chi^{2}$ for a given model, $\chi^{2}_{\mathrm{min}}$ is the
$\chi^{2}$ of the best-fit model, and $N_{\mathrm{dof}}$ is the number of
degrees of freedom. In order to calculate the number of sigma a model has
improved (or degraded) the fit compared to the baseline model we define
$N_{\sigma}^{\mathrm{rel}}=\operatorname{sgn}(\chi^{2}_{\mathrm{model}}-\chi^{2}_{\mathrm{BL}})S^{-1}\sqrt{|\chi^{2}_{\mathrm{model}}-\chi^{2}_{\mathrm{BL}}|}$,
where
$\chi^{2}_{\mathrm{min}}=\min(\chi^{2}_{\mathrm{model}},\chi^{2}_{\mathrm{BL}})$
in $S$.. Third, models which produce more than $4.74$ neutrinos above
$10^{15.9}$ eV are rejected at $99\%$ confidence level (CL) [36] as they
violate constraints from IceCube [22, 37]. Finally, we consider limits on the
gamma-ray flux at GeV–TeV energies, from Fermi-LAT [38, 39], and at EeV
energies from Auger [40, 41, 42], but find that no models compatible with
other multimessenger constraints are capable of violating them. This
combination of constraints limits $m\leq+6$.

It is possible that some realizations of this model will have parameters which
imply a large source with a strong magnetic field. In this case, pions and
muons produced in the environment suffer significant synchrotron losses before
decaying, effectively cutting-off the resulting neutrino spectrum. We find
that excluding such model realizations does not change our results. However
if, in reality, UHECR environments are in a regime where synchrotron losses
are significant, only the results of Section III.1 would be applicable.

##  III Results

###  III.1 Cosmogenic-only case

For a fixed source evolution we find the flux of neutrinos at $1$ EeV,
$\phi_{18}$, to have a strong correlation with the observed proton fraction
above $30$ EeV, $f^{\mathrm{obs}}_{p}$, so that $\phi_{18}\propto
f^{\mathrm{obs}}_{p}$ (see Fig. 1), as was reported by [8]. To capture the
dispersion of this correlation we find the maximum and minimum values of the
observed proton fraction-to-$1$ EeV neutrino flux ratio, $r_{p\nu,18}$, among
all models compatible with multimessenger data. By construction, all models
must then obey $r^{\mathrm{min}}_{p\nu,18}\leq r_{p\nu,18}\leq
r^{\mathrm{max}}_{p\nu,18}$. For a fixed source evolution, this fact allows a
constraint to be placed on either $f^{\mathrm{obs}}_{p}$ or $\phi_{18}$ if the
other quantity is known according to:

$\displaystyle f^{\mathrm{obs,min}}_{p}$
$\displaystyle=\min\left(r_{p\nu,18}^{\mathrm{min}}\tilde{\phi}_{18},F_{p}^{\mathrm{obs,max}}\right)~{},$
(6) $\displaystyle f^{\mathrm{obs,max}}_{p}$
$\displaystyle=\min\left(r_{p\nu,18}^{\mathrm{max}}\tilde{\phi}_{18},F_{p}^{\mathrm{obs,max}}\right)~{},$
(7) $\displaystyle\phi_{18}^{\mathrm{min}}$
$\displaystyle=\min\left(\frac{\tilde{f_{p}}^{\mathrm{obs}}}{r_{p\nu,18}^{\mathrm{max}}},\Phi^{\mathrm{max}}_{18}\right)~{},$
(8) $\displaystyle\phi_{18}^{\mathrm{max}}$
$\displaystyle=\min\left(\frac{\tilde{f_{p}}^{\mathrm{obs}}}{r_{p\nu,18}^{\mathrm{min}}},\Phi^{\mathrm{max}}_{18}\right)~{},$
(9) where $\displaystyle\tilde{f_{p}}^{\mathrm{obs}}$
$\displaystyle=\min\left(f_{p}^{\mathrm{obs}},F_{p}^{\mathrm{obs,max}}\right)$
(10) and $\displaystyle\tilde{\phi}_{18}$
$\displaystyle=\min\left(\phi_{18},\Phi^{\mathrm{max}}_{18}\right)$ (11)

are the observed proton fraction and $1$ EeV neutrino flux truncated at their
maximum realizable values compatible with multimessenger data,
$F_{p}^{\mathrm{obs,max}}$ and $\Phi^{\mathrm{max}}_{18}$ 444A similar
analysis could be carried out by considering the integral UHE photon flux
above $10$ EeV, as this observable also correlates well with both the observed
proton fraction and the neutrino flux. However, due to the short attenuation
length of UHE photons, the results are very sensitive to the assumed distance
to the nearest source. For this reason we leave such an analysis for future
work..

To be conservative we subtract off the baseline source population’s
contribution to $f_{p}^{\mathrm{obs}}$ and $\phi_{18}$, but in principle
either of these values may have a “floor” set by the baseline source
population depending on the true evolution of these sources. With these
constraints in hand one can determine the range of source evolutions

$\displaystyle\Delta m=m_{\mathrm{max}}-m_{\mathrm{min}}$ (12)

which satisfy ${f_{p}^{\mathrm{obs,min}}\leq f_{p}^{\mathrm{obs}}\leq
f_{p}^{\mathrm{obs,max}}}$ and
${\phi_{18}^{\mathrm{min}}\leq\phi_{18}\leq\phi_{18}^{\mathrm{max}}}$ for a
particular $(f_{p}^{\mathrm{obs}},\phi_{18})$. These constraints are shown in
Fig. 2 assuming the Sibyll2.3c [44] HIM. Analogous figures showing the results
when assuming the EPOS-LHC [45] HIM can be found in Appendix B.

It is clear that for the cosmogenic-only case, the dispersion of the
correlation between the observed proton fraction and $1$ EeV neutrino flux is
small (as evidenced by Fig. 1 and the small value of $\Delta m$ over most of
the parameter space in Fig. 2). This enables UHECR and neutrino measurements
to jointly measure the evolution of such a population of UHE proton sources.
For example, Fig. 2 shows if AugerPrime [31] measures the proton fraction
above $30$ EeV to be $3\%$, then assuming the Radio Echo Telescope for
Neutrinos (RET/̄N) [27] measures the $1$ EeV neutrino flux to be $10^{-8.5}$
GeV$/$cm${}^{2}/$s$/$sr, the source evolution will be constrained $\Delta
m\lesssim 3$ and $m>3$.

(a)  (b)

Figure 3: Same as Fig. 2 but for the general case where protons may or may not
have significant interactions in the source environment, and for measuring the
neutrino flux at $1$ EeV (left) and $10$ EeV (right). Dashed contours
indicating the upper-bound on the source evolution’s power law index are also
shown for the $1$ EeV case. Additional $90\%$ CL limit forecasts for $10$ EeV
neutrino sensitivity from [23] are shown in the top panels for a variety of
ongoing and future neutrino experiments [24, 25, 26, 27, 28, 29, 23, 30, 46]
(dotted colored lines).

Fig. 2 also shows that a $1$ EeV neutrino detection alone will be able to
constrain the source evolution of such UHE proton sources. For example, Fig. 2
shows that if the Radio Neutrino Observatory in Greenland (RNO/̄G) [25, 26]
detects a $1$ EeV neutrino then source evolutions with $m<3$ will be excluded.
This is simply because these source evolutions are not capable of producing a
large enough $1$ EeV neutrino flux to be detectable by RNO/̄G without
violating current CR composition constraints.

Finally, we comment that Fig. 2 can be used to provide benchmark sensitivities
for future UHECR and neutrino experiments. Measurement, or constraint, on the
flux of one messenger places an upper-bound on the flux of the other
messenger. For example, Fig. 2 shows if the observed proton fraction is
constrained to be less than $1\%$ then that would imply the $1$ EeV neutrino
flux is less than $10^{-8}$ GeV$/$cm${}^{2}/$s$/$sr. A similar statement for
the observed proton fraction is possible under very mild assumptions about the
source evolution. For example, the $1$ EeV neutrino flux were constrained to
be less than $10^{-10}$ GeV$/$cm${}^{2}/$s$/$sr then the observed proton
fraction must be less than $5\%$ for positive source evolutions.

###  III.2 General case

In the more general case, where we allow for the possibility of a significant
number of interactions in the environment host to the accelerator, the
correlation between the flux of neutrinos at $1$ EeV and the observed proton
fraction above $30$ EeV is weaker. This weaker correlation amounts to a wider
dispersion and, therefore, a larger range of $r_{p\nu,18}$ values. The
resulting constraints for this more general case are shown in Fig. 3a. In
particular, whereas in the cosmogenic-only case one can always effectively
measure the source evolution, in the general case this is only possible for
certain $(f_{p}^{\mathrm{obs}},\phi_{18})$ combinations. More generally, it
may only be possible to set an upper- or lower-bound on the value of $m$ using
the $1$ EeV neutrino flux.

However, several planned and proposed neutrino experiments in the near future
will have peak sensitivity in the $10$ EeV range, rather than the $1$ EeV
range. Given the large dispersion of the $\phi_{18}$–$f^{\mathrm{obs}}_{p}$
correlation in the general case, it is worthwhile to explore how these higher
energy neutrino observatories will be able to provide insight into the
evolution of UHECR sources. Therefore, we also consider the correlation
between the $10$ EeV neutrino flux, $\phi_{19}$, and the observed proton
fraction. Similar to the $1$ EeV case, we find that
$f^{\mathrm{obs}}_{p}\propto\phi_{19}$ but with a large dispersion. Defining
the observed proton fraction-to-$10$ EeV neutrino flux ratio, $r_{p\nu,19}$,
we can constrain the realizable range of $f^{\mathrm{obs}}_{p}$ and
$\phi_{19}$ compatible with multimessenger data analogously to equations
(6)-(9). The constraints based on the $10$ EeV neutrino flux can be found in
Fig. 3b.

From Fig. 3b we see that it is only possible to constrain the source evolution
for some combinations of $f_{p}^{\mathrm{obs}}$ and $\phi_{19}$ – similar to
the situation in Fig. 3a. However even in the general case, where UHE protons
have a significant number of interactions in the source environment,
measurement of the source evolution may be possible by combining the observed
proton fraction with measurement of the $1$ EeV and $10$ EeV neutrino fluxes.
For example, let’s assume the observed proton fraction was measured to be
$1\%$, the $1$ EeV neutrino flux were measured to be $10^{-10}$
GeV$/$cm${}^{2}/$s$/$sr, and the $10$ EeV neutrino flux were measured to be
$10^{-9.5}$ GeV$/$cm${}^{2}/$s$/$sr. Then from Fig. 3 we see that these
measurements would allow us to infer that $m\lesssim 3$ (driven by the $1$ EeV
neutrino flux measurement) and $m\gtrsim 0$ (driven by the $10$ EeV neutrino
flux measurement). However, it is important to note that not all combinations
of these observables yield strong constraints on the source evolution.

Tantalizingly, both Figs. 2 and 3 show that the best-fit models assuming
astrophysical source evolutions predict an observed proton fraction of
$\gtrsim 1\%$ regardless of the particular astrophysical scenario. These best-
fit models also predict a $1$ EeV neutrino flux that will be detectable by the
next generation of neutrino experiments. This prediction suggests that
discovery of such a UHE proton component and its neutrino flux – and therefore
a measurement of the evolution of its sources – by the next generation of
UHECR and neutrino experiments is a realistic possibility.

##  IV Summary

In this study we have considered the prospects for constraining the evolution
of a population of pure-proton sources by combining UHECR and neutrino data.
Neither of these messengers can determine the source evolution alone. However,
we have found that near-future UHECR and neutrino detectors could
realistically place strong constraints on the evolution of such a population.

In the case of a purely cosmogenic flux of neutrinos, near-future detectors
will constrain the source evolution as long as the proton fraction above $30$
EeV is $\gtrsim 10^{-4}$ and the neutrino flux at $1$ EeV is $\gtrsim
10^{-10}$ GeV/cm2/s/sr – a requirement favored by the best-fits to the UHECR
spectrum and composition data we find. In this case, the $1$ EeV neutrino flux
and the observed proton fraction can be combined to constrain the source
evolution to a narrow range of possibilities.

In the case that source interactions result in a significant astrophysical
neutrino flux, more information may be required to constrain the source
evolution. We have shown that by combining the observed proton fraction with
measurements of the $1$ EeV and $10$ EeV neutrino flux, future detectors may
be able to constrain the source evolution’s power-law index $m$ to a limited
interval. Even if the neutrino flux is only measured at one of these energies,
an upper- or lower-bound may still be placed on $m$.

Importantly, even if source interactions are significant, best-fit models
still predict that this proton component and its secondary neutrinos at $1$
EeV will be detectable by the next generation of UHECR and neutrino
experiments for many of the evolutions often considered for the sources of
UHECRs.

Our results underscore the complementarity of neutrino and UHECR detectors, as
well as, the need for a next-generation of detectors for both of these
messengers. Perhaps epitomizing the strength of multimessenger astrophysics,
our results show how combining neutrino and UHECR observations provides access
to a quantity inaccessible by either of these messengers alone.

###### Acknowledgements.

We thank Foteini Oikonomou, David Seckel, and Domenik Ehlert for useful
feedback on our analysis. The research of MSM is supported by the NSF MPS-
Ascend Postdoctoral Award #2138121. The research of SW is supported by NSF
Awards #2111232 and #2033500 and NASA grants 80NSSC20K0925, 80NSSC22K1519 and
80NSSC21M0116.

## References

  * Yushkov [2020] A. Yushkov (Pierre Auger), PoS ICRC2019, 482 (2020).
  * Aab  _et al._ [2014a] A. Aab _et al._ (Pierre Auger), Phys. Rev. D 90, 122005 (2014a), arXiv:1409.4809 [astro-ph.HE] .
  * Abraham  _et al._ [2010] J. Abraham _et al._ (Pierre Auger), Phys. Rev. Lett. 104, 091101 (2010), arXiv:1002.0699 [astro-ph.HE] .
  * Bellido [2018] J. Bellido (Pierre Auger), PoS ICRC2017, 506 (2018).
  * Aab  _et al._ [2014b] A. Aab _et al._ (Pierre Auger), Phys. Rev. D 90, 122006 (2014b), arXiv:1409.5083 [astro-ph.HE] .
  * Abdul Halim  _et al._ [2022] A. Abdul Halim _et al._ (Pierre Auger), (2022), arXiv:2211.02857 [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] .
  * van Vliet  _et al._ [2019] A. van Vliet, R. Alves Batista, and J. R. Hörandel, Phys. Rev. D 100, 021302 (2019), arXiv:1901.01899 [astro-ph.HE] .
  * Møller  _et al._ [2019] K. Møller, P. B. Denton, and I. Tamborra, JCAP 05, 047, arXiv:1809.04866 [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] .
  * 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] .
  * Note [1] Other parametrizations of the source photon field, including a broken power-law, were explored but did not have a significant effect on our results, in agreement with the results of [10, 47].
  * Robertson  _et al._ [2015] B. E. Robertson, R. S. Ellis, S. R. Furlanetto, and J. S. Dunlop, Astrophys. J. Lett. 802, L19 (2015), arXiv:1502.02024 [astro-ph.CO] .
  * Stanev [2008] T. Stanev, (2008), arXiv:0808.1045 [astro-ph] .
  * Kistler  _et al._ [2008] M. D. Kistler, H. Yuksel, J. F. Beacom, and K. Z. Stanek, Astrophys. J. Lett. 673, L119 (2008), arXiv:0709.0381 [astro-ph] .
  * Yuksel  _et al._ [2008] H. Yuksel, M. D. Kistler, J. F. Beacom, and A. M. Hopkins, Astrophys. J. Lett. 683, L5 (2008), arXiv:0804.4008 [astro-ph] .
  * Aab  _et al._ [2020a] A. Aab _et al._ (Pierre Auger), Phys. Rev. D 102, 062005 (2020a), arXiv:2008.06486 [astro-ph.HE] .
  * Aab  _et al._ [2020b] A. Aab _et al._ (Pierre Auger), Phys. Rev. Lett. 125, 121106 (2020b), arXiv:2008.06488 [astro-ph.HE] .
  * Verzi [2020] V. Verzi (Pierre Auger), PoS ICRC2019, 450 (2020).
  * Abreu  _et al._ [2013] P. Abreu _et al._ (Pierre Auger), JCAP 02, 026, arXiv:1301.6637 [astro-ph.HE] .
  * Aab  _et al._ [2019] A. Aab _et al._ (Pierre Auger), JCAP 10, 022, arXiv:1906.07422 [astro-ph.HE] .
  * Aartsen  _et al._ [2018] M. G. Aartsen _et al._ (IceCube), Phys. Rev. D 98, 062003 (2018), arXiv:1807.01820 [astro-ph.HE] .
  * Ackermann  _et al._ [2022] M. Ackermann _et al._ , JHEAp 36, 55 (2022), arXiv:2203.08096 [hep-ph] .
  * Allison  _et al._ [2020] P. Allison _et al._ (ARA), Phys. Rev. D 102, 043021 (2020), arXiv:1912.00987 [astro-ph.HE] .
  * Aguilar  _et al._ [2021a] J. A. Aguilar _et al._ (RNO-G), JINST 16 (03), P03025, arXiv:2010.12279 [astro-ph.IM] .
  * Aguilar  _et al._ [2021b] J. A. Aguilar _et al._ (RNO-G), PoS ICRC2021, 001 (2021b).
  * Prohira  _et al._ [2020] S. Prohira _et al._ , Phys. Rev. Lett. 124, 091101 (2020), arXiv:1910.12830 [astro-ph.HE] .
  * Brown  _et al._ [2021] A. M. Brown, M. Bagheri, M. Doro, E. Gazda, D. Kieda, C. Lin, N. Otte, I. Taboada, and A. Wang, PoS ICRC2021, 1179 (2021).
  * Wissel  _et al._ [2020] S. Wissel _et al._ , JCAP 11, 065, arXiv:2004.12718 [astro-ph.IM] .
  * Á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] .
  * Aab  _et al._ [2016] A. Aab _et al._ (Pierre Auger), (2016), arXiv:1604.03637 [astro-ph.IM] .
  * Olinto  _et al._ [2021] A. V. Olinto _et al._ (POEMMA), JCAP 06, 007, arXiv:2012.07945 [astro-ph.IM] .
  * Hörandel [2021] J. R. Hörandel (GCOS), PoS ICRC2021, 027 (2021), arXiv:2203.01127 [astro-ph.HE] .
  * Note [2] This was done for computational efficiency, as we found that also refitting the source parameters of the baseline model had a negligible effect on our results.
  * Note [3] We follow the PDG [48, 49] defining the number of sigma from the best-fit as $N_{\sigma}^{\prime}=S^{-1}\sqrt{\chi^{2}_{\mathrm{model}}-\chi^{2}_{\mathrm{min}}}$, where $S=\sqrt{\chi^{2}_{\mathrm{min}}/N_{\mathrm{dof}}}$ is the scale factor introduced to enlarge the uncertainties to account for a $\chi^{2}_{\mathrm{min}}/N_{\mathrm{dof}}>1$ [49], $\chi^{2}_{\mathrm{model}}$ is the $\chi^{2}$ for a given model, $\chi^{2}_{\mathrm{min}}$ is the $\chi^{2}$ of the best-fit model, and $N_{\mathrm{dof}}$ is the number of degrees of freedom. In order to calculate the number of sigma a model has improved (or degraded) the fit compared to the baseline model we define $N_{\sigma}^{\mathrm{rel}}=\operatorname{sgn}(\chi^{2}_{\mathrm{model}}-\chi^{2}_{\mathrm{BL}})S^{-1}\sqrt{|\chi^{2}_{\mathrm{model}}-\chi^{2}_{\mathrm{BL}}|}$, where $\chi^{2}_{\mathrm{min}}=\mathop{min}\displaylimits(\chi^{2}_{\mathrm{model}},\chi^{2}_{\mathrm{BL}})$ in $S$.
  * Feldman and Cousins [1998] G. J. Feldman and R. D. Cousins, Phys. Rev. D 57, 3873 (1998), arXiv:physics/9711021 .
  * Aartsen  _et al._ [2021] M. G. Aartsen _et al._ (IceCube), Nature 591, 220 (2021), [Erratum: Nature 592, E11 (2021)], arXiv:2110.15051 [hep-ex] .
  * Ackermann  _et al._ [2015] M. Ackermann _et al._ (Fermi-LAT), Astrophys. J. 799, 86 (2015), arXiv:1410.3696 [astro-ph.HE] .
  * Ackermann  _et al._ [2016] M. Ackermann _et al._ (Fermi-LAT), Phys. Rev. Lett. 116, 151105 (2016), arXiv:1511.00693 [astro-ph.CO] .
  * Rautenberg [2021] J. Rautenberg (Pierre Auger), PoS ICRC2019, 398 (2021).
  * Abreu  _et al._ [2021a] P. Abreu _et al._ (Pierre Auger), PoS ICRC2021, 373 (2021a).
  * Abreu  _et al._ [2022] P. Abreu _et al._ (Pierre Auger), Astrophys. J. 933, 125 (2022), arXiv:2205.14864 [astro-ph.HE] .
  * Note [4] A similar analysis could be carried out by considering the integral UHE photon flux above $10$ EeV, as this observable also correlates well with both the observed proton fraction and the neutrino flux. However, due to the short attenuation length of UHE photons, the results are very sensitive to the assumed distance to the nearest source. For this reason we leave such an analysis for future work.
  * Fedynitch  _et al._ [2019] A. Fedynitch, F. Riehn, R. Engel, T. K. Gaisser, and T. Stanev, Phys. Rev. D 100, 103018 (2019), arXiv:1806.04140 [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] .
  * Abarr  _et al._ [2021] Q. Abarr _et al._ (PUEO), JINST 16 (08), P08035, arXiv:2010.02892 [astro-ph.IM] .
  * Fiorillo  _et al._ [2021] D. F. G. Fiorillo, A. Van Vliet, S. Morisi, and W. Winter, JCAP 07, 028, arXiv:2103.16577 [astro-ph.HE] .
  * Zyla  _et al._ [2020] P. A. Zyla _et al._ (Particle Data Group), PTEP 2020, 083C01 (2020).
  * Rosenfeld [1975] A. H. Rosenfeld, Ann. Rev. Nucl. Part. Sci. 25, 555 (1975).
  * Note [5] For $Y=X_{\rm max}$ we compared the sensitivities for two-component normal and Gumbel distributions with the same merit factor and the same ratio of standard deviations, $\sigma(Y_{\mathrm{He}})/\sigma(Y_{p})$. The resulting sensitivities are of similar magnitude, about 40% better for the normal distribution than for the skewed Gumbel distribution for proton and helium primaries.
  * Abreu  _et al._ [2021b] P. Abreu _et al._ (Pierre Auger), PoS ICRC2021, 324 (2021b).
  * Coleman  _et al._ [2023] A. Coleman _et al._ , Astropart. Phys. 147, 102794 (2023), arXiv:2205.05845 [astro-ph.HE] .
  * Arbeletche and de Souza [2020] L. B. Arbeletche and V. de Souza, Astropart. Phys. 116, 102389 (2020), arXiv:1903.03174 [astro-ph.HE] .
  * Note [6] A non-finite value of $r_{\mathrm{size}}$ indicates that the highest energy CRs considered in our model enter the quasi-ballistic diffusion regime before escaping the source.
  * Aartsen  _et al._ [2020] M. G. Aartsen _et al._ (IceCube), Phys. Rev. Lett. 125, 121104 (2020), arXiv:2001.09520 [astro-ph.HE] .
  * Bird  _et al._ [1995] D. J. Bird _et al._ , Astrophys. J. 441, 144 (1995), arXiv:astro-ph/9410067 .
  * Matthews [2022] J. Matthews (Telescope Array), Results from the Telescope Array, Talk at the ”Cosmic Rays in the Multimessenger Era” workshop, APC Paris (2022).
  * Abreu  _et al._ [2021c] P. Abreu _et al._ (Pierre Auger), Eur. Phys. J. C 81, 966 (2021c), arXiv:2109.13400 [astro-ph.HE] .

##  Appendix A Estimation of UHECR experimental sensitivity to
$f^{\mathrm{obs}}_{p}$

To estimate the ability of ongoing and future UHECR experiments to constrain
$f^{\mathrm{obs}}_{p}$ we consider two possible cases. First, we consider the
optimistic case where measurement of the proton flux is background-free (i.e.
the experiment can perfectly separate protons from observed heavier nuclei).
In this case, the strongest limit would be set if no proton events are
detected, then the $90\%$ CL upper-limit is given by

$\displaystyle f^{\mathrm{obs,\text{zero
bg.}}}_{p}=\frac{\mathrm{FC}(0,0)}{N_{\mathrm{evts}}}~{},$ (13)

where $\mathrm{FC}(0,0)=2.44$ is $90\%$ CL Feldman-Cousins upper-limit for
zero observed and background events, and $N_{\mathrm{evts}}$ corresponds to
the total number of CR events detected above $30$ EeV.

More generally, we consider the case where the proton flux cannot be perfectly
separated from heavier nuclei, but a considerable overlap exists in the
distributions of the mass-sensitive experimental variable $Y$ (e.g. the shower
maximum $X_{\mathrm{max}}$). Conservatively, we assume a background of helium
events, since this is the nucleus which is most difficult to separate from
protons. We define the proton-fraction sensitivity as the minimum fraction
with which the null-hypothesis (pure helium flux) can be rejected at a
confidence level of $90\%$ CL. This minimum fraction is determined by
repeatedly sampling $Y$ distributions of $N_{\mathrm{evts}}$ helium events. To
each of these simulated data sets we fit a two-component (proton and helium)
model. The 90% quantile of the obtained proton fraction distribution for a
pure helium composition defines then the sensitivity for the proton fraction.
Figure 4 shows the resulting sensitivity as a function of these two variables
assuming $Y$ is normally distributed 555For $Y=X_{\rm max}$ we compared the
sensitivities for two-component normal and Gumbel distributions with the same
merit factor and the same ratio of standard deviations,
$\sigma(Y_{\mathrm{He}})/\sigma(Y_{p})$. The resulting sensitivities are of
similar magnitude, about 40% better for the normal distribution than for the
skewed Gumbel distribution for proton and helium primaries..

For this purpose, we obtain the expected number of events above $30$ EeV for
each experiment by multiplying its expected exposure by the integral CR flux
above $30$ EeV according to the Auger spectrum model [51]. Experimental
exposures, $\mathcal{E}$, were taken from [52].

Our assumptions about the mass sensitive variable $Y$ vary depending on the
experiment. For POEMMA we use $X_{\mathrm{max}}$ as the mass-sensitive
variable $Y$, distributed according to a generalized Gumbel distribution with
parameter values from [53]. For AugerPrime and GCOS we assume that
$\sigma(Y_{\mathrm{He}})/\sigma(Y_{p})=\sigma(X_{\mathrm{max,He}})/\sigma(X_{\mathrm{max,}p})\simeq
0.71$ to determine the separation of the proton and helium distributions for a
given merit factor. For AugerPrime we use the published proton-helium merit
factors [31], while for GCOS we consider a high- and low-resolution design
with proton-helium merit factors of $0.7$ (optimistic) and $0.3$
(conservative), respectively. In both cases we assume $Y$ follows a normal
distribution.

Figure 4: Estimate of the $90\%$ CL proton sensitivity in the presence of
helium background above $3\times 10^{19}$ eV for a given proton-helium merit
factor, $f_{\mathrm{MF}}$, and observed number of events, $N_{\mathrm{evts}}$.
Further details can be found in the text. The values of $f_{\mathrm{MF}}$ and
$N_{\mathrm{evts}}$ assumed for each of the experiments considered are also
indicated.

##  Appendix B Results for EPOS-LHC

Figure 5: Same as Fig. 2 but using EPOS-LHC as the HIM. The Auger measurement
[4, 2] of the observed proton fraction above $30$ EeV is shown for EPOS-LHC
(vertical light blue line, central value is solid and $1\sigma$ errors are
dashed).

Figure 5 shows the range of power-law indices allowed for a particular
combination of the observed proton fraction and the $1$ EeV neutrino flux for
the cosmogenic-only case, when assuming the EPOS-LHC HIM. Compared with Fig.
2, EPOS-LHC generally allows for a larger proton fraction, due to the fact
that it infers the composition to be lighter from air shower data. This
difference in interpretation of air shower data also leads Auger to infer a
non-zero proton fraction of $\sim 5\%$ [4]. The best-fit models for
astrophysically-informed source evolutions also favor larger proton fractions
at Earth.

The most noticeable difference between Fig. 2 and Fig. 5, though, is the much
smaller range of $\Delta m$ values in the EPOS-LHC case. This is due to the
narrower dispersion in the $f_{p}^{\mathrm{obs}}-\phi_{18}$ correlation,
driven by the fact that fits to UHECR data assuming EPOS-LHC generally have a
poorer quality and that we require $\chi^{2}/ndf<5$. This results in models
assuming EPOS-LHC effectively being more constrained than those assuming
Sibyll2.3c.

Figure 6 shows the range of power-law indices allowed for particular
combinations of the observed proton fraction with the $1$ EeV and $10$ EeV
neutrino flux for the general case, where protons may or may not have
significant interactions in the source environment, for EPOS-LHC. Similar to
the general case under Sibyll2.3c (see Fig. 3) not all combinations of these
observables lead to constraints on the source evolution, due to the large
dispersion in $r_{p\nu,18}$ and $r_{p\nu,19}$. However, for some combinations
of these observables it is possible to place an upper- or lower-bound on $m$,
and combining measurements of all three observables can result in strong
constraints on the source evolution in some cases.

(a)  (b)

Figure 6: Same as Fig. 3 but using EPOS-LHC as the HIM. The Auger measurement
[4, 2] of the observed proton fraction above $30$ EeV is shown for EPOS-LHC
(vertical light blue line, central value is solid and $1\sigma$ errors are
dashed).

##  Appendix C Best-fit baseline model parameters

Table 1 shows the best-fit source parameters for baseline models assuming
either a SFR, GRB, or AGN source evolution. The parameters in this table are
as follows: $\gamma_{\mathrm{inj}}$ the spectral index ($J\propto
E^{\gamma_{\mathrm{inj}}}$) at injection into the source environment;
$R_{\mathrm{max}}$ is the maximum rigidity of the injected CR spectrum, where
the spectrum is cutoff exponentially; $r_{\mathrm{esc}}$ is the ratio of the
escape and totalinteraction times for a $10^{19}$ eV iron nucleus;
$r_{g\gamma}$ is the ratio of the hadronic and photohadronic interaction times
for a $10^{19}$ eV iron nucleus; $R_{\mathrm{diff}}$ is the characteristic
rigidity scale of diffusion in the source’s turbulent magnetic field;
$r_{\mathrm{size}}$ is the ratio of the source’s size to the magnetic field’s
coherence length 666A non-finite value of $r_{\mathrm{size}}$ indicates that
the highest energy CRs considered in our model enter the quasi-ballistic
diffusion regime before escaping the source; $T$ is the black-body temperature
of the ambient photon field surrounding the source; and, $A_{\mathrm{inj}}$ is
the mass number of the CRs injected into the source environment (non-integers
represent the average mass due to a mixture of two consecutive mass numbers in
order for $A_{\mathrm{inj}}$ to be a continuous model parameter).

|  SFR |  GRB |  AGN
---|---|---|---
Parameter | Sibyll2.3c | EPOS-LHC | Sibyll2.3c | EPOS-LHC | Sibyll2.3c | EPOS-LHC
$\gamma_{\mathrm{inj}}$ | $-1.14$ | $-1.71$ | $-1.1$ | $-0.0$ | $-0.99$ | $-0.55$
$\log_{10}(R_{\mathrm{max}}/\mathrm{V})$ | $18.58$ | $18.76$ | $18.65$ | $18.48$ | $18.64$ | $18.58$
$\log_{10}{r_{\mathrm{esc}}}$ | $2.38$ | $3.11$ | $2.64$ | $1.89$ | $2.51$ | $2.33$
$\log_{10}{r_{g\gamma}}$ | $9.99$ | $1.02$ | $9.84$ | $1.26$ | $6.61$ | $1.7$
$\log_{10}(R_{\mathrm{diff}}/\mathrm{V})$ | $17.66$ | $14.0$ | $14.0$ | $14.15$ | $14.0$ | $14.01$
$\tanh(\log_{10}{r_{\mathrm{size}}})$ | $0.8$ | $1.0$ | $1.0$ | $1.0$ | $1.0$ | $1.0$
$T/\mathrm{K}$ | $1800$ | $1003$ | $5000$ | $6013$ | $5002$ | $4007$
$A_{\mathrm{inj}}$ | $32.52$ | $26.26$ | $32.0$ | $25.53$ | $32.0$ | $27.75$
Table 1: Best-fit source parameters for baseline models assuming an
observationally-informed source evolution. Definitions of the source
parameters are given in the text.

##  Appendix D Maximum UHE neutrino flux

Figure 7 shows the maximum neutrino flux realizable by our model while
remaining compatible with multimessenger constraints. The flux shown is the
total neutrino flux produced by the pure-proton source population alone. The
maximum neutrino flux is broken into two cases: 1) cosmogenic-only neutrinos
(dashed lines, corresponding to Section III.1) and 2) both cosmogenic
neutrinos and neutrinos produced inside the source environment (solid lines,
corresponding to Section III.2). At low energies the neutrino flux allowed by
our analysis exceeds the IceCube measurements since this analysis only
excluded models using constraints on the neutrino flux above $10$ PeV, where
no neutrinos have been observed. Neutrinos at lower energies do not effect the
results of our analysis.

As can be seen from Fig. 7, significant interactions in the source environment
primarily contribute to the neutrino flux in the $100$ PeV to $10$ EeV energy
range, below the main peak at $\sim 10$ EeV. Unsurprisingly, the overall
normalization of the flux increases for more positive source evolutions — with
SFR being the least positive and AGN the most positive. Interestingly, while
the AGN and GRB source evolutions are strong enough to saturate current
IceCube limits, the SFR evolution is unable to do so. Finally, as can be seen
comparing Figs. 7a and 7b, EPOS-LHC results in a slightly higher neutrino flux
due to its lighter inference on the UHECR composition data, allowing for a
slightly higher proton fraction.

(a)  (b)

Figure 7: The maximum realizable neutrino flux in each energy bin for models
compatible with multimessenger constraints (N.B. $99\%$ CL neutrino
constraints are used) assuming Sibyll2.3c (left) and EPOS-LHC (right). The
maximum flux for both the cosmogenic-only (dashed lines) and general case
(solid lines) are shown for three observationally-informed source evolutions.
Current $90\%$ CL neutrino limits from IceCube and Auger are shown, along with
measurements of the astrophysical neutrino flux [55, 37].

##  Appendix E Constraints on trans-GZK spectral recovery

Whether the cutoff observed by Auger and TA is truly the end of the UHECR
spectrum remains an open question. The model we have described here explores
the possibility of a pure-proton recovery of the spectrum above the observed
cutoff energy. To quantify when a model has a significant recovery compared to
expectation, we calculate the maximum of the ratio of the model spectrum to
the Auger model spectrum [51] above $10^{20.3}$ eV,

$\displaystyle r_{\mathrm{rec}}=\max_{E>10^{20.3}\text{
eV}}\left(\frac{J_{\mathrm{model}}}{J_{\mathrm{Auger}}}\right)~{}.$ (14)

We consider a model to have a significant recovery over expectation if
$r_{\mathrm{rec}}\geq 50$.

A significant flux of protons above $10^{20.3}$ eV will result in a
significant flux of neutrinos at $10$ EeV. We find that the maximum allowed
recovery in $E^{2}J$ above $10^{20.3}$ eV is well-correlated with the $10$ EeV
neutrino flux. This allows for neutrino flux measurements to constrain the
level of recovery in the UHECR spectrum.

Importantly, this connection relies crucially on the assumption that the
recovery includes a pure-proton component above $10^{20.3}$ eV. It is
reasonable to assume that if a recovery does occur, that its lowest-energy
component be protonic. However, our results are not applicable if this
component falls below $10^{20.3}$ eV or if the recovery is via a pure, heavy
component. In that case, the recovery could be much larger than would be
suggested by the flux of neutrinos at $10$ EeV.

Similarly, the level of possible recovery depends on the assumed distance to
the nearest source in the pure-proton population. However, because we assume a
continuous source distribution to $z=0$, the level of possible recovery in our
model is maximized allowing us to set a conservative upper-bound on the
recovery.

(a)  (b)

Figure 8: Lower limit on the $10$ EeV neutrino flux for various levels of
UHECR proton recovery above $10^{20.3}$ eV (in units of
$E_{0}^{2}J_{0}=\text{eV/km}^{2}\text{/sr/yr}$) as a function of source
evolution. The maximum realizable spectral recovery compatible with
multimessenger data is indicated by the thick black line. The dependence on
HIM is illustrated in the left (Sibyll2.3c) and right (EPOS-LHC) panels.

We find that a protonic recovery in the UHECR spectrum as large as $\gtrsim
10^{15.5}$ eV/km2/sr/yr is compatible with current multimessenger data. Figure
8 shows the minimum $10$ EeV neutrino flux compatible with various levels of
this recovery. Importantly, for positive source evolutions, the next
generation of neutrino detectors will be able to constrain this recovery. By
contrast, the next generation of UHECR observatories will not be able to probe
the peak of this component in general. However, they may be able to determine
whether the spectrum is beginning to recover, as can be seen in Fig. 9, in
some cases. The compatibility of a strong trans-GZK spectral recovery with
existing multimessenger data raises the tantalizing possibility that previous
events measured beyond the observed spectral cutoff by Fly’s Eye [56] and TA
[57] might originate from such a population of sources.

(a)  (b)

Figure 9: The UHECR spectrum (upper panels) and composition (lower panels) for
models maximizing the spectral recovery above $10^{20.3}$ eV for various
levels of $10$ EeV neutrino flux (colored lines, in units of
$E_{0}^{2}\phi_{0}=$ GeV$/$cm${}^{2}/$s$/$sr), assuming a SFR evolution. The
Auger fit to the UHECR spectrum ([51], black dashed line) is shown for
comparison. Results are shown for the Sibyll2.3c (left) and EPOS-LHC (right)
HIMs. Also shown are the Auger spectrum [58] and composition, as well as,
upper-limits on the spectrum at the highest energies (black points and upper-
limits). Projected $84\%$ CL upper-limits on the spectrum above $10^{20.3}$ eV
for GCOS are also shown (green upper-limits) based on a $10^{6}$ km2 sr yr
exposure, given in [52]. Predicted $\langle X_{\mathrm{max}}\rangle$ and
$\sigma\left(X_{\mathrm{max}}\right)$ values for pure-proton and pure-iron
spectra are shown for each HIM (gray lines).