figurec

IceCube Collaboration

In this Letter, we report a search for EHE$\nu$ using $12.6\text{\,}\mathrm{years}$ of data from the IceCube Neutrino Observatory. The data were taken between June 2010 and June 2023, corresponding to 4605 days of livetime. This is $50$% more exposure than that of the previous IceCube search [23]. Additionally, the event selection has been reoptimized, improving the effective area by $\sim$15$\%$ near $1\text{\,}\mathrm{EeV}$. The null observation of cosmogenic neutrinos places significant constraints on the cosmological evolution of UHECR sources and, moreover, the composition of UHECRs. In this work, we investigate the specific hypothesis of a proton-only composition of UHECRs with the GZK cutoff generating the observed suppression of UHECRs at $\mathrm{EeV}$ energies. The method was proposed in [9], and in a similar fashion applied to the neutrino measurement by the Pierre Auger Collaboration [24]. We find, at $90$% CL, that the observed fraction of UHECRs that are protons at Earth above $\simeq$30\text{\,}\mathrm{EeV}$$ cannot exceed $70$% if the source evolution is comparable to the star formation rate (SFR)—a general tracer of matter density in the Universe. These constraints are complementary to, and agree with, direct air-shower measurements [25, 26] by being insensitive to uncertainties associated with hadronic interaction models.

Data Sample—IceCube [27] is a neutrino detector at the South Pole. It consists of $5160$ digital optical modules (DOMs), distributed on $86$ strings, instrumenting a cubic kilometer of ice at depths between $1450$ and $2450\text{\,}\mathrm{m}$. Each DOM hosts a 10-inch photomultiplier tube [28] and readout electronics [29]. Charged particles produced in neutrino interactions give rise to Cherenkov light when propagating through the ice. Those Cherenkov photons are detected by DOMs and converted into photoelectrons ($\mathrm{PE}$). On top of the IceCube strings, a surface array called IceTop [30] measures cosmic-ray air showers. EHE$\nu$ events in IceCube are observed as either tracks—light depositions along the trajectory of a long-range $\mu$/$\tau$ produced in $\nu_{\mu}$/$\nu_{\tau}$ charged-current interactions—or cascades—approximately spherical light depositions arising from all-flavor neutral-current interactions and charged-current interactions of $\nu_{e}$.

This search aims to select cosmogenic neutrinos while rejecting backgrounds. Because $\gtrsim 100$ $\mathrm{TeV}$ neutrinos are absorbed by Earth [31], EHE$\nu$s are expected to be downgoing or horizontal at IceCube. The dominant background is downgoing atmospheric muon bundles produced in cosmic-ray air showers. This flux is modeled using corsika [32], with sibyll2.3c [33] as the hadronic interaction model and the cosmic-ray flux prediction from [34]. Further backgrounds arise from atmospheric and astrophysical neutrinos. Atmospheric neutrinos are produced by meson decays during cosmic-ray air showers. Their flux is divided into a conventional component [35] originating from pion and kaon decays, and a yet-unobserved prompt component [36] produced by heavier, short-lived mesons. All neutrinos—atmospheric, astrophysical, and cosmogenic—are simulated using the juliet code [37].

The $3\text{\,}\mathrm{kHz}$ event trigger rate in IceCube is dominated by atmospheric muons, while the cosmogenic neutrino flux is already constrained to $\ll 1$ event per year [23]. The signal-to-noise ratio is improved by employing an event selection based on quality cuts of high-energy events in combination with an IceTop veto. Full details of the event selection are presented in End Matter. In brief, cosmogenic signal events have extremely high energies and, therefore, produce large amounts of charge ($\mathrm{PE}$). As the atmospheric muon background is exclusively downgoing, we reject most backgrounds with a zenith-angle-dependent charge threshold. To do this, the event direction is reconstructed with a maximum-likelihood reconstruction using an infinite-length track hypothesis [38].

Energy loss profiles for single muons show large stochastic variations as the muons propagate. In contrast, in high-multiplicity muon bundles, these single-muon fluctuations partly average out. To leverage this, the energy loss profile of events is reconstructed along the track direction, and stricter charge requirements are imposed upon less stochastic events. Using stochasticity information improved the effective area by $15$% between $100\text{\,}\mathrm{PeV}$ and $1\text{\,}\mathrm{EeV}$ relative to the previous search. Lastly, IceTop is used to further reduce the background from atmospheric events, as described in [39]. The sample is divided into subsamples of tracks and cascades based on the reconstructed particle velocity, and their deposited energy and arrival direction are reconstructed with likelihood-based methods [40, 41] (cf. End Matter).

After the event selection, the expected atmospheric background is $$0.40$\pm 0.03$ events, and up to $\sim 5$ cosmogenic neutrinos are expected for the most optimistic model [10], consisting of $73$% tracks and $27$% cascades. The flux beyond $\mathrm{PeV}$ energies is not well constrained, and the expectation strongly depends on the assumed model. The astrophysical expectation ranges from $\sim 9$ events for an unbroken power law with a hard spectral index ($\gamma=2.37$) [42], down to $\sim 0.5$ events assuming a power law with a cutoff ($\gamma=2.39$, $E_{\mathrm{cutoff}}=$1.4\text{\,}\mathrm{PeV}$$ ¹¹1These best-fit parameters follow from the analysis scheme in [44], and they will be published in future work.). At the highest energies, $E_{\nu}>$100\text{\,}\mathrm{PeV}$$, this expectation is reduced to $0.9$ events and $3\text{\times}{10}^{-30}$ events, respectively.

For both astrophysical and cosmogenic neutrinos, a flavor ratio of $\nu_{\mathrm{e}}:\nu_{\mu}:\nu_{\tau}=1:1:1$ at Earth is assumed [45], as well as equal fluxes of neutrinos and antineutrinos.

Three events with $\mathrm{PeV}$ energies survive the event selection: a through-going track [46], an uncontained cascade [47], and a starting track [48].

Analysis and Results—To infer physics parameters, the data are fit using a binned Poisson likelihood in the space of reconstructed direction and energy, following the method in [23]. Being in the regime of small statistics, all hypothesis tests and limits are based on ensembles of pseudoexperiments. Confidence intervals are determined using the likelihood ratio test statistic [49].

Systematic uncertainties are treated similarly to [45]: They are varied in pseudoexperiments based on estimated priors. The effect of incorporating systematics is modest, with the differential limit weakened by $4$% at $100\text{\,}\mathrm{PeV}$, reducing toward $1$% at the highest energies. Details of the likelihood and systematics treatment are available in Supplemental Material [50].

Differential limit and model tests—The differential upper limit on the neutrino flux above $5\text{\times}{10}^{6}\text{\,}\mathrm{GeV}$ is depicted in Fig. 1 as a red line. The sensitivity, i.e., the limit in case of a null observation, and previous limits are also shown. The limit is weakened with respect to the sensitivity below $100\text{\,}\mathrm{PeV}$ due to the observed events. At the highest energies the improvement to the previous limit (IceCube 9yr [23]) is comparable to the increase in detector livetime as expected for a largely background-free analysis. At around 1 EeV, additionally, event selection enhancements make a sizable contribution. Notably, by a few hundred $\mathrm{PeV}$, the previous limit deviates from the sensitivity due to observed $\mathrm{PeV}$ neutrinos. In this analysis, although a similar number of $\mathrm{PeV}$ neutrinos were observed, improved energy reconstructions mean they are incompatible with a flux centered at $100\text{\,}\mathrm{PeV}$ and above.

The differential limit is compared to a representative variety of cosmogenic neutrino models as gray lines. Qualitatively, larger normalizations, higher maximum accelerating energies, and stronger source evolutions generate larger cosmogenic neutrino fluxes [12, 11, 59]. In contrast, when the injected cosmic-ray primaries are heavy nuclei, photodisintegration becomes the dominant process over photopion production and the neutrino flux is suppressed [6, 60, 13]. All model predictions shown in Fig. 1 (except Van Vliet et al. 2019 [20], abbreviated “vV2019” hereafter) assume a pure-proton composition with moderate source redshift evolutions comparable to the SFR. The maximum acceleration energy varies between ${10}^{11}\text{\,}\mathrm{GeV}$ and ${10}^{12}\text{\,}\mathrm{GeV}$.

For each aforementioned model, we performed a likelihood ratio test as described in Supplemental Material [50]; the results are in Table 1. Although three events were observed, the best-fit normalization for a cosmogenic flux component is zero for all tested models. This indicates the data can be sufficiently explained by astrophysical neutrinos. All tested cosmogenic models assuming a pure proton composition of UHECRs are rejected at $95$% CL. This indicates that regardless of the differences between those models, if the SFR is driving the source evolution of UHECRs, a proton-only composition can be excluded.

Proton fraction constraints—Given the measured UHECR flux, the nonobservation of neutrinos imposes constraints on the sources. This approach is complementary to many existing models, which focus on accurately describing the cosmic-ray energy spectrum and composition and, thus, also obtain an estimation of the accompanying cosmogenic neutrino flux [15, 19, 65, 66].

The crpropa package [67] is used to model cosmogenic fluxes (following vV2019 [20]). In the simulation, protons and secondary neutrinos are propagated to Earth including energy losses from photopion production and pair production on the CMB and EBL [68], neutron decay and cosmological adiabatic losses. Identical sources are distributed homogeneously and isotropically with a power-law injection spectrum $\Phi(E)\propto E^{-\alpha}\exp(-E/E_{\mathrm{max}})$ with spectral index $\alpha\in[$1.0$,$3.0$]$ and exponential cutoff at $E_{\mathrm{max}}\in[$4\text{\times}{10}^{10}\text{\,}\mathrm{GeV}$,${10}^{14}\text{\,}\mathrm{GeV}$]$.

Two different models for cosmological source evolution are tested:

with $z^{\prime}=1.5$ up to $z_{\mathrm{max}}=4$ [20], and a more conservative model of $\mathrm{SE}_{2}(z)=(1+z)^{m}$ with $z_{\mathrm{max}}=2$, where $m$ denotes the source evolution parameter. The simulation is normalized to the all-particle cosmic-ray flux measured by Auger at $10^{10.55}\,$$\mathrm{GeV}$. We normalize to the highest-energy data point below the observed GZK suppression, such that the cosmic-ray flux at the suppression energy is saturated. This defines the flux corresponding to a proton fraction at Earth ($f_{p}$) of $100$% above energies of $\simeq$30\text{\,}\mathrm{EeV}$$. The reference energy impacts the resulting neutrino fluxes. A systematic shift to the reference energy on the order of the systematic energy scale of Auger of $\pm$14$\%$ [25] results in a $5$% shift of the overall neutrino flux.

Figure 2 shows the construction of the $f_{p}$ constraints. The light-colored histograms show the simulated proton flux saturating the Auger measurement and the secondary neutrino flux. The source parameters $\alpha$ and $E_{\mathrm{max}}$ are chosen to minimize the integral neutrino energy flux to obtain a conservative prediction for a given value of $m$. The relatively wide range for $\alpha$ is motivated by both experimental and theoretical work, e.g., the Auger Collaboration [69] showing that $\alpha$ between 1 and 2 are allowed. Allowing this wider range gives a slightly more conservative result than bounding $\alpha$ at 2. The range for $E_{\mathrm{max}}$ is bracketed by the “GZK-cutoff” energy at the low end, and by a value much higher than the observed cosmic rays on the other. In practice, when marginalizing, the minimum value for $E_{\mathrm{max}}$ is always chosen. The flux shown in the figure is in tension with IceCube data, and, thus, $f_{p}$ can be constrained based on the determined upper limit. As suggested in vV2019 [20], $f_{p}$ can be determined by comparing the predicted neutrino flux with the experimental limit at $1\text{\,}\mathrm{EeV}$. However, we instead perform a model test, which improves the sensitivity.

This procedure is repeated for different values of the source evolution parameter $m$, and the resulting constraints are shown in Fig. 3 for the source evolution models $\mathrm{SE}_{1}(z)$ and $\mathrm{SE}_{2}(z)$. The value of $m$ for UHECR sources is unknown, but here we focus on the range in which $f_{p}$ can be constrained by this analysis. For instance, given the source evolution is comparable to the SFR or stronger, $f_{p}$ is constrained to be below about $70$%. Alternatively, due to the degeneracy between $f_{p}$ and $m$, the results can be interpreted as an upper bound on the source evolution of $m\lesssim 3$ for proton-dominated UHECRs, strengthening the claim of the previous analysis [39].

The predicted neutrino fluxes are dominated by distant cosmic-ray sources, from which high-energy cosmic rays are not expected to survive. The model presented here assumes that cosmic-ray sources are distributed homogeneously within the Universe. This is true at large distances, but due to Earth’s position within the Local Supercluster and Local Sheet, the local density of sources is enhanced. This leads to a relative reduction of distant sources and, thus, of the expected neutrino flux. Including a model of the local overdensity of sources based on the star formation rate of local galaxies [25, 71] weakens the neutrino fluxes, and the corresponding proton fraction constraints, by about $3$% ($m=6.0$) to $4$% ($m=2.0$). Additionally, a recent study by Auger shows that up to $5$% of UHECRs above $40\text{\,}\mathrm{EeV}$ can be associated with Centaurus A as the dominant local source [69]; in this case, the $f_{p}$ constraint becomes weaker by the same fraction.

For context, constraints already exist from direct air-shower measurements. The Auger $X_{\mathrm{max}}$ data [25] and the TA event isotropy [26] both favor small proton fractions. However, the data allow [24, 72], and some analyses find [66, 65, 73], a proton fraction as high as $\sim$10$\%$. Air-shower-based composition measurements, while powerful, are dependent on hadronic interaction models and the associated uncertainties. As such, independent measurements that do not rely on air-shower observables directly are highly complementary and necessary. In particular, our result in this study does not rely on observables from air showers and is, therefore, insensitive to the uncertainties associated with hadronic interaction models. Although the estimation of the atmospheric muon or neutrino background has a dependence on hadronic interaction models, the influence on the derived constraints is negligible.

Auger has also used their neutrino data to constrain $f_{p}$ [24], similar to Fig. 3. Here, we substantially improve the constraints. In particular, a direct comparison with the $\mathrm{SE}_{2}$ model with $z_{\mathrm{max}}=2$ can be made, where the resulting exclusion contour—the orange shaded region in Fig. 3—is shifted to smaller values of $m$ relative to the Auger result (the black dotted line) by $\sim 1$. For example, at $m\sim 3.4$, which is comparable to the SFR, we find $f_{p}\lesssim 70\%$ at 90%CL, where the Auger result is fully compatible with unity. We note that the IceCube result achieves this improvement despite making very conservative modeling choices, e.g., marginalizing over $\alpha$ and $E_{\mathrm{max}}$.

AGN model constraints—In addition, we tested the active-galactic-nuclei (AGN) model from [63], instead of a cosmogenic flux model. The astrophysical neutrino flux described in this model cannot be explained by the observed sub-PeV astrophysical neutrinos (cf. Fig. 1). These UHE astrophysical neutrinos are indistinguishable from cosmogenic neutrinos event by event. The neutrino emission is based on observed photon fluxes, using phenomenological parameters like the cosmic-ray loading factor $\xi_{\mathrm{CR}}$. The modeled flux scales linearly with $\xi_{\mathrm{CR}}$, so the limit (cf. Table 1) can be interpreted as an upper limit of $\xi_{\mathrm{CR}}\leq 1.4$ and $\xi_{\mathrm{CR}}\leq 62$ for assumed CR spectral indexes of $\gamma=$2.0$$ and $2.3$, respectively. That the resulting MRFs are $<1$ indicates that inner jets of AGN are unlikely to be a dominant source for UHECRs in this model scenario.

KM3-230213A — Recently, KM3NeT published a $\sim$220\text{\,}\mathrm{PeV}$$ neutrino candidate [62]. The inferred diffuse flux, also shown in Fig. 1, assumes an $E^{-2}$ ranging from $72\text{\,}\mathrm{PeV}$ to $2.6\text{\,}\mathrm{EeV}$ and significantly exceeds the limits presented in this work. With the exposure of this analysis, this flux leads to an expectation of $\sim 70$ events, inconsistent with our nonobservation at $>10\sigma$; a transient source hypothesis could reduce this tension [74, 62, 75, 76]. Considering a joint fit between IceCube, Auger, and KM3NeT, the tension in the diffuse hypothesis is significantly reduced [77]. After repeating the joint fit with the IceCube exposure presented here, the probability of the joint fit resulting in one observed event in KM3NeT (with $\mu_{\mathrm{KM3}}=0.01$ expected events) and no events in both Auger ($\mu_{\mathrm{A}}=0.3$) and IceCube ($\mu_{\mathrm{IC}}=0.68$) is $\sim$$0.35$%. The corresponding goodness-of-fit $p$ value determined by the saturated Poisson likelihood test [78] is $0.4$% ($2.9$$\sigma$).

The impact on $f_{p}$ constraints depends on the neutrino’s origin. If produced in a neutrino source environment, the constraints would be unaffected. If cosmogenic [79], a combined analysis will weaken the inferred limits.

Summary—The nonobservation of neutrinos with energies well above $10\text{\,}\mathrm{PeV}$ in $12.6\text{\,}\mathrm{years}$ of IceCube data places the most stringent limit on cosmogenic neutrino fluxes to date, reaching a neutrino flux of $E^{2}\Phi_{\nu_{e}+\nu_{\mu}+\nu_{\tau}}\simeq${10}^{-8}\text{\,}\mathrm{GeV}\text{\,}{\mathrm{cm}}^{-2}\text{\,}{\mathrm{s}}^{-1}\text{\,}{\mathrm{sr}}^{-1}$$. Additionally, we provide the strongest constraints on the composition of UHECRs obtained by neutrino astronomy, disfavoring proton-only UHECRs if their sources are evolving with the SFR or stronger.