Modeling the spectrum and composition of ultrahigh-energy cosmic rays with two populations of extragalactic sources

We fit the ultrahigh-energy cosmic-ray (UHECR, E≳0.1\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$E\gtrsim 0.1$$\end{document} EeV) spectrum and composition data from the Pierre Auger Observatory at energies E≳5·1018\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$E\gtrsim 5\cdot 10^{18}$$\end{document} eV, i.e., beyond the ankle using two populations of astrophysical sources. One population, accelerating dominantly protons (1\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$^1$$\end{document}H), extends up to the highest observed energies with maximum energy close to the GZK cutoff and injection spectral index near the Fermi acceleration model; while another population accelerates light-to-heavy nuclei (4\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$^4$$\end{document}He, 14\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$^{14}$$\end{document}N, 28\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$^{28}$$\end{document}Si, 56\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$^{56}$$\end{document}Fe) with a relatively low rigidity cutoff and hard injection spectrum. A significant improvement in the combined fit is noted as we go from a one-population to two-population model. For the latter, we constrain the maximum allowed proton fraction at the highest-energy bin within 3.5σ\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\sigma $$\end{document} statistical significance. In the single-population model, low-luminosity gamma-ray bursts turn out to match the best-fit evolution parameter. In the two-population model, the active galactic nuclei is consistent with the best-fit redshift evolution parameter of the pure proton-emitting sources, while the tidal disruption events could be responsible for emitting heavier nuclei. We also compute expected cosmogenic neutrino flux in such a hybrid source population scenario and discuss possibilities to detect these neutrinos by upcoming detectors to shed light on the sources of UHECRs.


Introduction
Identifying the sources of ultrahigh-energy cosmic rays (UHECRs, E 0.1 EeV) is one of the outstanding problems in astroparticle physics [1,2]. Active Galactic Nuclei (AGNs) residing at the centers of nearby radio-galaxies are a e-mail: saikatdas@rri.res.in (corresponding author) b e-mail: srazzaque@uj.ac.za c e-mail: nayan@rri.res.in considered to be a potential candidate source class of UHECR acceleration [3][4][5][6][7]. Studies involving the origin of TeV γ -rays assert blazars as ideal cosmic accelerators [8][9][10][11]. A recent analysis by the Pierre Auger Observatory has found a possible correlation between starburst galaxies and the observed intermediate scale anisotropy in UHECR arrival directions, with a statistical significance of 4σ in contrast to isotropy [12][13][14]. There are also propositions of other transient highenergy phenomena like gamma-ray bursts (GRBs) [15][16][17][18][19][20], tidal disruption events (TDEs) of white dwarfs or neutron stars [21][22][23][24], as well as, pulsar winds [25,26] which can reach the energy and flux required to explain the observed UHECR spectrum. Nevertheless, a direct correlation of these known source catalogs, derived from X-ray and γ -ray observations, with an observed UHECR event is yet to be made [27][28][29][30]. The different source classes allow an extensively wide range of UHECR parameters to be viable in the acceleration region. UHECRs produce neutrinos and γ -rays on interactions with the cosmic background photons during their propagation over cosmological distances. The current multimessenger data can only constrain UHECR source models and provide hints towards plausible accelerator environments [31,32], rejecting the possibility of a pure proton composition at the highest energies [33][34][35][36][37]. Deflections in Galactic and extragalactic magnetic fields pose an additional challenge in UHECR source identification.
The Pierre Auger Observatory (PAO) in Malargüe, Argentina [38] and the Telescope Array (TA) experiment in Utah, United States [39] are attaining unprecedented precision in the measurement of UHECR flux, composition, and arrival directions from 0.3 EeV to beyond 100 EeV using their hybrid detection technique [40,41]. On incidence at the Earth's atmosphere, these energetic UHECR nuclei initiate hadronic cascades which are intercepted by the surface detector (SD), and the simultaneous fluorescence light emitted by the Nitrogen molecules in the atmosphere is observed using the fluorescence detector (FD). This extensive air shower (EAS) triggered by the UHECRs is recorded to measure the maximum shower-depth distribution (X max ) [42]. However, even with the large event statistics observed by PAO, the mass composition is not as well constrained as the spectrum and anisotropy up to ∼ 100 EeV [43]. The first two moments of X max , viz., the mean X max , and its fluctuation from showerto-shower σ (X max ) serves the purpose of deducing the mass composition. The standard shower propagation codes, eg., corsika [44], conex [45], etc., depend on the choice of a hadronic interaction model and photodisintegration crosssection, which are extrapolations of the hadronic physics to the ultrahigh-energy regime. Uncertainties in these models propagate to uncertainties in the reconstruction of the masscomposition of observed events. Lifting the degeneracy in the mass composition will be essential to constrain the source models.
The current LHC-tuned hadronic interaction models viz., sybill2.3c [46], epos-lhc [47], and qgsjet-II.04 [48] differ in their inherent assumptions and thus lead to different inferences of the mass composition using the same observed data. Current estimates from PAO predict that the relative fraction of protons decreases with increasing energy above 10 18.3 eV for all three models. For the first two models, N dominates at 10 19.6 eV, while for the third model, the entire contribution at the highest energy comes from He. The ankle at E ≈ 10 18.7 eV corresponds to a mixed composition with He dominance and lesser contributions from N and H, except for qgsjet-II.04 which suggests a zero N fraction [49]. The ankle is often inferred as a transition between two or more different populations of sources, leading to a tension between the preference of Galactic or extragalactic nature of the subankle spectrum. Based on the observed anisotropy and light composition, some UHECR models invoke increased photohadronic interactions of UHECRs in the environment surrounding the source. The magnetic field of the surrounding environment can confine the heavier nuclei with energies higher than that corresponding to the ankle, while they undergo photo-disintegration/spallation to produce the light component in the sub-ankle region [50][51][52]. This requires only a single class of UHECR sources that accelerate protons and nuclei. However, it is also possible to add a distinct light nuclei population of extragalactic origin that can explain the origin of the sub-ankle spectrum [53][54][55]. A purely protonic component, in addition to a Milky Way-like nuclear composition, has also been studied [55]. The proton fraction in the UHECR spectrum for various source models can be constrained through composition studies and compliance to multimessenger data [56,57].
In this work, first, we perform a combined fit of spectrum and composition data at E 5 · 10 18 eV measured by PAO [58], to find the best-fit parameters for a single-population of extragalactic UHECR sources injecting a mixed compo-sition of representative elements ( 1 H, 4 He, 14 N, 28 Si, 56 Fe). The best-fit 1 H abundance fraction is found to be zero in this case, conceivable within our choice of the photon background model, photodisintegration cross-section, and hadronic interaction model. Next, we show that within the permissible limit of current multimessenger photon and neutrino flux upper limits [59,60], the addition of a purely protonic ( 1 H) component up to the highest-energy bin can significantly improve the combined fit of spectrum and composition. We consider this component originates from a separate source population than the one accelerating light-to-heavy nuclei and fit the region of the spectrum above the ankle, i.e., E 5 · 10 18 eV. The best-fit values of the UHECR parameters are calculated for both the populations, allowing for a one-to-one comparison with the single-population case. We study the effect of variation of the proton injection spectral index, which is not done in earlier studies and indicate the maximum allowed proton fraction at the highest-energy bin up to 3.5σ statistical significance. We calculate the fluxes of cosmogenic neutrinos that can be produced by these two populations. We also explore the prospects of their observation by upcoming detectors, and probe the proton fraction at the highest-energy of the UHECR spectrum. Lastly, we take into account the redshift evolution of the two source populations, which is found to further improve the combined fit. We interpret the credibility of the best-fit redshift distributions in light of known candidate classes. We explain our model assumptions and simulation setup in Sect. 2 and present our results for both single-population and two-population models in Sect. 3. We discuss our results and possible source classes in light of the two-population model in Sect. 4 and draw our conclusions in Sect. 5.

UHECR propagation and shower depth distribution
UHECRs propagate over cosmological distances undergoing a variety of photohadronic interactions. These interactions lead to the production of secondary particles, viz., cosmogenic neutrinos and photons. The dominant photopion production of UHECR protons on the cosmic microwave background (CMB) via delta resonance occurs at ≈ 6.8 × 10 19 eV, producing neutral and charged pions (π 0 , π + ) with 2/3 and 1/3 probability, respectively. The neutral pions decay to produce γ -rays (π 0 → γ γ ), while the charged pions decay to produce neutrinos (π + → μ + +ν μ → e + +ν e +ν μ +ν μ ). Neutrinos can also be produced through other pγ processes and neutron beta decay (n → p + e − + ν e ). Bethe-Heitler interaction of UHECR protons of energy ≈ 4.8 × 10 17 eV with CMB photons can produce e + e − pairs. The e + and e − produced through various channels can iteratively produce high-energy photons by inverse-Compton scattering of cosmic background photons or synchrotron radiation in the extragalactic magnetic field (EGMF). The produced photons can undergo Breit-Wheeler pair production. All these interactions also hold for heavier nuclei ( A Z X , Z > 1), in addition to photodisintegration. The interactions may also occur with the extragalactic background light (EBL), having energy higher than the CMB, with cosmic-rays of lower energy. Besides, all particles lose energy due to the adiabatic expansion of the universe. We consider ΛCDM cosmology with the parameter values H 0 = 67.3 km s −1 Mpc −1 , [61]. While cosmic rays are deflected by the Galactic and extragalactic magnetic fields, the neutrinos travel unaffected by matter or radiation fields, and undeflected by magnetic fields.
The observed spectrum depends heavily on the choice of injection spectrum. We consider all elements are injected by the source following the spectrum given by, This represents an exponential cutoff power-law function, where K i and α are the abundance fraction of elements and spectral index at injection. A 0 and E 0 are arbitrary normalization flux and reference energy, respectively. A similar spectrum has been considered in the combined fit analysis by the PAO [43]. The broken exponential cutoff function is written as, We use the CRPropa 3 simulation framework to find the particle yields obtained at Earth after propagating over extragalactic space from the source to the observer [62]. We find the best-fit values of the UHECR parameters α, rigidity cutoff (R cut ) and K i for both one-population and two-population models. The normalization depends on the source model and the source population. The spectrum of EBL photons and its evolution with redshift is not as well known as for CMB. We use a latest and updated EBL model by Gilmore et al. [63] and talys 1.8 photodisintegration cross-section [64].
We use the parametrizations given by PAO based on the Heitler model of EAS to calculate the mean depth of cosmicray air shower maximum X max and its dispersion from the first two moments of ln A [65,66].
where X max p is the mean maximum depth of proton showers and f E is a parameter which depends on the energy of the UHECR event, where ξ , D, and δ depend on the specific hadronic interaction model. σ 2 ln A is the variance of ln A distribution and σ 2 sh is the average variance of X max weighted according to the ln A distribution, where σ 2 p is the X max variance for proton showers depending on energy and three model-dependent parameters. In this work, we use the updated parameter values 1 obtained from the conex simulations [45], for one of the post-LHC hadronic interaction models, sybill2.3c.

Results
We perform a combined fit of our UHECR source models to the spectrum and composition data measured by PAO [49,67], for one-population and two-population model of the UHECR sources. The fit region corresponds to energies above the ankle, i.e., E 5 · 10 18 eV in the spectrum, as well as, composition. We calculate the goodness-of-fit using the standard χ 2 formalism, where the subscript j corresponds to any of the three observables, viz., spectrum, X max , or σ (X max ). To find the best-fit cases, we minimize the sum of all the χ 2 j values. Here y obs i (E) is the measured value of an observable in the i−th energy bin corresponding to a mean energy E and y mod i (E; a M ) is the value obtained numerically. a M are the best-fit values of M parameters varied in the simulations. σ i are the errors provided by PAO. We denote the spectral fit as χ 2 spec and the composition fit as χ 2 comp . The latter represents the goodnessof-fit considering X max and σ (X max ) simultaneously. In the following subsections, we demonstrate the one-population model in Sect. 3.1, the transition due to the addition of an exclusive proton injecting class in Sect. 3.2, and finally the effects of redshift distribution in Sect. 3.3.

One-population model
We start by considering a single population of extragalactic sources up to a redshift z = 1, injecting a mixed composition of representative elements 1 H, 4 He, 14 N, 28 Si, and 56 Fe following an injection spectrum given by Eq. (1). The elements are injected with energy between 0.1 and 1000 EeV. The combined fit analysis done by PAO argues that only particles originating from z 0.5 are able to reach Earth with E > 5·10 18 eV [43,68,69]. Indeed, in our case, the contribution at the spectral cutoff comes from 56 Fe. Hence, the sources which are located further in the distance than z max = 1 are unable to contribute to the spectrum above the ankle (≈ 10 18.7 eV) [see, eg., Appendix C of 31]. This is because, as the distance of such heavy nuclei injecting sources increases, the rate of photodisintegration also gradually increases, thus decreasing their survival rate at the highest energies. Moreover, it was found that increasing z max has no effect on the best-fit parameters found with z max = 1 [32]. The source distribution is assumed to be uniform over comoving distance. In a later subsection, we check the effects of a non-trivial redshift evolution for one-population model. We scan the parameter space by varying the rigidity cutoff log 10 (R cut /V) between [18.0, 18.5] with a grid spacing of 0.1 and the injection spectral index α between [− 1.5, 1.0] with a grid spacing of 0.1. For each set of values {α, log 10 (R cut /V)}, we find the best-fit abundance fraction of the injected elements. The number of physical parameters varied is 7 and we consider the normalization to be an additional free parameter. Hence the number of degrees of freedom (d.o.f) is N d = 33 − 7 − 1 = 25 in this model, since the fitting is done to a total of 33 data points. All the parameter values for the best-fit case of the single-population model are listed in Table 1.
We see that the best-fit 1 H fraction turns out to be zero, and a non-zero 56 Fe component is unavoidable in this case. Indeed from the best-fit spectrum, shown in the upper left panel of Fig. 1, the contribution from Z = 1 component above 5 · 10 18 eV is infinitesimal. Since the heavier nuclei must come from nearby sources, for them to survive at the highest energies, the maximum rigidity, in this case, suggests that the cutoff in the spectrum originates from maximum acceleration energy at the sources. The fit, however, corresponds to a negative injection spectral index, which is difficult to explain by either the existing particle acceleration models or by sufficient hardening due to photohadronic interactions in the environment surrounding the source. The slope of the simulated X max plot (cf. Fig. 1), in comparison to data, suggests that the addition of a light element above 10 19 eV can improve the fit. Motivated by these aforementioned characteristics of the combined fit, it is impulsive to add the contribution from another source population and check the effects on the spectrum and composition.

Two-population model
We consider a discrete extragalactic source population injecting 1 H following the spectrum of Eq. (1). We refer to this as the source-class I (abbv. Cls-I). This pure-proton component has a distinct rigidity cutoff R cut,1 , and injection spectral index α 1 2, such that the spectrum extends up to the highest-energy bin of the observed UHECR spectrum. The normalization A 1 = A p is fixed by the condition the observed spectrum and E h is the mean energy of the highest-energy bin. f H is an additional parameter that takes care of the proton fraction in the highest-energy bin of the UHECR spectrum.
Another population (source-class II, abbv. Cls-II) injects light-to-heavy nuclei, viz., 4 He, 14 N, 28 Si, and 56 Fe, as we have already seen that for a mixed composition at injection, the contribution of 1 H abundance tends to be zero above the ankle energy. Cls-II also follows the spectrum in Eq. (1) with rigidity cutoff R cut,2 and injection spectral index α 2 , and the abundance fraction at injection given by K i ( i K i = 100%). The normalization A 2 in this case is a free parameter which is adjusted to fit the spectrum and composition. As in the single-population model, here too, we set the maximum redshift of the sources to z max = 1. Athough the anisotropy of UHECR arrival directions suggest that the observed spectrum depends on the position distribution of their sources, a definitive source evolution model is difficult to find. The rigidity cutoff and the injection spectral index will vary widely with the variation of evolution function and its exponent. We first consider that both the source populations are devoid of redshift evolution, i.e., m = 0 in the (1 + z) m type of source evolution models. Afterwards, in the next subsection, we present the m = 0 cases for one-and two-population models.
The cumulative contribution of Cls-I and Cls-II is used to fit the UHECR spectrum and composition for fixed values of f H . We vary f H from 1.0 to 20.0%, at intervals of 0.5% between 1.0 and 2.5% and at intervals of 2.5% between 2.5 and 20.0% to save computation time. α 1 is varied through the values 2.2, 2.4, and 2.6, inspired by previous analyses with light elements fitting the UHECR spectrum [70,71]. We vary log 10 (R cut,1 /V) between the interval [19.5, 20.2] at grid spacings of 0.1, and log 10 (R cut,2 /V) between [18.22, 18.36] 1 , f H }, we find the best-fit values of log 10 (R cut,1 /V), log 10 (R cut,2 /V), α 2 , and composition K i at injection of Cls-II; that minimizes the χ 2 tot of the combined fit. Due to increased number of parameters, we set the precision of composition K i to 0.25%. These parameter sets are listed in Table 2. For α 1 = 2.2 and 2.4, the χ 2 tot value monotonically increases with f H beyond the best-fit value, while for α 1 = 2.6, an alternating behaviour is obtained. The best-fits are found at f H = 1.5%, 2.5%, and 2.0%, respectively for α 1 = 2.2, 2.4, and 2.6. For all the cases, a significant improvement in the combined fit is evident compared to the one-population model. It is worth pointing out that the minimum of χ 2 comp and χ 2 spec do not occur simultaneously and the variation in the best-fit value of log 10 (R cut,2 /V) is insignificant. In the top right and bottom panels of Fig. 1, we show the best-fit It is instructive to compare the all-flavor neutrino fluxes resulting from the two-population model with the current 90% C.L. differential flux upper limits imposed by 9-years of IceCube data [60]. The hard spectral index and lower maximum rigidity in case of one-population model leads to a neutrino spectrum much lower than the current and upcoming neutrino detector sensitivities. This is shown in the top left panel of Fig. 2 along with the current sensitivity by PAO [72,73] and that predicted for 3-years of observation by GRAND [74,75] and POEMMA [76,77]. We also present the allowed range of neutrino flux from Cls-I and Cls-II in the two-population model for f H = 1.0-20.0%. The cosmogenic neutrino flux from Cls-I is within the reach of the proposed GRAND sensitivity. The all-flavor integral limit for GRAND implies an expected detection of ∼ 100 neutrino events within 3-years of observation for a flux of ∼ 10 −8 GeV cm −2 s −1 sr −1 . This implies that with a further increase in exposure time, GRAND should be able to constrain our two-population model parameters if f H 10%.
As we find the best-fit H fraction is zero in Table 1, K H is a redundant parameter in this case. Scanning the parameter space excluding the latter will result in the same values of the remaining 6 parameters and thus, the resulting model coincides with that of Cls-II in Table 2. Thus, for a Δχ 2 calculation between the one-population and twopopulation model, we consider the number of parameters in the former to be 6 and not 7. The difference in the number of parameters varied between one-population and twopopulation model is one, i.e., R cut,1 . A smooth transition from the two-population model to one-population model can be done by setting R cut,1 = 0. This necessarily implies that f H = 0 and there remains no α 1 . Based on the values obtained from, we estimate the maximum allowed proton fraction at 3.5σ confidence level (C.L.) in the highest-energy bin. For α 1 = 2.2 this corresponds to ≈ 12.5%, α 1 = 2.4 corresponds to ≈ 15.0%, and for α 1 = 2.6 it turns out to be ≈ 17.5%. However the maximum | Δχ 2 |, which also indicates the most significant improvement in contrast to one-population model, is found for α 1 = 2.2, as shown in Fig. 3. The 2.6σ and 3.5σ C.L. are also indicated.

Redshift evolution of sources
In the preceding study with flat redshift evolution of the two populations of extragalactic sources, we see that the contribution of 1 H from the light-to-heavy nuclei injecting sources, to the combined fit of energy spectrum and mass composition beyond the ankle, is infinitesimal. Whereas, the pure proton spectrum from Cls-I maintains a steady contribution up to the GZK energies superposed on the Peters cycle pattern [78], resulting from Cls-II. We carry out a systematic analysis over plausible strengths of redshift evolution of the source classes. We assume the source distribution evolves with red- First, we find out the best-fit value of m in the one-population model assuming a mixed composition at injection comprising of 4 He, 14 N, 28 Si, and 56 Fe. The combined fit improves with comparison to the flat evolution case, but not substantially. The resulting spectrum and composition fit are shown in the top left panel of Fig. 4. The composition fit is found to be more significant than the flat evolution case. In this case too, we see the contribution from protons, resulting in the photodisintegration of heavier elements, is sub-dominant at E 10 18.7 eV. The redshift evolution index m is varied in the range −6 m +6 at intervals of 1.0 and the corresponding best-fit values of cutoff rigidity, injection spectral index, and the composition are calculated. The best-fit parameters and the fit statistics are indicated in Table 3. The number of d.o.f is 25. The minimum χ 2 is obtained for the fit corresponding to m = +2. This indicates a wide range of candidate classes, eg., low-luminosity GRBs (LL GRBs) where the UHECR nuclear survival is possible inside the source/jet [20]. However, their redshift evolution is not well known but expected to follow that of long GRBs, given by ψ(z) ∝ (1 + z) 2.1 for 0 < z < 3 [79].
For the two-population case, we need to take into account two values of the redshift evolution index m 1 and m 2 , respectively for Cls-I and Cls-II. In case of high-luminosity γray sources, the dynamical timescales are larger than nuclear interaction timescales, inside the acceleration region. The relativistic jet provides suitable environment for heavier nuclei to dissociate via interactions with ambient matter and radiation. Hence, they are ideal candidate for 1 H injection. We identify our Cls-I with AGNs, injecting predominantly protons. The redshift evolution of AGNs follow the function ψ(z) ∝ (1 + z) 3.4 for z < 1.2 and X-ray luminosity in the range L X ∼ 10 43 − 10 44 erg/s [80]. An even higher luminosity might be required to accelerate UHECR protons up to 10 20 eV [81]. In principle, one can consider even higher luminosity AGNs, but the number density decreases sharply with luminosity. The redshift evolution of medium-high luminosity AGNs (L X ∼ 10 44 − 10 45 ) is given by ψ(z) ∝ (1 + z) 5.0 for z < 1.7. Radio-loud quasers with bolometric γ -ray luminosity 10 47 erg/s and a number density of 10 −5 −10 −4 Mpc −3 can meet the energy requirements for UHECR acceleration [82]. Hence, we vary m 1 through 3, 4, and 5. While for m 2 , we consider a wide range of values spanning from positive to negative, viz., +2, 0, −3, and −6, to find the best-fit region. Once again, we fix the injection spectral index of protoninjecting sources to α 1 = 2.2, 2.4, and 2.6, which is now a more physically motivated choice for AGNs. As before, we vary the proton fraction ( f H ) at the highest energy bin from 1.0 to 2.0% at intervals of 0.5%, and from 2.5 to 10.0% at intervals of 2.5%. We find the best-fit value of R cut,1 , R cut,2 , α 2 , and the fractional abundance of elements at injection (K i ) for Cls-II. For this case, we vary log 10 (R cut,2 /V) with a precision of 0.1.
We represent the goodness-of-fit for the best-fit case corresponding to each set of {α 1 , f H , m 1 , m 2 } values in Fig. 5, distinctly for the combined fit, spectrum fit, and the composition fit from top to bottom, respectively. The combined fit improves as we go to more and more negative values of m 2 and lower values of α 1 . We see the best-fit occurs for m 1 = +3, m 2 = −6, α 1 = 2.2 and f H = 1.5%. The details for this set are given in Table 4. The best-fit spectrum and composition are displayed on the top right panel of Fig. 4. However, it is interesting to note from Fig. 5 that the best-fit composition and best-fit spectrum cases are not coincident. The best-fit composition (χ 2 comp = 10.29) is obtained for m 1 = +5, m 2 = 0, α 1 = 2.2 and f H = 1.5%, whereas the best-fit spectrum (χ 2 spec ) occurs at m 1 = +5, m 2 = −6, α 1 = 2.6 and f H = 2.5%. We also calculate the neutrino fluxes originitaing from the best-fit one-population and twopopulation models, after considering redshift evolution of the source classes. In the bottom left panel of Fig. 4, the neutrino flux increases in the case of one-population model, owing to the positive redshift evolution (m = 2). While, in case of two-population model, the flux from heavy nuclei injecting sources is greatly reduced (m 2 = −6), and the cumulative neutrino flux distribution is dominated by that from protons (m 1 = 3). The shaded region indicates the flux range enclosed by f H = 1.0−20.0%, α 1 = 2.2, and is within the flux upper limit imposed by 9-yr of IceCube data. Even a small fraction of proton can yield a neutrino flux which is within the reach of 3-yr extrapolated sensitivity of the proposed GRAND detector.
The preference over large negative values of m 2 can be attributed to specific source classes, such as tidal disruption events (TDEs) [21,83]. The event rate of TDEs depend on the number density of SMBH as a function of redshift. The bestfit empirical model indicates a negative redshift evolution   [84]. TDEs forming relativistic jets can be the powerhouse of UHECR acceleration, but their event rate severely constrains the UHECR flux [85], thus requiring a mixed or heavy composition at injection. Metal-rich composition consisting of a significant Si and Fe fraction is required to explain the spectrum with a population of TDE [24]. In our case too, a high fraction of Fe is required to explain the spectrum at the highest energies, as indicated in Table 4. However, the survival of UHECR nuclei depends on the specific outflow model. For luminous jetted TDEs like Swift J1644+57 [86,87], which reaches a bolometric luminosity L bol 10 48 erg/s in the high state, UHECR acceleration becomes difficult via internal shock model, but is allowed for TDEs with lower luminosities. Forward/reverse shock models were also found in accordance with heavy nuclei injection [88].
The best-fit obtained in two-population model with nontrivial redshift evolution (m 1 , m 2 =0) is better than the flat evolution case and also compared to the one-population Table 4 Best-fits to UHECR spectrum and composition for two-population model (m 1 = +3, m 2 = −6) model with redshift evolution. A smooth transition can be made from two-population model to one-population model by setting R cut,1 = 0 and m 1 = 0. So the difference in the number of parameters varied is two. Figure 6 shows the | Δχ 2 | values for two d.o.f between the one-population and two-population cases, as a function of proton fraction at the highest-energy bin. As before, we constrain the maximum allowed proton fraction at 3.5σ confidence level, which turns out to be between 7.5 and 15% for different values of the proton injection spectral index considered. We also calculate the correlation between the fit parameters for the specific case of m 1 = 3, m 2 = −6, α 1 = 2.2, and f H = 1.5%, i.e., for the best-fit case SRE-1 listed in Table 4. We vary the cutoff rigidity log 10 Figure 7 shows the 1σ , 2σ , and 3σ C.L. contours for 25 d.o.f. The variation of cutoff rigidity is not shown, since they were found to be insensitive to variation of other parameters. We see a hard injection spectral index α 2 ≈ 1.6 is preferred, which conforms with the recent predictions by the PAO [43], and analysis done by other works [31,32]. The composition is also in accordance with those predicted from latest measurements, implying a progressively heavier composition at higher energies. The 1σ region in composition space corresponds to a high value of Fe fraction, which is indeed needed for TDEs to have significant contribution in the UHECR spectrum. Another candidate class which can represent our Cls-II is the lowluminosity (L γ < 10 44 erg/s) and high synchrotron peaked BL Lacertae objects. They possess a negative redshift evolution and are predicted to be more numerous than their highluminosity counterpart. However, a direct detection of these low-luminosity objects are difficult, and the current 4LAC catalog consists ∼ 20 such sources.

Discussions
The composition fit corresponding to the one-population model, especially the departure of simulated X max and σ (X max ) values from the data, leaves a substantial window for improvement. The addition of a light nuclei component up to the highest observed energies shall alleviate the mismatch. We exploit this possibility in our work by adding a distinct source population injecting 1 H that extends up to the highest observed energies. Earlier works have considered a pure-protonic component with an assumed steep injection spectral index [54] to explain the region of the spectrum below the ankle. We do not fit the UHECR spectrum below the ankle, and the proton spectrum considered in our work contributes directly to the improvement in composition fit at the highest energies. A relatively hard spectrum (α 1 = −1) in addition to a Milky Way-like nuclear composition is considered in Ref. [55], extending up to the highest energies. We do not assume any fixed abundance fraction for the lightto-heavy nuclei injecting sources and calculate the best-fit values within the resolution adopted.
Reference [57] have proposed an interaction-model independent method to probe the allowed proton fraction for E p 30 EeV, constrained by the cosmogenic neutrino flux upper limits at 1 EeV. Thus, they do not take the composition of primary cosmic rays into account, inferred from air shower data. They have considered a generalized redshift evolution function of the proton injecting sources, parametrized by the evolution index m. In our work, we fit the composition data X max , σ (X max ), and the energy spectrum simultaneously to infer the proton fraction in a two-population scenario.
Here, we find a significant improvement in the combined fit to spectrum and composition data, when adding an extragalactic source population emitting UHECRs as protons. For our choice of steep proton injection indices (α 1 ), the goodness-of-fit is found to be comparable to each other. We also consider the injection index (α 2 ), maximum rigidity (R cut,2 ), and composition fractions (K i ) of the second population injecting light-to-heavy nuclei to be variables and find the corresponding best-fit values. The corresponding improvement in the combined fit is found to be 3σ in some cases.  We have also surveyed our results for a wide range of source redshift evolution. Such an analysis is already done earlier for a single source population for a mixed composition of injected elements [31,32,89]. In our analysis, we find that, although a positive evolution index is preferred in onepopulation model, the best-fit value changes sign on going to two-population model. However, with increasing values of z max , the variation of m can significantly affect the neutrino spectrum. We have kept the contributing sources within z 1 in view of the fact that particles originating at higher redshifts will contribute below the ankle, which we do not fit here. Thus within the minimal requirements of this model, our neutrino spectrum can be considered as a conservative lower bound in the two-population scenario.
The resultant neutrino spectrum in two-population model at E 0.1 EeV is dominated by that from pure-protons. Even a small fraction of protons at the highest energy is capable of producing a significant flux of neutrinos. This is expected because of the maximum energy considered for proton-injecting sources. Even for low f H , the values of E max are very close to GZK cutoff energy, where the resonant photopion production occurs, leading to pion-decay neutrinos. The double-humped feature of the neutrino spectrum is a signature of interactions on the CMB and EBL by cosmic rays of different energies. The higher energy peak produced from protons possesses the highest flux, and the detection of these neutrinos at ∼ 3 · 10 18 eV will be a robust test of the presence of a light component at the highest energies, thus also constraining the proton fraction. For E < 0.1 EeV, the neutrinos from Cls-II becomes important with peaks at ∼ 1 PeV and ∼ 40 PeV. Hence, the cumulative neutrino spectrum (Cls-I + Cls-II) exhibits three bumps for α 1 = 2.2 (see Fig. 2). But gradually with increasing values of α 1 , the lower energy peak of Cls-I becomes significant, diminishing the "three-peak" feature until neutrinos from protons dominate down to ∼ 1 PeV for α 1 = 2.6 Fig. 7 Correlation between fit parameters for the best-fit case corresponding to m 1 = 3, m 2 = −6, α 1 = 2.2, and f H = 1.5%, i.e., for the best-fit case SRE-1 listed in Table 4. It can be seen that a high fraction of Fe is required at injection along with a hard injection spectral index. The diagonal plots represent the posterior probability distribution and red dots in others indicate the central values. The 1σ , 2σ , and 3σ standard deviations are shown by dark to light-colored shading We present the upper limit on the maximum allowed proton fraction in two-population model at ≈ 1.4×10 20 eV. This is based on the improvement in the combined fit compared to the one-population model, up to 3.5σ statistical significance. For a higher C.L., the proton fraction is even lower at the highest-energy bin. However, a non-zero proton fraction is inevitable. It is studied earlier that the flux of secondary photons increases with an increasing value of α 1 [33]. If a single population injecting protons is used to fit the UHECR spectrum, the resulting cosmogenic photon spectrum saturates the diffuse gamma-ray background at ∼ 1 TeV for α 1 = 2.6, m = 0 [71]. In our two-population model, the proton fraction at the highest energies is much lower than the total observed flux. This ensures the resulting photon spectrum from Cls-I is well within the upper bound imposed by Fermi-LAT [59]. For Cls-II injecting heavier nuclei, the main energy loss process is photodisintegration, contributing only weakly to the cosmogenic photon flux. Hence the two-population model, which we invoke in our study, is in accordance with the current multimessenger data.
The choice of the hadronic interaction model for our analysis is based on the interpretation of air shower data by the PAO [58,66]. It is found that qgsjet-II.04 is unsuitable compared to the other two models and leads to inconsistent interpretation of observed data [49]. Also, for our choice of photodisintegration cross-section, i.e. talys 1.8, the hadronic model sybill2.3c yields superior fits [32]. In general, the sybill2.3c model allows for the addition of a higher fraction of heavy nuclei, compared to others, at the highest energies. Indeed in Table 2, it is seen that the lowest-χ 2 cases correspond to high K Fe , which increases monotonically with α 1 . The requirement of Fe abundance in one-population model is much lower than in the case of two-population model. For the latter, the cutoff in the cosmic ray spectrum cannot be solely explained by the maximum acceleration energy of iron nuclei at the sources, but also, must be attributed to photopion production of UHECR protons on the CMB to some extent.
In going from one-population to the two-population model, the injection spectral index of the population injecting heavier elements changes sign from negative to positive, making it easier to accept in the context of various astrophysical source classes. Young neutron stars, eg., can accelerate UHECR nuclei with a flat spectrum, α 2 ∼ 1 [90]. Particle acceleration in magnetic reconnection sites can also result in such hard spectral indices [see for eg., 91]. Luminous AGNs and/or GRBs are probably candidates for Cls-I, accelerating protons to ultrahigh energies [15]. The Cls-II injecting light-to-heavy nuclei suggests the sources to be compact objects or massive stars with prolonged evolution history, leading to rich, heavy nuclei abundance in them. In particular the high negative redshift evolution and substantial Fe fraction allows us to identify the Cls-II with TDEs. The problem in the case of a highly luminous object is, although heavier nuclei may be accelerated in the jet, they interact with ambient matter and radiation density in the environment near the sources [92]. To increase the survivability of UHECR nuclei, less luminous objects such as LL GRBs [93] are preferred.

Conclusions
Based on the spectrum and composition data measured by PAO, a combined fit analysis with a single-population of extragalactic sources suggest that the composition fit at the highest energy deserves improvement. The slope of the simulated X max curve implies that fitting the highest-energy data points with contribution from only 56 Fe will diminish the abundance of lighter components 28 Si, 14 N, and 4 He. This will in turn decrease the flux near the ankle region, thus resulting in a bad fit. Addition of another light component of extragalactic origin, preferably pure proton, extending up to the highest-energy bin can resolve this problem. From a critical point of view, this solution is not unique, but definitely a rectifying one. The combined fit improves significantly and we present the maximum allowed proton fraction at the highest-energy bin of spectrum data corresponding to > 3σ statistical significance. An additional population of extragalactic protons has also been suggested in Ref. [94], in the context of fitting the UHECR spectrum.
There are observational indications that different astrophysical source populations likely contribute to the UHECR data. A plausible hot spot around the nearby starburst galaxy M82 [95,96] in the TA data and an intermediate-scale anisotropy around the nearest radio galaxy Cen-A in the Auger data already suggest possibility of two types of source populations. The Auger data, however, do not show any small-scale anisotropy, suggesting that the majority of UHECR sources are distributed uniformly in the sky. The recent 3σ correlation of an observed high-energy muon neutrino event detected by IceCube with a flaring blazar TXS 0506+056 at a moderate redshift of 0.34 is consistent with this scenario [97,98].
The two generic source classes of UHECRs studied here by us is also representative of the scenario described above. High luminosity AGNs or GRBs could contribute a pure proton component that is significant at the highest UHECR energies. The resulting cosmogenic neutrino spectrum can be detected by future experiments with sufficient exposure and the proton fraction in the highest energy UHECR data can be tested.