Constraints on the proton fraction of cosmic rays at the highest energies and the consequences for cosmogenic neutrinos and photons

Over the last decade, observations have shown that the mean mass of ultra-high-energy cosmic rays (UHECRs) increases progressively toward the highest energies. However, the precise composition is still unknown, and several theoretical studies hint at the existence of a subdominant proton component up to the highest energies. Motivated by the exciting prospect of performing charged-particle astronomy with ultra-high-energy (UHE) protons we quantify the level of UHE-proton flux that is compatible with present multimessenger observations and the associated fluxes of neutral messengers produced in the interactions of the protons. We study this scenario with numerical simulations of two independent populations of extragalactic sources and perform a fit to the combined UHECR energy spectrum and composition observables, constrained by diffuse gamma-ray and neutrino observations. We find that up to of order $10\%$ of the cosmic rays at the highest energies can be UHE protons, although the result depends critically on the selected hadronic interaction model for the air showers. Depending on the maximum proton energy ($E_\text{max}^\text{p}$) and the redshift evolution of sources, the associated flux of cosmogenic neutrinos and UHE gamma rays can significantly exceed the multimessenger signal of the mixed-mass cosmic rays. Moreover, if $E_\text{max}^\text{p}$ is above the GZK limit, we predict a large flux of UHE neutrinos above EeV energies that is absent in alternate scenarios for the origin of UHECRs. We present the implications and opportunities afforded by these UHE proton, neutrino and photon fluxes for future multimessenger observations.


Introduction
Ultra-high-energy cosmic rays (UHECRs), charged particles of astrophysical origin with energy above ∼ 10 18 eV, are the most energetic cosmic messengers and, as such, probes of the most extreme astrophysical environments.Because of extragalactic and Galactic magnetic fields, their sources remain elusive, even after years of high-precision observation by the latest generation of UHECR detectors, in particular the Pierre Auger Observatory (Auger) and the Telescope Array (TA).
Observations suggest that the composition of UHECRs is surprisingly pure, with each accelerated nuclear species only dominant in a very narrow band of the UHECR spectrum, and the entire spectrum is produced through a carefully tuned combination of the individual peaks (e.g.[1][2][3][4]).The combination of a smooth increase of the average mass and pure composition at all energies implies that the population variance of sources must be remarkably low ( [5]; see also [6]).Under these circumstances, the observed flux cutoff at E CR ≳ 50 EeV is generally predicted to be an effect of the maximum particle energy reachable at the cosmic accelerators.Within this "Peters cycle" [7,8] model of cosmic-ray acceleration with rigiditydependent maximum energy, no light cosmic rays (CRs) are expected at the highest energies.
Nevertheless, the existence of protons or light nuclei at the highest energies, where there are no measurements of composition-sensitive observables with the fluorescence detectors of Auger and TA, cannot be ruled out at present.A very interesting possibility would be the existence of an additional proton-dominated component at the highest energies.Such a flux cannot be easily explained by reprocessing of accelerated UHECR within the source as proposed for extragalactic protons below the ankle, see e.g.[1], but must originate from a secondary population of independent sources that exclusively accelerates protons to ultrahigh energies or where heavier nuclei are efficiently disintegrated before escaping the source region.Motivations for an additional source population come from the expected differences between possible UHECR accelerators, e.g.active galactic nuclei [9] or gamma-ray bursts [10].Such a proton flux does not need to be produced by astrophysical processes necessarily, but could also originate from the decay of heavy dark matter (e.g.[11,12]).
Circumstantial evidence for an additional proton component is provided by an apparent flattening of the increase in observed UHECR mass at E CR ≳ 30 EeV, as reported in an analysis of Auger surface detector data [13,14].This feature could indicate a flux of UHE protons with different spectral index to the bulk of the UHECRs, either from a secondary source population or from a single nearby source [15], but it could also originate from a natural mass limit of the mixed UHECR flux.
Similar two-component models have been previously studied, either in the context of the transition region between Galactic and extragalactic cosmic rays below 10 18.7 eV [16][17][18][19], or similar to the present paper at the highest energies [20,21].
Compared to Ref. [21] we consider a much wider range of proton source parameters.In particular, we study three distinct scenarios for the spectrum injected by the UHE proton sources, of which only one is considered in [21].We also take into account the production of cosmogenic photons, both in the GeV-TeV band and at ultra-high energy, and analyse quantitatively the impact of multimessenger constraints on the UHECR source parameters.
While we were finalising this article, another study by Muzio et al. [22] on a subdominant population of UHE proton sources appeared.Unlike our work, they consider only mono-elemental injection by medium/heavy-composition cosmic ray sources and include insource photohadronic interactions (see also [20] for an earlier work using the same model).They assume a particular blackbody-like photon field within the source region resulting in photohadronic interactions of the cosmic rays, whereas we do not consider interactions in the source environment.Our results are therefore more general but the inferred parameters should be understood as effective parameters of the cosmic rays after they escape the source environment.In the case of astrophysical environments with small optical depth to hadronic and photohadronic interactions, such as low-luminosity gamma-ray bursts, BL Lac objects, and radio lobes of jetted active galactic nuclei, our results closely resemble the spectra produced inside the astrophysical sources.In addition, [22] assume the same redshift evolution of the emissivity for both UHECR source populations and only optimise the contributions of both source populations sequentially.Our results are complementary and more general in terms of the source parameters of both source populations which we optimise simultaneously.
UHE protons, should they exist, are of significant interest for "UHECR astronomy" due to their high rigidity and consequently weak deflections in magnetic fields.Additionally, if they are accelerated to energies beyond ∼ 10 19.7 eV, the cross-section for photo-pion production on CMB photons is enhanced due to the ∆-resonance.This effect, known as the Greisen-Zatsepin-Kuzmin (GZK) limit [23,24], leads to strong attenuation of UHE protons above this energy if they are produced in sources more distant than ∼ 100 Mpc (see e.g.[8]) and the abundant production of charged and neutral pions.The subsequent decay of these pions will result in a large flux of high-energy neutrinos and gamma rays.
In this paper we quantify the maximum flux of UHE protons compatible with current observations of UHECR spectrum and composition, considering multimessenger constraints from gamma rays and neutrinos.We investigate two separate scenarios for the maximum proton energy; (i) a high-E p max and (ii) a low-E p max scenario.A brief overview of the model is provided in Section 2. Injection and propagation of the cosmic rays are simulated with the Monte-Carlo framework CRPropa 3 [25,26], taking into account the interaction with the cosmic microwave background and extragalactic background light [27].The best-fit source parameters are obtained in Section 3 by comparing the model predictions with existing observations, and in Section 4 we discuss the expected multimessenger signal.A specific, exotic scenario with flux recovery beyond the GZK cutoff is presented in Section 5. Finally, we discuss our results in the context of similar existing studies in Section 6, and conclude in Section 7 that current UHECR data is compatible with a significant contribution by this additional proton component of up to 15% at 20 EeV.The precise value depends critically on the choice of the hadronic interaction model for air shower modelling and the maximum proton energy.

Methods
The primary, mixed-composition, UHECR sources (MIX) are modelled following the effective parametrisation introduced in [2] but with minor modifications detailed in [5].We assume the acceleration to be universal in particle rigidity1 , following a "Peters cycle", with a power-law source spectrum and an exponential cutoff at the highest energies.Sources within the MIX population are assumed as identical, with a volumetric emission rate for the five injected elements A ∈{ 1 H, 4 He, 14 N, 28 Si, 56 Fe}.Here Q E 0 A is the local (z = 0) emission rate at a normalisation energy E 0 ≪ E p max in erg −1 Mpc −3 yr −1 , and γ is the spectral index which is ≈ 2 for diffusive shock acceleration.The source emissivity, i.e. the luminosity density, can be derived from the emission rate as where we have chosen E min = 10 17.8 eV.
The predicted flux at Earth for an observed nuclear mass A ′ and energy E ′ , and for a redshift evolution n(z) of the source population emissivity 3) The last term translates the injected spectrum at the sources to the observed spectrum after propagation and is obtained via Monte-Carlo simulations with CRPropa.In general, the population-emissivity redshift evolution n(z) is composed of the evolution of per-source luminosities and the density evolution of the source population.In our analysis, we do not attempt to distinguish the difference between these effects and describe the evolution with a (broken) power law with z 0 = 1.5 and z max = 4 [28].Sources at z ≳ 1 have a negligible impact on the observed UHECR flux because of attenuation effects, but they play an important role for the expected multimessenger signal of co-produced neutrinos and low-energy gamma rays.A more conservative estimate of the cosmogenic neutrino flux is obtained if these high-redshift sources, which cannot be constrained by the cosmic-ray fit, are ignored.
For the additional population of UHE pure-proton sources (PP), we are particularly interested in the predicted flux of cosmogenic neutrinos at E ν ≈ 1 EeV since this corresponds to the peak sensitivity interval of many existing and planned neutrino experiments.If these neutrinos are produced in the interactions of cosmic rays with photon fields, they typically receive ∼ 5% of the primary CR energy [8], which implies that the relevant energy is E CR ≈ 20 EeV.We define this value as the reference energy at which we evaluate the contribution of the PP UHE protons to the observed flux of UHECRs.Properties of the pure-proton sources are described by the independent set of parameters E PP max , γ PP , m PP , and L PP 0 .The interactions of UHECRs with cosmic background photons lead to the production of secondary photons and neutrinos, with the strength of this "cosmogenic" multimessenger signal depending predominantly on the cosmic-ray composition, injection spectral index and source distance.We compare our model predictions for the UHECR spectrum and composition with publicly available data by Auger [29,30].Since the composition cannot be observed directly, the mean, ⟨X max ⟩, and standard deviation, σ(X max ), of the depth of the air-shower maximum are used as proxy observables, and the conversion is performed with the hadronic interaction models Epos-LHC [31] and Sibyll2.3c[32].To minimise the impact of a possible contribution of other sources dominating the observed flux below the ankle we limit our analysis to E CR ≥ 10 18.7 eV.However, spectral points at lower energies are included as upper limits and scenarios with excessive sub-ankle flux are rejected.
The best-fit source parameters are determined in a two-step fitting process.We discretise the parameter space in maximum energy/rigidity, spectral index and redshift evolution for both source classes and sample a large number of possible combinations of these parameters.For each of these possible source configurations, we then use the Levenberg-Marquardt algorithm2 to find the injection fractions f A of the MIX sources and emissivities [L 0 , L PP 0 ] of both source populations that minimise the χ 2 differences between our model predictions for the UHECR spectrum and composition and the Auger data points.If a reasonable fit (χ 2 < 250) is found for a particular combination of source parameters then adjacent points are also evaluated in an iterative process.This is different to the two-step approach of [22] since we optimise the contributions of both source populations at the same time.
Constraints on the source parameters derived from a comparison of the predicted cosmogenic flux of gamma rays and neutrinos with observations and upper limits are taken into account with additional ∆χ 2 -penalty terms.For observed fluxes, such as the Fermi-LAT IGRB [33] and parts of the IceCube HESE neutrino flux [34] we consider a simple one-sided χ 2 penalty that only contributes if the predicted flux exceeds observations.For upper-limit points with a low number of, or zero, events per bin, e.g. the Auger UHE neutrino [35] and UHE gamma-ray limits [36,37] we use the Poisson likelihood χ 2 [38] but the penalty is only applied if the predicted number of events in a bin exceeds the observed number.The relevant data sets are: ∆χ 2 ν : IceCube HESE flux [34] & Auger UHE neutrino limit [35], (2.5) ∆χ 2 γ : Fermi-LAT IGRB flux [33] , Auger hybrid UHE gamma-ray limit [36] & Auger SD UHE gamma-ray limit [37].
We exclude possible source configurations where the combined multimessenger penalty exceeds the level of two sigma, i.e. when ∆χ 2 ν +∆χ 2 γ > 4.However, in the plots of the cosmogenic neutrino and gamma-ray fluxes we only include the rejections by the respective messenger.

Fit with an Additional Proton Component
We investigate two different scenarios in terms of the proton maximum energy, assuming Epos-LHC as hadronic interaction model.Results for Sibyll2.3c are shown in Appendix A. In both scenarios, we find the redshift evolution of the PP number density to be unconstrained by cosmic-ray observations alone.Since the PP flux is pure protons, interactions during propagation do not affect the observed composition.However, propagation effects soften the distribution and attenuate the original UHE proton flux.Stronger redshift evolutions require harder injection spectra and higher source emissivity.
The two models described in the following, with maximum energy of the protons fixed to 10 5 EeV and 10 EeV respectively, represent the most characteristic scenarios identified during a scan of E PP max (see Appendix B).
3.1 Two-Source-Class Dip Model (2SC-dip) We are particularly interested in scenarios that produce a large flux of UHE neutrinos and gamma rays.This requires proton energies sufficiently above the GZK limit to enable copious photo-pion production on CMB photons, and we, therefore, choose E PP max = 10 23 eV for our first scenario.The best-fit properties of both source populations are listed in Table 1, 2nd column, and the predicted spectrum and composition at Earth are shown in Figure 1, left.The preferred maximum rigidity, spectral index and redshift evolution of the mixed-composition source population are compatible with the values obtained for the single-population model within uncertainties, and the additional protons provide a relatively constant contribution of approx.5 − 10% between the ankle and the end of the GZK cutoff (Figure 2, teal band).We find that the overall shape effectively corresponds to the predictions from the classical "proton-dip" explanation of the UHECR flux [39].While this model is inconsistent with current measurements of the UHECR composition and the high-energy neutrino flux [40], our results show that it can still be relevant if the total proton contribution remains subdominant to the primary, mixed-composition, cosmic-ray flux.We refer to the presented source model as the "dip" or 2SC-dip (two-source-class dip) model.
Here, the proton sources are required to exhibit a soft injection spectrum (see Appendix C), which could be a distinguishing feature of this additional source population in the observed flux, provided that reliable event-by-event mass reconstruction becomes available in the future.Softer spectra than suggested by the best fit are disfavoured since the associated sub-ankle flux would exceed observational limits.For hard spectra, γ PP ≲ 2, the additional protons only contribute at energies around the GZK cutoff and the possibilities for improving the fit over the entire energy range are consequently limited.The combination of both effects results in a clearly localised preferred spectral index of the proton sources.

Two-Source-Class Best-Fit Model (2SC-uhecr)
An alternative scenario is presented by proton sources with energies comparable to the standard, mixed-composition, cosmic-ray sources.For this model, we set E PP max = 10 EeV.At the best fit (Table 1, 3rd column), the improvement over the dip-model is ∆χ 2 ≈ −15 but very hard proton spectra are required (see Appendix C).The predicted PP proton spectrum at Earth exhibits a peak-like shape reminiscent of the individual, peaked, mass groups originating from the mixed composition sources (Figure 1, right).However, due to the choice of E PP max , the peak energy is shifted upward by approximately an order of magnitude compared to the mixed-population proton peak.Compared to the 2SC-dip model, the best-fit observed proton fraction at 20 EeV is significantly larger, up to 15%, but the contribution is limited to a small energy interval and becomes negligible below the ankle (Figure 2, brown band).
While this scenario, the "UHECR best-fit" model (2SC-uhecr), provides a significant improvement in the cosmic-ray fit, it comes at the cost of extremely hard proton injection spectra, and the expected cosmogenic neutrino and UHE gamma-ray signal associated with the protons is reduced due to the sub-GZK maximum proton energies.With the injection spectrum of the additional protons similar to the bulk of the cosmic rays, separation of the two components will be difficult even if event-by-event mass reconstruction were available.
Table 1.Best-fit parameters for the single-and two-population source models with EPOS-LHC used as the hadronic-interaction model describing air-shower development.The 1σ uncertainties include the penalty factor for the total best-fit quality proposed in [41].The "1SC" scenario is the benchmark model with only a single population of sources injecting mixed-composition cosmic rays."Population 1" refers to the baseline source class that injects a mixed cosmic-ray flux of protons to iron, and "Population 2" denotes pure-proton sources.The best fit of UHECR spectrum and composition is given in the "CR" column, and the best fit after including neutrino and gamma-ray limits in the "CR + MM" columns.For the 2SC-uhecr model, the cosmic-ray best fit is compatible with existing multimessenger limits.Confidence intervals that extend to the edges of the sampled parameter range are indicated by an asterisk.3.77 +0.06

Model
68.5 +3.9However, the predicted existence of two separate proton bumps in the cosmic-ray spectrum is a distinguishing feature of this model.

Multimessenger Signal
In the following, we discuss the predicted multimessenger signal produced through interactions with the CMB and the Extragalactic Background Light during the propagation of the cosmic rays.We focus on the 2SC-dip "proton-dip" model which predicts a large flux of cosmogenic neutrinos and UHE gamma rays.The multimessenger signal of the 2SC-uhecr model is briefly discussed at the end.

2SC-dip
Photons, electrons, and positrons produced with PeV-EeV energies in photohadronic interactions of the UHE protons interact with cosmic photon fields, leading to the development of electromagnetic cascades and reprocessing to lower energies.In the scenario of low-E max and mixed-composition cosmic-ray sources only, most of the gamma-ray signal is expected at GeV-PeV energies since the CR energies are insufficient for large interaction cross-sections with CMB photons.In this energy range (Figure 3, left), the predicted gamma-ray flux associated with the PP protons in our model is at a similar level to the flux expected from the mixed cosmic rays.Depending on the exact choice of source parameters, the combined gamma-ray flux of both populations can saturate the upper limit imposed by the re-scaled 3Fermi-LAT flux at ∼ 700 GeV, however, the tension is not statistically significant.Most of the gamma-ray flux at E γ ≳ 100 GeV is produced by the mixed-composition cosmic rays.At lower energies, the cosmogenic gamma rays are safely below the observed diffuse background flux.
The situation is more promising at ultra-high energies where the signal from the ordinary, mixed cosmic rays is expected to be very small.By construction, the protons injected at the PP sources have typical energies E γ > 10 18 eV and consequently large cross-sections for photo-pion production on the abundant CMB photons.The predicted UHE gamma-ray flux from the protons is therefore orders of magnitude above the flux produced by the mixed cosmic rays (Figure 3, right).It correlates inversely with the PP spectral index -harder injection spectra result in more cosmogenic UHE photons.As indicated previously, hard injection spectra generally require strongly positive redshift evolutions to soften the observed spectrum.Present limits by Auger and TA are not constraining, even in the most optimistic scenario within 3σ uncertainties, however, the difference is not more than a factor of a few and it is clear that future detectors -such as GRAND200k and AugerPrime -will provide strong constraints for the viable PP spectral index and redshift evolution.
The expected flux of cosmogenic neutrinos (Figure 4) is not well constrained by the  1, 2nd column) is indicated by a solid line.The 1, 2, 3σ contours, under the condition that ∆χ 2 γ < 4, are indicated by brown bands in decreasing intensity for the contribution from the additional protonic UHECRs, and by blue bands for the gamma-ray flux from the regular, mixed cosmic rays.These intervals do not include the best-fit penalty factor of [41].Observations include the Fermi-LAT [33] and HAWC [43] diffuse gamma-ray background in the GeV-PeV range, the 95% upper limits at UHE of Auger [36,37] and TA [44], the optimistic 3-year sensitivity of the planned GRAND200k [45], and a combination of the latest Auger SD limit with the projected AugerPrime exposure for 10 years of observations under the assumption of 100% photon selection efficiency and zero background.cosmic-ray fit alone and can vary by approx.a factor of 1000 within the 99.7% confidence interval.In the most pessimistic case, when the redshift evolution of the proton sources is strongly negative, the neutrino flux produced by PP protons is subdominant to the neutrinos from the default CR population at all energies E ν ≲ 1 EeV and the UHE flux is small.On the other hand, for strong redshift evolutions, the expected neutrino flux saturates the flux observed by IceCube in the few-PeV energy range and exceeds significantly the limits above 10 PeV and at UHE.This includes the source configuration corresponding to the best UHECR spectrum and composition fit.By requiring that the neutrino limits are not violated (∆χ 2 ν < 4) we can constrain the properties of the proton sources to γ PP ≳ 1.6 , m PP ≲ 4 , and L PP 0 ≲ 10 44.5 erg Mpc 3 yr .
Irrespective of the total level, the predicted neutrino flux exhibits a characteristic doublebump profile, with the first peak at E ν ≈ 5 PeV from photo-pion production of the cosmicray protons on the extragalactic background light, and the second peak at E ν ≈ 1 EeV from photo-pion production on the less energetic, but more abundant, CMB photons.Due to the soft spectrum of the UHE protons, both peaks are present at the same time and the UHE neutrino limits can be used to constrain the contribution of this cosmogenic neutrino flux to the observed IceCube HESE flux at 1.3 PeV to f PP HESE ≲ 20%.

2SC-uhecr
In the "UHECR best fit" model, the maximum proton energy is below the required level for photo-pion production with the bulk of CMB photons and the expected multimessenger signal is low.UHE gamma rays are at least three orders of magnitude below existing limits  1, 2nd column).The shaded confidence intervals include the additional χ 2 penalty from the existing neutrino limits.The IceCube HESE flux [34], upper limits from IceCube [46] and Auger [35,47], and predicted sensitivities of planned detectors [45,[48][49][50] are shown as a reference.
and at GeV-TeV energies, the contribution is subdominant compared to the cosmogenic photons from the MIX cosmic rays.The total contribution to the Fermi-LAT IGRB is < 50% even in the most optimistic scenario, although the upper limit in the highest energy bin is approximately saturated.While the neutrino signal of the UHE protons at the best fit is subdominant to the neutrinos from the mixed-composition cosmic rays, the shape of the neutrino spectrum is of particular interest.Unlike for the 2SC-dip model, few protons are present at lower energies and the low-energy peak originating from interactions with EBL photons is therefore absent.Only the peak from photo-pion production on the CMB remains.In this scenario, the observed IceCube neutrino flux at PeV energies and below, and the possible UHE neutrino flux are decoupled.It is possible, for strongly positive redshift evolutions of the proton sources, to produce a large neutrino flux at UHE with a negligible contribution to the IceCube HESE flux.Redshift evolutions stronger than m PP ≈ 4 can be excluded by the current UHE neutrino limits of IceCube and Auger.

Exotic Flux Recovery Scenario (2SC-rec)
A combination of the 2SC-dip and 2SC-uhecr models is provided by a proton source population with large maximum energy, E PP max = 10 5 EeV, as in the "proton-dip" model, and hard injection spectrum, γ PP = 1, similar to the "UHECR" model.The quality of the UHECR fit (Table 2) is reduced compared to the other two models and approaches the baseline singlesource-class model.Compared to the dip model, the potential for fit improvement is limited since the protons contribute only at the highest energies, while the position of the observed proton peak is at too-high energies to provide an improvement of similar magnitude as in the 2SC-uhecr model.However, an interesting feature in the form of a "flux recovery" at trans-GZK energies can be observed (Figure 5).We refer to this third model as the "recovery" or 2SC-rec model.
A recovery is only possible if the nearest source(s) is(are) located within the GZK volume at no more than ∼ 20 Mpc [e.g.8] as otherwise, the GZK cutoff provides a natural suppression  1 but for the "flux recovery" 2SC-rec scenario.Auger 90% upper limits above 10 20.4 eV were derived assuming an energy-independent exposure of 60 400 km 2 yr sr [29].Expected 90% upper limits for GCOS (40k) after 10 years of operation (ϵ ∼ 10 6 km 2 yr sr [51]) are shown in purple.
of the observable flux above ∼ 10 19.7 eV.Such a spectral recovery is not necessarily connected to a large UHE neutrino signal.In addition to a high E PP max and hard proton source spectra, the latter also requires strong redshift evolution of the source emissivity which is not a prerequisite for a CR flux recovery.However, the observation of neutrinos with energy above 10 19 eV by future extremely-UHE neutrino detectors such as PUEO [52] would provide a strong hint for the existence of a sizeable UHECR flux recovery beyond the GZK cutoff.
We have noted in Section 4 for large E PP max that hard proton source spectra are excluded by existing neutrino limits, under the condition that only source configurations within 3σ of the best fit to the UHECR spectrum and composition under the 2SC-dip model are considered.If this limitation is lifted, such as for the 2SC-rec model, we can identify scenarios where the predicted neutrino flux is sufficiently below existing limits, i.e. ∆χ 2 ν < 4.This constrains the redshift evolution of the proton sources to m PP ≲ 3 (≲ 2 if gamma-ray limits are included).
In contrast to the 2SC-dip model, the combination of hard injection spectrum, high maximum proton energy and uniform source distribution with minimum distance z min = 10 −3 , results in an increased flux of cosmogenic UHE gamma rays.We find that for spectral indices harder than γ PP ≲ 1 all possible realisations of the source model are excluded by the existing Auger UHE photon limits [36,37].We conclude that the joint consideration of neutrino and UHE gamma-ray limits severely constrains the allowed proton injection spectrum and, by extension, the maximum allowed flux recovery from this second source population above the GZK cutoff.This motivates our choice of γ PP = 1 as benchmark spectral index for the 2SC-rec model.The spectrum and composition corresponding to the maximum recovery allowed by the cosmic-ray fit and multimessenger constraints are shown in Figure 5.The preferred source parameters stay unchanged except for the PP redshift evolution and luminosity density.Finally, we comment that the projected sensitivity of the proposed Global Cosmic Ray Observatory (GCOS) [51] would place strong constraints on the  1 but for the extreme 2SC-recovery model.The best-fit parameters of the mixed-composition sources are given for the CR-only fit, but the preferred values for the CR+MM scenario are compatible within quoted uncertainties.Our results for this scenario can look similar to what is expected from propagation models involving Lorentz invariance violation (LIV) (compare with [53][54][55], for example), which can suppress the photopion interaction at high energies and strengthen the flux recovery.If a recovery is observed, it is therefore prudent to investigate whether this is due to LIV or a 2SC-dip/recovery scenario.Clear differences between LIV models and a 2SC-dip/recovery scenario are the expected arrival directions as well as the expected cosmogenic neutrino and photon fluxes.

Discussion and Comparison to Previous Studies
Here we discuss the implications of fitted parameters and compare our findings to past works investigating a subdominant proton flux at the highest energies.
A similar conclusion in terms of the allowed UHE proton fraction for Epos-LHC versus Sibyll2.3c was reached in the recent paper by Muzio et al. [22].The best fit obtained in this study (not explicitly discussed in their paper) is qualitatively similar to our 2SCuhecr model with hard proton injection spectrum and low maximum energies4 (see also [20]).Their reported, best-fit observed proton fraction of the integral cosmic-ray flux above 30 EeV for typical astrophysical source evolutions is 5 − 10% and 2 − 3% for Epos-LHC and Sibyll2.3crespectively.These values are compatible with our preferred integral fractions, F p (≥ 30 EeV) = 10.2 +1.4 −1.5 % / 3.1 +2.6 −0.9 %.A direct comparison is difficult, however, since the authors assumed a mono-elemental injection of Silicon-like nuclei at the MIX sources and also included in-source photohadronic interactions, and the predicted cosmic-ray flux at Earth is not provided.
A solution to the two-population model with source parameters similar to our 2SC-dip best fit was found by Das et al. [21] who report a best-fit proton fraction of the observed flux in the highest energy bin of approx. 1 − 3%.This corresponds to 20 − 25% at E ref = 20 EeVa contribution that we found to be in strong tension with the observed UHECR composition, in particular the variance of shower maxima.They assumed a maximum source distance of z max = 1 which results in a conservative estimate of the associated flux of cosmogenic neutrinos.We consider, instead, sources out to redshift z max = 4, based on the approximate redshift evolution of source emissivities for probable astrophysical sources of UHECRs, resulting in a neutrino flux exceeding their prediction by more than an order of magnitude at E ν = 1 EeV.This enables us to use existing upper limits on the UHE neutrino flux by Auger and IceCube to significantly constrain the redshift evolution of the PP source population emissivity.Unlike the present work, the expected flux of cosmogenic gamma-rays is not discussed extensively in [21].The authors do not recover the best fit of the UHECR spectrum and composition that we identify in our 2SC-uhecr model since they only consider soft spectral indices of the proton sources, γ PP ≥ 2.2, and no mention is made of a possible flux recovery beyond the GZK cutoff.
A recent study by the Auger Collaboration [19] also investigated the co-existence of several source populations to explain the entire UHECR spectrum and composition above and below the ankle simultaneously.While the focus of that paper is somewhat different, their scenario 1 resembles our proton-bump model and the best-fit parameters are generally compatible.Mild disagreement can be identified for the preferred spectral index of the proton sources, which they predict to be much softer, and the redshift evolution of the mixed-composition sources, which they predict to be substantially stronger.Both of these are likely related the lower limit on the energy range in [19].They require a larger proton flux below the ankle to explain the entire observed flux while we only use those data points as upper limits.However, their result depends on the assumptions for the sub-ankle nitrogen flux component which they have fixed ad-hoc.
Important information about the potential sources of the UHE pure-proton flux can be gained from the total emissivity L PP 0 (luminosity density) required by the UHECR fit.Although the cosmic-ray emissivity of astrophysical objects is generally not known, other observable properties such as gamma-ray and X-ray emissivities can be used for relative calibration.For a summary of population emissivities see [56].Assuming equipartition of the available energy budget into gamma rays / X-rays and cosmic rays, we observe that all typically considered source classes (gamma ray bursts, tidal disruption events, starburst galaxies, active galactic nuclei, BL Lacertae, flat-spectrum radio quasars, and radio galaxies) can satisfy the emissivity required of the pure-proton sources in the 2SC-uhecr and 2SC-dip models, although gamma-ray bursts and tidal disruption events are marginally challenged in the latter scenario.For the extreme 2SC-recovery model, only the entire AGN population and the population of all BL Lacs can easily meet the required emissivity.GRBs and TDEs, in contrast, are excluded unless their cosmic-ray emissivity exceeds the observed gamma-ray emissivity by at least a factor of ten.FSRQs and radio galaxies sit close to the minimum luminosity density required by the cosmic-ray fit.
Given the hard spectrum and high maximum energy, it might be challenging that the UHE proton flux predicted by the 2SC-rec model is produced by astrophysical accelerators.An alternative explanation for the spectrum could be provided by the decay of hypothetical super-heavy dark matter (SHDM) with masses up to the Planck mass [57][58][59][60][61][62][63][64].These heavy particles can be produced gravitationally during the early stages of the Universe, e.g. as part of the reheating epoch from a hypothesised, decaying inflaton field, or from coherent oscillations of this field before the inflation phase [65][66][67].If they never reached thermal equilibrium after production and the lifetime is larger than the age of the Universe then these heavy relics can provide a possible explanation for observed DM densities [63].Similar to the original proton-dip model [39], "top-down" scenarios of decaying SHDM are disfavoured as the single origin of the observed UHECR flux [68,69], and it was shown that decaying SHDM cannot explain the detected high-energy IceCube neutrino events if a hadronic decay channel is considered [70][71][72].Still, a subdominant contribution to the observed UHECR flux, and a possible flux recovery due to very hard decay spectra are not fully excluded.Crucially, existing upper limits on the post-GZK cosmic-ray flux provide only weak constraints on the allowed flux recovery, and UHE photon limits prove superior for M DM < 10 14 GeV [64].
We do not investigate the SHDM scenario further, however, we wish to point out several key differences compared to our assumed source model.If the additional protons are produced in the decay of super-heavy dark matter, a substantial anisotropy in arrival directions and extremely local production of the observed UHECRs should be expected since the signal is predicted to be dominated by dark matter in the Milky Way with a particular clustering around the Galactic centre [73].This is in sharp contrast to our proposed continuous distribution of sources in redshift up to z max = 4 and minimum source distance of ∼ 4 Mpc.Consequently, in the SHDM scenario, the expected flux of cosmogenic neutrinos and low-energy gamma rays is severely reduced.In addition, we only consider the cosmogenic production of neutrinos and gamma rays while in the SHDM model the multimessenger signal is likely dominated by production during the decay of the dark matter.

Summary and Conclusions
In this work, we have investigated the possible existence, and allowed parameter space, for an additional, proton-dominated component of UHECRs, produced by an independent astrophysical source population.We have presented the maximum contribution of such a population to the UHECR flux at Earth, taking into account the fit to the UHECR spectrum and composition-sensitive observables.In addition, we have derived predictions for the spectral shape and redshift evolution of the independent UHE-proton population model as well as the expected secondary neutrino and photon fluxes produced by UHECR interactions and their detectability.
This analysis was performed for two distinct choices of the maximum proton energy.For sources with maximum energy far beyond the GZK limit (2SC-dip model), the proton spectrum at Earth reproduces the predictions of the classic "proton-dip" model [39], albeit with the proton flux subdominant to the contribution of the principal, mixed-composition cosmic rays.If instead maximum energies below 10 19.7 eV are assumed (2SC-uhecr model), the cosmic-ray fit is improved by ∆χ 2 ≈ −15 but the source spectrum must be hard and the associated multimessenger signature is generally small.In both scenarios, the redshift evolution of the proton sources cannot be constrained by the cosmic-ray fit alone.
We find that the maximum proton contribution to the observed, diffuse UHECR flux depends strongly on the choice of hadronic interaction model for the interpretation of the extensive air showers, and on the maximum proton energy.With Sibyll2.3c a proton fraction of ≲ 1% is expected at 20 EeV in the 2SC-dip model and the improvement over the baseline model is negligible.Under the 2SC-uhecr model, a contribution of 2 − 5% is predicted with a minor 1.1σ significance compared to the baseline one-population model.Assuming Epos-LHC instead, for the 2SC-dip model, approximately 8% of the UHECR flux is expected to be protons, with the contribution nearly constant over the entire energy range above the ankle.For the 2SC-uhecr model, where E p max = 10 EV, the contribution to the observed UHECR flux peaks around E ref ≈ 20 EeV at up to 15%, but the relative proton fraction decreases rapidly for energies away from the peak and the source spectra are required to be hard.The improvement of the two-population model over the baseline single-population scenario is 2.2σ (2SC-dip) and 3.7σ (2SC-uhecr).
We demonstrated that for our fiducial high-E p max model a distinguishing feature of the independent UHE proton component is a soft spectral index (γ = 2.5 ± 0.3), which can be tested by AugerPrime or other facilities with event-by-event mass determination capabilities.In addition, the cosmogenic neutrino and UHE photon fluxes produced by this component are substantial and dominate over those from the mixed-composition population.Current neutrino upper limits from IceCube and Auger already weakly constrain the available parameter space for the proton population from the fit to the UHECR data alone.
Finally, as an "exotic" scenario, we have considered proton sources with high maximum energy E p max ≫ 10 20 eV and hard spectral index.We find that existing limits on the neutrino and UHE gamma-ray flux constrain the proton spectral index to γ p ≳ 1 and therefore provide an upper limit on the possible cosmic ray flux beyond the GZK cutoff.However, a significant recovery is still allowed.Values on the order of 10 EeV are preferred (corresponding to the 2SC-uhecr model); however, trans-GZK maximum energies cannot be rejected at appreciable significance (including the penalty for the quality of the global best fit [41]).

B Proton Maximum Energy
Adding another free parameter in the form of the maximum proton energy E PP max to our fit is computationally prohibitive for our regular parameter pixelisation used in Section 3.However, if the resolution is reduced then a scan over E PP max is possible.The results are shown in Figure 8.We note that E PP max ≲ 3 EeV can be rejected with more than 4σ confidence.At this energy, the observed proton flux produced by the second population becomes coincident with the protons from the default sources (primary or from disintegration).We thus obtain no improvement of the fit compared to the single-population scenario.
In addition, the fit is asymptotically insensitive to the maximum proton energy for E PP max ≳ E GZK .This justifies our choice of 10 23 eV for the 2SC-dip scenario as a representative case for extremely-UHE proton sources.

C Proton Source Parameter Space
Results of the source-parameter scan are shown in Figure 9 as a 2D surface plot over the PP injected spectral index γ PP src and redshift evolution m PP for both the 2SC-dip and 2SC-uhecr models.

Figure 1 .
Figure 1.Predicted spectrum and composition at Earth for the investigated scenarios, with Epos-LHC as hadronic interaction model.Left: "proton-dip" (2SC-dip).Right: "UHECR" best fit (2SCuhecr).Best-fit parameter values are listed in Table 1.Dashed lines indicate the contributions of the separate mass groups from the mixed-composition sources, with [A min , A max ].The additional protons from the second population are shown as a solid, orange line.Coloured bands indicate the 68% uncertainties.

Figure 2 .
Figure 2. Contribution of the PP protons to the observed, differential UHECR flux as a function of energy, within 1σ of the best fit to CR spectrum and composition (see Figure 1).

Figure 3 .
Figure 3. Predicted cosmogenic gamma-ray signal for the "proton-dip" model (2SC-dip), with Epos-LHC as hadronic interaction model, in the GeV-TeV (left) and EeV (right) energy range.The photon flux for each source class corresponding to the UHECR best fit (Table1, 2nd column) is indicated by a solid line.The 1, 2, 3σ contours, under the condition that ∆χ 2 γ < 4, are indicated by brown bands in decreasing intensity for the contribution from the additional protonic UHECRs, and by blue bands for the gamma-ray flux from the regular, mixed cosmic rays.These intervals do not include the best-fit penalty factor of[41].Observations include the Fermi-LAT[33] and HAWC[43] diffuse gamma-ray background in the GeV-PeV range, the 95% upper limits at UHE of Auger[36,37] and TA[44], the optimistic 3-year sensitivity of the planned GRAND200k[45], and a combination of the latest Auger SD limit with the projected AugerPrime exposure for 10 years of observations under the assumption of 100% photon selection efficiency and zero background.

Figure 4 .
Figure 4. Same as Figure 3 but for the predicted cosmogenic neutrinos in the 2SC-dip (left) and 2SCuhecr model (right).The maximum allowed flux within 3σ of the best CR fit but without including the multimessenger penalty is shown as a dashed line of the respective colour.Solid lines indicate the neutrino flux corresponding to the cosmic-ray best fit without multimessenger constraints (Table1, 2nd column).The shaded confidence intervals include the additional χ 2 penalty from the existing neutrino limits.The IceCube HESE flux[34], upper limits from IceCube[46] and Auger[35,47], and predicted sensitivities of planned detectors[45,[48][49][50] are shown as a reference.

Figure 7 .
Figure 7. Same as Figures 3 and 4 for the "proton-dip" model (2SC-dip) but with Sibyll2.3c as hadronic interaction model.Left: UHE gamma rays.Right: neutrinos.The jagged upper limit of the UHE gamma-ray flux is a result of limited statistics in the numerical simulation.

Figure 8 .
Figure 8. Fit quality (∆χ 2 ) compared to the best fit as function of the proton maximum energy E PP max .Values on the order of 10 EeV are preferred (corresponding to the 2SC-uhecr model); however, trans-GZK maximum energies cannot be rejected at appreciable significance (including the penalty for the quality of the global best fit[41]).

Figure 9 .
Figure 9. Fit quality for the 2SC-dip (left) and 2SC-uhecr (right) model, marginalised onto γ PP src ×m PP space.The best fit is marked with a white cross and contour lines indicate the one (green), two (orange) and three (red) sigma confidence intervals.

Table 2 .
Same as Table

Table 3 .
Same as Table1but for Sibyll2.3c as hadronic interaction model.The 2SC-dip best fit is excluded by the neutrino limits at ∆χ 2 ν ≈ 4, however, compatibility is obtained for m PP = 6 → 5.For the 2SC-uhecr scenario, the CR best fit is again compatible with the multimessenger constraints.