The Astrophysical Journal, 928:50 (14pp), 2022 March 20

© 2022. The Author(s). Published by the American Astronomical Society.

OPEN ACCESS

Received 2021 November 19; revised 2022 January 14; accepted 2022 January 18; published 2022 March 25

Original content from this work may be used under the terms of the Creative Commons Attribution 4.0 licence. Any further distribution of this work must maintain attribution to the author(s) and the title of the work, journal citation and DOI.

Introduction

The field of neutrino astronomy has gained momentum since the IceCube collaboration discovered a diffuse flux of astrophysical neutrinos in multiple detection channels [9, 3, 13]. Prominent examples of this continuous journey are the identification of a first joint source of high-energy gamma rays and neutrinos, TXS 0506+056 [6, 11], increasing hints for the emission of high-energy neutrinos from the radio galaxy NGC1068 [7], and detection of a particle shower at the Glashow resonance energy [1].

Measuring the total observed flux strength and energy spectrum of high-energy astrophysical neutrinos complements direct searches for neutrino sources and is important to understand the processes behind the acceleration and propagation of high-energy cosmic rays.

The origin of high-energy neutrinos has been predicted from nonthermal Galactic and extragalactic sources, as reviewed in Learned & Mannheim [48], Becker [23], Halzen & Klein [42], and Ahlers & Halzen [20]. As the diffuse flux detected by IceCube is close to isotropically distributed and does not follow the Galactic plane, Galactic sources have been disfavored, while still being discussed as a subdominant contribution; see Becker Tjus & Merten [24] for a summary. The prompt phase of gamma-ray bursts (GRBs) has also been excluded as a dominant neutrino source by the dedicated IceCube analyses [15]; [10].

However, prominent possible sources of neutrinos exist, for example, choked GRBs in dense environments [58]; [26]. A second promising source of high-energetic extragalactic sources is active galactic nuclei (AGNs; Murase et al. 2014; Kimura et al. 2015; Liu et al. 2018), including BL Lac objects (Tavecchio & Ghisellini 2015; Padovani et al. 2015). Tidal disruption events (TDEs) are a third promising source class of high-energy neutrinos and ultra-high-energy cosmic rays [34]; [31]; [40]; [51]; [59]; [27]; [60]. Here, particle acceleration could be driven by either a hidden jet, a hidden subrelativistic wind, emission from a hot corona above the accretion disk, or radiatively inefficient accretion flows (see [43] for a review). Starburst galaxies have also been considered a possibly contributing source class [50]; [63]. These models have in common that they expect neutrinos to be produced when a power-law-distributed population of cosmic rays interacts with gas or photon fields in the source or its vicinity in order to produce pions and kaons, which in turn, decay into neutrinos. These neutrinos would follow the same power law as the initial cosmic-ray population, however, effects like an energy-dependent cross section, breaks, or spectral features in the accelerated cosmic-ray spectrum can lead to deviations from a pure power law. In this paper, we present an improved measurement of the energy spectrum of astrophysical muon neutrinos, including models beyond the single power law and tests of a selection of the abovementioned models. Also, we update the sample of up-going muon tracks (θzenith>85∘\theta_{\mathrm{zenith}}>85^\circ) originating from the northern celestial hemisphere [3].

Data Sample

The IceCube Neutrino Observatory is a gigaton-scale Cherenkov detector embedded in the Antarctic ice at the geographical South Pole [10]. Its fundamental building blocks are 5160 digital optical modules (DOMs). These spherical detection units each include a 10″ (25.4 cm) photomultiplier tube suited to detect weak light signals [14], read-out, and digitization electronics [18], and are positioned on 86 vertical cable strings that are arranged in a hexagonal grid. In this analysis, the data from eight strings forming the denser instrumented region of the DeepCore array [16] are excluded in favor of a more homogeneous detector geometry. Neutrinos are detected indirectly via the Cherenkov light emitted by charged relativistic secondary particles. These emerge from deep inelastic neutrino–nucleon interactions which occur in the instrument and the surrounding ice or the bedrock below IceCube. The final detector configuration (IC86) was completed in 2010 December. Recording of data already occurred prior to finalization with partial configurations. Data from two such partial configurations, IC59 and IC79, are included in this analysis, with the number indicating the number of active strings.

The presented study is based on a sample of through-going and starting muon tracks detected between 2009 May 20 and 2018 December 31. Muons can travel through the ice for hundreds of thousands of meters, producing track-like signatures in the detector. The event selection is identical to the one presented in [3]. It focuses on events with extended track-like signatures and filters out cascade-like events, which are events where the light emission region is small compared to the detector string spacing of 125 m. Such events can, for example, arise from neutral current interactions or electron–neutrino-induced charged current interactions. Neutrino-induced muon tracks are dominated by charged current interactions of muon neutrinos. A small additional energy-dependent contribution arises from muonic tau decays in tau neutrino interactions which contributes roughly 5% [62]. This is included in the signal expectation of the analysis. The contribution of tracks from Glashow resonance interactions [1] is negligible in this analysis. By restricting the field of view of the detector to zenith angles Notes. See text for a description of the Pass-2 re-calibration campaign that was conducted for all data after 2010 May.

The livetime for the season IC79 has been corrected by a factor 0.94 to account for the lower expected trigger rate of the partial detector configuration.

θzenith>85∘\theta_{\mathrm{zenith}} > 85^\circ, the overwhelming background of atmospheric muons from cosmic-ray-induced air showers is successfully suppressed (purity >99.8%> 99.8\%) as these muons are stopped in the Earth or the ice overburden before reaching IceCube. The direction of the selected muon-track events is obtained using the MPE algorithm with tabulated ice properties [21]; [56] and the muon energy is reconstructed using the truncated energy (DOM’s method) algorithm described in [17] which takes the stochastic energy losses of high-energy muons into account, resulting in an energy resolution of 0.22 in log⁡10Eμ\log_{10} E_\mu [17]. The resulting proxy energy Eμ,proxyE_{\mu,\mathrm{proxy}} is inherently only a lower limit on the true muon and primary neutrino energies, as the neutrino–nucleon interactions can occur outside the detection volume.

Compared to the previous iteration of this analysis based on six years of data, more than 300,000 new events were detected and analyzed; see Table 1. In addition, the full sample of events was reprocessed within the scope of IceCube’s Pass-2 campaign [8] to apply the latest detector calibrations to archival data as well. In this campaign, the same event filtering, selections, and reconstructions are consistently applied to all data after 2010 May (IC79−IC86-2018), while IC59 is still treated separately in the analysis. More details can be found in [62].

Data-taking SeasonZenith Range (deg)Effective Livetimeᵃ (years)Number of EventsPass-2
IC5990 − 1800.953821411No
IC79−IC86-201885 − 1808.181651377Yes
Total9.135672788

Table 1. Summary of the Event Sample

Five new events with reconstructed energies Eμ,proxy≳200E_{\mu,\mathrm{proxy}} \gtrsim200 TeV are observed in the additional data-taking period. This energy approximately corresponds to a zenith-averaged signalness S≳0.5S \gtrsim0.5:

S(Eμ)=Φsignal(Eμ)Φsignal(Eμ)+Φbackground(Eμ)≳0.5.(1)S(E_\mu) = \frac{\Phi_{\mathrm{signal}}(E_\mu)}{\Phi_{\mathrm{signal}}(E_\mu) + \Phi_{\mathrm{background}}(E_\mu)} \gtrsim0.5. \tag*{(1)}

So, for these events, the probability of the particle belonging to the differential astrophysical signal flux Φsignal\Phi_{\mathrm{signal}} exceeds 0.5.

The event with the highest reconstructed energy from this additional period is the horizontal track event with ID II.32 (Table 7) which has an energy of more than Eμ,proxy≳1.1E_{\mu,\mathrm{proxy}} \gtrsim1.1 PeV, and has also been reported as a real-time alert [10]. The event with the overall highest reconstructed energy has the table ID II.22, and is described in detail in [3]. The threshold Eμ,proxy≳200E_{\mu,\mathrm{proxy}} \gtrsim200 TeV to report individual events does not change significantly with the updated best-fit parameterization for the energy spectrum discussed in Section 4.4 The modified calibrations, however, which include a −4%-4\% shift of the charge corresponding to a single photoelectron, do lead to changes of the reconstructed energies and directions of some of the events that have been reported in [3], shown in the rightmost three columns of Table 7. The impact on the reported quantities is small, especially for the directional reconstructions (most directions change by less than 0.1∘0.1^\circ and stay well within their uncertainty ranges), but for single events with special topology, deviations in reconstructed energy can occur. Notably, such deviations arise for data taken in 2010 (IC79), where the information from the DeepCore strings is now excluded. As a consequence, six of the previously reported events fall below the threshold of Eμ,proxy≳200E_{\mu,\mathrm{proxy}} \gtrsim200 TeV. Additionally, one event (ID6 in [3], IC79-season) does not pass the event selection anymore, because the IC79 season data as used in [3] were based on an individual dedicated event selection, which has now been unified with the treatment of later seasons where the event selection includes additional measures to reject cascade-like events more efficiently. These had not been applied to the IC79 data previously. For a detailed list of all events and the respective changes, see Table 7 in the supplementary material.

Table of Observed Events with Eμ,proxy>200E_{\mu,\mathrm{proxy}}>200 TeV, i.e., with a Signalness S(Eμ)=ϕsignal(Eμ)ϕsignal(Eμ)+ϕbackground(Eμ)≳0.5S(E_\mu)=\frac{\phi_{\mathrm{signal}}(E_\mu)}{\phi_{\mathrm{signal}}(E_\mu)+\phi_{\mathrm{background}}(E_\mu)}\gtrsim0.5

IDMJDEnergy (TeV)SignalnessDecl. (deg)R.A. (deg)ΔEnergy\Delta_{\mathrm{Energy}} (TeV)ΔDecl.\Delta_{\mathrm{Decl.}} (deg)ΔR.A.\Delta_{\mathrm{R.A.}} (deg)
I.4a^{a}55370.7150………−110……
I.6b^{b}55421.5…………………
I.13a^{a}55722.4180………−30……
I.15a^{a}55896.9110………−190……
I.17a^{a}56063.0180………−20……
I.19a^{a}56211.8200………−10……
I.28a^{a}57049.5190………−20……
II.1IC5955056.74800.781.23−0.22+0.181.23^{+0.18}_{-0.22}29.51−0.38+0.4029.51^{+0.40}_{-0.38}………
II.2IC5955141.12500.5211.74−0.38+0.3211.74^{+0.32}_{-0.38}298.21−0.57+0.53298.21^{+0.53}_{-0.57}………
II.355355.53500.7222.63−2.15+2.7222.63^{+2.72}_{-2.15}346.29−2.11+2.33346.29^{+2.33}_{-2.11}10−0.951.36
II.455387.52000.5321.22−1.73+2.5321.22^{+2.53}_{-1.73}307.31−2.50+2.77307.31^{+2.77}_{-2.50}−300.220.35
II.555464.93900.7413.48−0.24+0.3213.48^{+0.32}_{-0.24}266.26−0.39+0.55266.26^{+0.55}_{-0.39}−700.08−0.03
II.655478.46000.8411.03−0.36+0.3011.03^{+0.30}_{-0.36}331.13−0.48+0.39331.13^{+0.39}_{-0.48}−60−0.060.05
II.755497.38500.880.34−0.14+0.160.34^{+0.16}_{-0.14}88.86−0.20+0.2988.86^{+0.29}_{-0.20}−100−0.16−0.09
II.855513.64600.793.17−0.27+0.313.17^{+0.31}_{-0.27}285.56−0.32+0.48285.56^{+0.48}_{-0.32}−600.02−0.39
II.955589.62100.551.16−0.29+0.201.16^{+0.20}_{-0.29}307.79−0.51+0.40307.79^{+0.40}_{-0.51}−300.130.08
II.1055702.84200.7620.07−0.93+0.9320.07^{+0.93}_{-0.93}234.93−0.98+0.98234.93^{+0.98}_{-0.98}120−0.23−0.20
II.1155764.22200.566.02−5.54+6.776.02^{+6.77}_{-5.54}315.53−4.10+2.47315.53^{+2.47}_{-4.10}100.73−0.13
II.1255911.36300.8418.49−1.43+1.8018.49^{+1.80}_{-1.43}36.56−1.54+1.2936.56^{+1.29}_{-1.54}−30−0.61−0.09
II.1356146.22500.591.60−0.25+0.231.60^{+0.23}_{-0.25}329.68−0.31+0.47329.68^{+0.47}_{-0.31}−100.03−0.42
II.1456226.67000.8728.16−0.60+0.4328.16^{+0.43}_{-0.60}169.98−0.71+0.96169.98^{+0.96}_{-0.71}−500.120.37
II.1556470.16600.8514.17−0.87+1.0814.17^{+1.08}_{-0.87}93.74−0.96+0.8393.74^{+0.83}_{-0.96}−10−0.290.36
II.1656521.84200.76−2.87−0.52+0.46-2.87^{+0.46}_{-0.52}223.77−0.45+0.42223.77^{+0.42}_{-0.45}201.57−1.12
II.1756579.92100.5510.28−0.46+0.3510.28^{+0.35}_{-0.46}32.92−0.57+0.7432.92^{+0.74}_{-0.57}−1800.08−0.02
II.1856666.58300.8833.02−0.42+0.4333.02^{+0.43}_{-0.42}293.12−0.99+0.67293.12^{+0.67}_{-0.99}−200.20−0.17
II.19c^{c}56757.12400.5981.22−5.86+7.7281.22^{+7.72}_{-5.86}2.11−47.87+185.552.11^{+185.55}_{-47.87}+160……
II.2056800.03000.6417.90−1.12+1.4717.90^{+1.47}_{-1.12}349.50−2.69+2.69349.50^{+2.69}_{-2.69}−100−0.150.11
II.2156817.63400.701.31−0.73+0.841.31^{+0.84}_{-0.73}106.26−1.72+2.20106.26^{+2.20}_{-1.72}00.02−0.00
II.2256819.244000.9911.45−0.15+0.1811.45^{+0.18}_{-0.15}110.65−0.58+0.46110.65^{+0.46}_{-0.58}−1000.030.02
II.2357157.92300.5912.14−0.43+0.4512.14^{+0.45}_{-0.43}91.49−0.63+0.8091.49^{+0.80}_{-0.63}−10−0.04−0.11
II.2457217.93000.6326.36−1.81+1.5826.36^{+1.58}_{-1.81}326.29−1.08+1.23326.29^{+1.23}_{-1.08}00.260.79
II.2557246.83700.736.17−0.42+0.456.17^{+0.45}_{-0.42}328.27−0.61+0.75328.27^{+0.75}_{-0.61}−100.17−0.13
II.2657269.82200.5728.08−0.42+0.4528.08^{+0.45}_{-0.42}133.77−0.42+0.72133.77^{+0.72}_{-0.42}00.08−0.23
II.2757312.72200.5819.95−1.94+2.3219.95^{+2.32}_{-1.94}197.53−2.08+2.05197.53^{+2.05}_{-2.08}−100.05−0.07
II.2857340.97300.8712.71−0.62+0.5612.71^{+0.56}_{-0.62}76.16−1.11+1.1576.16^{+1.15}_{-1.11}−100.11−0.14
II.2957478.63700.7315.48−0.64+0.5515.48^{+0.55}_{-0.64}151.22−0.53+0.47151.22^{+0.47}_{-0.53}−10−0.12−0.11
II.3057672.13800.741.16−2.60+5.371.16^{+5.37}_{-2.60}32.08−4.41+4.1732.08^{+4.17}_{-4.41}50−25.44d^{d}22.38d^{d}
II.31c^{c}57951.84400.7725.16−1.15+1.0025.16^{+1.00}_{-1.15}208.39−0.98+1.43208.39^{+1.43}_{-0.98}………
II.3258063.812000.937.44−0.24+0.277.44^{+0.27}_{-0.24}340.14−0.53+0.52340.14^{+0.52}_{-0.53}………
II.33e^{e}58141.72900.628.61−0.39+0.578.61^{+0.57}_{-0.39}76.33−1.89+1.8876.33^{+1.88}_{-1.89}………
II.34c^{c}58205.15500.827.44−1.27+1.377.44^{+1.37}_{-1.27}307.09−4.00+2.22307.09^{+2.22}_{-4.00}………
II.35e^{e}58264.33900.748.27−4.79+1.498.27^{+1.49}_{-4.79}210.67−6.35+5.56210.67^{+5.56}_{-6.35}………

Table 7.

In order to calculate the expected number of events at the detector, a large sample of events is simulated, taking into account the following propagation steps: the propagation of neutrinos through matter, primary neutrino–nucleon interactions (CSMS cross section: [30]), propagation of secondary particles through the ice (PROPOSAL code: [46]), and the production and propagation of Cherenkov photons in the detector medium (Chirkin 2015), CLSim code: [29]. They are modeled in detail before the final detector electronics and data acquisition are simulated. In addition, these simulations are repeated with altered assumptions on systematic uncertainties, e.g., scattering or absorption of the Cherenkov photons in the natural glacier and in the man-made boreholes.

Analysis Method

The measurement of the energy spectrum of astrophysical muon neutrinos is performed via a forward-folding fit: the experimentally observed events are compared to the simulation-based expectation in a two-dimensional histogram as a function of reconstructed zenith angle and reconstructed muon energy. By maximizing the Poisson likelihood from (2), the best-fitting flux hypothesis is obtained for the observed data DD, which consist of a number nn of observed events per bin:

L(D∣θ⃗,ξ⃗)=∏binNbinspPoisson(nbin,μbin(θ⃗,ξ⃗)).(2)\mathcal{L}(D|\vec{\theta},\vec{\xi}) = \prod_{\mathrm{bin}}^{N_{\mathrm{bins}}} p_{\mathrm{Poisson}}(n_{\mathrm{bin}},\mu_{\mathrm{bin}}(\vec{\theta},\vec{\xi})). \tag*{(2)}

The expected number of events at the detector, μ\mu, is modeled as a function of signal parameters θ⃗\vec{\theta} and nuisance parameters ξ⃗\vec{\xi}. The latter absorb systematic uncertainties of the measurement; see Section 3.1.

The two-dimensional analysis histogram consists of 50 energy bins equally spaced in $10^2 < \log_{10}(E_{\mu,\mathrm{proxy}}/\mathrm{GeV}) < 10^7 and 33 bins spaced in reconstructed cos⁡(Zenith)\cos(\mathrm{Zenith}) ranging from 85° to 180°. The binning is identical to the one used in [3] with Nbins=1650N_{\mathrm{bins}} = 1650.

Four flux components are considered and the expected number of events is given as a sum of these components per bin. First and second, atmospheric backgrounds are considered: conventional atmospheric neutrinos, emerging from the decay of pions and kaons in cosmic-ray-induced air showers and prompt atmospheric neutrinos from the decay of heavy charmed hadrons in the same showers. For both contributions, the prediction is updated with respect to [3], which used a prediction from [45]. The software Matrix Cascade Equation solver, MCEq ([35]), is employed assuming the flux model Gaisser-H4a ([38]) for primary cosmic rays and the hadronic interaction model Sibyll 2.3c ([36]). Third, a subdominant contamination from atmospheric muons is expected and modeled based on a large sample of cosmic-ray air showers. The resulting model for this background flux is denoted as the muon template. It is simulated with the CORSIKA package ([44]), assuming the mentioned primary cosmic-ray and hadronic interaction models. Fourth, and most important, a signal component of astrophysical neutrinos is considered. It represents the cumulative flux from all sources of high-energy astrophysical neutrinos, which we assume to be isotropic.

A wide range of parameterizations for the energy spectrum of the astrophysical component are investigated in Sections 4 and 5. The expected rate of events as a function of zenith and energy is visualized in Figure 1 for all four components assuming a single power-law spectrum of astrophysical neutrinos. Figure 2 further illustrates the potential to measure astrophysical neutrinos at high energies by showing the expected ratio of astrophysical signal over background.

Expected event distributions and data-to-model ratios for reconstructed zenith angle and muon energy proxy

Figure 1. Single power-law model: best-fit distributions and one-dimensional projection on reconstructed zenith angle and muon energy. The experimental data (black dots) are shown together with the best-fit expectation from simulations. Data taken in the IC59 detector configuration are kept in a separate analysis histogram. The conventional atmospheric component (purple) dominates the total flux for all zenith angles. Except for the highest energies, the line is thus hidden below the overall sum (black). The astrophysical component (red) is modeled as a single power law. The prompt component is drawn at nominal prediction for visualization (green-dashed) although a zero best-fit normalization is obtained. The best-fit expectation for the remaining background of muons is shown by the orange line. The central 68% range of the best-fit expectation is drawn as a gray band. It is obtained by varying all fit parameters according to their joint posterior distribution. The orange band additionally shows the statistical uncertainty of the simulated data.

Muon energy proxy versus cosine of the zenith angle, showing the ratio of signal over background

Figure 2. The color indicates the ratio of the signal (astrophysical) over the background (atmospheric and muonic) component expectation assuming the best-fit energy spectrum obtained in Section 4. The binning is equivalent to the two-dimensional histogram used in the analysis. For energies above Eμ,proxy≳50E_{\mu,\mathrm{proxy}} \gtrsim50 TeV, the number of observed events per bin is indicated, in the bins below this threshold; the number of events is proportional to the size of the black squares, with the maximum number of data events in a single bin being 3204.

Updated Treatment of Systematic Uncertainties

In order to absorb systematic uncertainties of the measurement, nuisance parameters are introduced in the fit. The first class of systematic uncertainties covers the detector response and reconstruction quality of IceCube. Optical properties of the natural ice and the ice in the man-made boreholes are considered as well as the optical detection efficiency. Compared to the last iteration of this analysis, the model of optical properties in the natural ice has been updated (Spice 3.2.1: [55]) and the impact of the ice in the boreholes (unified hole-ice model: [32]) is now taken into account as an additional source of uncertainty. For each of these systematic uncertainties, dedicated simulations have been performed and the impact on reconstructed energy and zenith have been parameterized.

The second class of systematic uncertainties arises from the flux predictions of atmospheric neutrinos. Since the absolute normalization of the primary cosmic-ray flux and the yield of neutrinos from cosmic-ray-induced air showers are not known precisely, the absolute flux normalizations ϕconv\phi_{\mathrm{conv}} and ϕprompt\phi_{\mathrm{prompt}} are free fit parameters. Additionally, the parameters ΔCR\Delta_{\mathrm{CR}} (CR spectral index shift) and λCR Model\lambda_{\mathrm{CR\ Model}}, which is a prior-constrained linear interpolation between the Gaisser-H4a and GST-4gen

Note. The astrophysical normalization is given in units of Cunits=10−18 GeV−1cm−2s−1sr−1C_{\mathrm{units}} = 10^{-18}\,\mathrm{GeV}^{-1}\mathrm{cm}^{-2}\mathrm{s}^{-1}\mathrm{sr}^{-1}. Confidence intervals (68%) are constructed from one-dimensional profile likelihood scans employing Wilks’ theorem [65].

flux models ([38]; [39]), are introduced analogously to Aartsen et al. ([3]) in order to cover uncertainties on the exact spectral shape of the primary cosmic-ray flux. This shape affects all atmospheric flux components alike. Furthermore, we update the treatment of atmospheric flux uncertainties by adopting the scheme from Barr et al. ([22]). By varying the relevant parameters HH (pions), WW (kaons), YY (kaons), and ZZ (kaons) within their estimated prior ranges and repeating the calculation of atmospheric fluxes, uncertainties from hadronic interaction models are covered by this suite of parameters. The parameters affecting the spectral shapes of the atmospheric fluxes are correlated with the absolute normalizations. For more details, see A.1 and Figure 6 in the supplementary material. Finally, the absolute normalization of the subdominant muon contamination is included as a parameter and allowed to vary within its estimated prior range as well, which is a Gaussian with a width of 0.5.

Pearson correlation coefficients between the signal and nuisance parameters

Figure 6. Pearson correlation coefficients between the signal and nuisance parameters are shown for the parameters of the single power-law fit.

Results: Single Power-law Model

The result of the likelihood fit is shown as a one-dimensional projection by the solid black line in Figure 1 together with the experimental data. At energies above Eμ,proxy≳100E_{\mu,\mathrm{proxy}} \gtrsim100 TeV, the excess of astrophysical neutrinos above the atmospheric components is clearly visible. Overall, the data are well described by the sum of atmospheric components and an astrophysical component following the standard paradigm of a single power-law energy spectrum. Figure 7 in the supplementary material shows the

Statistical pull per bin and pull density distribution

Figure 7. The upper figure shows the statistical pull per bin between the experimental data and the MC expectation assuming the best-fit energy spectrum obtained in sec:4. The lower figure shows the pull density distribution for the 1048 analysis bins containing data events.

statistical pull for all bins in the two-dimensional histogram indicating no obvious mismatches. Taking the systematic and statistical uncertainty of the best-fit expectation into account, a χ2\chi^2/ (degrees of freedom) for the single power-law fit is calculated to be 1.0, resulting in a pp-value of 50% and confirming that the fit result is a viable description of the measured data. The corresponding best-fit parameters of the astrophysical flux are listed in Table 2 and the profiled likelihood landscape of the two astrophysical signal parameters is shown in Figure 3. The sensitive energy range of the astrophysical measurement is determined by comparing the per-bin likelihood values of the best-fit hypothesis to the values obtained when repeating the fit assuming a background hypothesis. The true neutrino energy distribution is then weighted with these likelihood differences, and the central 90% range of the obtained distribution is Eν=15E_\nu= 15 TeV–5 PeV. This energy range extends to lower energies than previous measurements, where this energy range extended from Eν=200E_\nu= 200 TeV–8 PeV ([3]). This change is driven by the updated modeling of the conventional atmospheric flux in this energy region.

Profile likelihood landscape as a function of spectral index and astrophysical normalization

Figure 3. Single power-law model: profile likelihood landscape as a function of spectral index and astrophysical normalization. The best-fit parameters are marked as a white triangle. The turquoise square (Stettner & IceCube Collaboration [61], γ=2.28−0.05+0.08\gamma= 2.28^{+0.08}_{-0.05}, Φ=1.44−0.24+0.25\Phi= 1.44^{+0.25}_{-0.24}), the orange star (Haack & IceCube Collaboration [41], γ=2.19−0.10+0.10\gamma= 2.19^{+0.10}_{-0.10}, Φ=1.01−0.23+0.26\Phi= 1.01^{+0.26}_{-0.23}), and the pink circle (Aartsen et al. [3], γ=2.13−0.13+0.13\gamma= 2.13^{+0.13}_{-0.13}, Φ=0.90−0.27+0.30\Phi= 0.90^{+0.30}_{-0.27}) mark results of previous measurements.

Astrophysical Norm. ϕastro./Cunits\phi_{\mathrm{astro.}}/C_{\mathrm{units}}1.44−0.26+0.251.44^{+0.25}_{-0.26}
Spectral Index γSPL\gamma_{\mathrm{SPL}}2.37−0.09+0.092.37^{+0.09}_{-0.09}

Table 2. Single Power-law Model: Best-fit Parameters Assuming a Single Power-law Energy Spectrum.

Compared to the previous analysis by Aartsen et al. ([3]), a slightly softer spectral index of γSPL=2.37−0.09+0.09\gamma_{\mathrm{SPL}} = 2.37^{+0.09}_{-0.09} is obtained. Figure 3 shows the best-fit points of the previous measurements, and the changes between them are listed here as an overview. The measurements from Aartsen et al. ([3]) and Haack & IceCube Collaboration ([41]) are based on the same event selection and analysis method, with two years of additional data included in the latter. The changes between Haack & IceCube Collaboration ([41]) and the results from Stettner & IceCube Collaboration ([61]) are mostly driven by the updated atmospheric background prediction. Additionally, the Barr parameters are introduced, resulting in more fit freedom at medium energies, and the subdominant effects of neutrino oscillations are considered. The updates to the detector simulations (Pass-2) also occurred in between these analyses, but have a negligible effect on the fit result and uncertainty. Changes between the result reported in [61] and this analysis are a new event simulation including the effects of the hole-ice [32] and the addition of the subdominant muon component. A purely atmospheric hypothesis can be excluded with very high significance at a 5.6σ5.6\sigma level. Note that this significance is smaller than previously reported results because of the updated, more conservative treatment of systematic uncertainties as well as the observed softening of the spectral index.

Prompt Atmospheric Neutrinos

We find that the normalization of prompt atmospheric neutrinos is constrained less strongly compared to the last publication. This is related to the observed softening of the spectral index of the astrophysical neutrino flux, resulting in an overall more similar shape of the two flux components. For the calculation of a prompt limit, the dominating astrophysical flux in the regions below the sensitive energy range of this analysis poses a fundamental limitation to our ability to quantitatively constrain the subdominant flux of prompt neutrinos.

The best fit of the prompt normalization is still zero, independent of the different assumed astrophysical flux parameterizations discussed in Section 5. The prompt flux model prediction from MCEq (Sibyll 2.3c; H4a) is used here compared to the ERS prediction [33] in the last publication [3], these two models predict very similar fluxes. A more recent model [25] predicts a prompt atmospheric flux component with about a factor of three smaller normalization.

It has been checked in detail that the measurement of the astrophysical component is not impacted by the prompt flux normalization: for example, we find a spectral index of γSPL=2.33\gamma_{\mathrm{SPL}}=2.33, well within the quoted uncertainty range, if the likelihood fit is repeated with the prompt component fixed to its nominal prediction (fprompt=1.0f_{\mathrm{prompt}}=1.0). Also, the observed suppression of the prompt component occurs at similar strength when the likelihood fit is repeated while excluding energies above ≈15\approx15 TeV, confining it to an energy region dominated by conventional atmospheric neutrinos and well below the sensitive energy range for astrophysical neutrinos. The exact reason for this nonobservation remains an open question and an updated limit on the prompt flux normalization with respect to [3] is not computed here.

Results: Beyond the Single Power Law

Power-law energy spectra are well motivated from the assumed acceleration mechanisms of cosmic rays, but they extrapolate over large energy ranges and potential structures in the energy spectrum can thus not be identified. In this section, a number of parameterizations beyond a single power law are compared to the experimental data. Parameterizations with more than three signal parameters are, however, not considered because the statistics of observed events with high signalness is too low to constrain more fit parameters.

Power Law with Cutoff

The first natural extension to the single power law would be a cutoff in the energy spectrum, e.g., introduced if an astrophysical source of cosmic rays (and neutrinos) reaches its maximum energy. The flux parameterization given in (3) extends the single power law with an exponentially decaying term. The cutoff neutrino energy EcutoffE_{\mathrm{cutoff}} is consequently added as a third signal parameter in the likelihood fit:

Φastro.νμ+νˉμ(Eν)=ϕcutoff×(Eν100 TeV)−γcutoff×e−EνEcutoff.(3)\Phi_{\mathrm{astro.}}^{\nu_\mu+\bar{\nu}_\mu}(E_\nu)=\phi_{\mathrm{cutoff}}\times\left(\frac{E_\nu}{100\ \mathrm{TeV}}\right)^{-\gamma_{\mathrm{cutoff}}}\times e^{-\frac{E_\nu}{E_{\mathrm{cutoff}}}}. \tag*{(3)}

Table 3 lists the obtained fit result using this parameterization: while the astrophysical normalization does not change strongly, a hard spectral index of γcutoff=2.0−0.28+0.22\gamma_{\mathrm{cutoff}}=2.0^{+0.22}_{-0.28} and a cutoff energy of Ecutoff=1.25−0.56+1.72E_{\mathrm{cutoff}}=1.25^{+1.72}_{-0.56} PeV is found. Compared to the single power-law hypothesis, the fit improves by 2ΔLLH=4.242\Delta\mathrm{LLH}=4.24. The probability to randomly achieve any such improvement by introducing EcutoffE_{\mathrm{cutoff}} corresponds to a pp-value of p(>2ΔLLH∣SPL)=6.1%\mathrm{p}(>2\Delta\mathrm{LLH}|\mathrm{SPL})=6.1\%, which is calculated from pseudo-experiments obtained from Monte Carlo simulations.

ParameterBest-fit value
Astrophysical Norm. ϕcutoff/Cunits\phi_{\mathrm{cutoff}}/C_{\mathrm{units}}1.64−0.36+0.391.64^{+0.39}_{-0.36}
Spectral Index γcutoff\gamma_{\mathrm{cutoff}}2.0−0.28+0.222.0^{+0.22}_{-0.28}
Cutoff Energy Ecutoff/PeVE_{\mathrm{cutoff}}/\mathrm{PeV}1.25−0.56+1.721.25^{+1.72}_{-0.56}
Significance over SPL2ΔLLH=4.242\Delta\mathrm{LLH}=4.24
p(>2ΔLLH∣SPL)=6.1%\mathrm{p}(>2\Delta\mathrm{LLH}|\mathrm{SPL})=6.1\%

Table 3. Single Power Law with Cutoff: Best-fit Parameters. Note. Confidence intervals (68%) are constructed from one-dimensional profile likelihood scans employing Wilks’ theorem.

Log-parabola Model

Similarly to the cutoff hypothesis, the log-parabola model, widely used in gamma-ray astronomy, extends the single power law and allows for curvature of the spectrum. Its parameterization is given in (4) and Table 4 lists the obtained best-fit parameters of the likelihood fit with an astrophysical component following this model. Again, a hard spectral index of αLogParab.=2.03−0.31+0.22\alpha_{\mathrm{LogParab.}}=2.03^{+0.22}_{-0.31} is obtained, and the best fit of the curvature parameter is βLogParab.=0.45−0.22+0.29\beta_{\mathrm{LogParab.}}=0.45^{+0.29}_{-0.22}. Compared to the single power-law hypothesis, which corresponds to βLogParab.=0\beta_{\mathrm{LogParab.}}=0, the description of the experimental data is improved by 2ΔLLH=6.822\Delta\mathrm{LLH}=6.82. Analogously to the treatment for the cutoff hypothesis, this can be translated to a pp-value of p(>2ΔLLH∣SPL)=1.3%\mathrm{p}(>2\Delta\mathrm{LLH}|\mathrm{SPL})=1.3\%.

ParameterBest-fit value
Log-parabola Norm. ϕLogParab./Cunits\phi_{\mathrm{LogParab.}}/C_{\mathrm{units}}1.79−0.38+0.401.79^{+0.40}_{-0.38}
Spectral Index αLogParab.\alpha_{\mathrm{LogParab.}}2.03−0.31+0.222.03^{+0.22}_{-0.31}
Curvature parameter βLogParab.\beta_{\mathrm{LogParab.}}0.45−0.22+0.290.45^{+0.29}_{-0.22}
Significance over SPL2ΔLLH=6.822\Delta\mathrm{LLH}=6.82
p(>2ΔLLH∣SPL)=1.3%\mathrm{p}(>2\Delta\mathrm{LLH}|\mathrm{SPL})=1.3\%

Table 4. Log-parabola Model: Best-fit Parameters. Note. Confidence intervals (68%) are constructed from one-dimensional profile likelihood scans employing Wilks’ theorem.

Φastro.νμ+νˉμ(Eν)=ϕLogParab.×(Eν100 TeV)−αLogParab.−βLogParab.log⁡(Eν100 TeV)(4)\Phi_{\mathrm{astro.}}^{\nu_\mu+\bar{\nu}_\mu}(E_\nu)=\phi_{\mathrm{LogParab.}}\times\left(\frac{E_\nu}{100\ \mathrm{TeV}}\right)^{-\alpha_{\mathrm{LogParab.}}-\beta_{\mathrm{LogParab.}}\log\left(\frac{E_\nu}{100\ \mathrm{TeV}}\right)} \tag*{(4)}

Note. Note that all piece-wise normalizations are optimized simultaneously in the fit, i.e., correlations between the segments are fully taken into account. The given 68.27% uncertainty ranges are obtained from one-dimensional profile likelihood scans. Pieces 1 and 5 have been added to cover the full energy range here; upper limits (90% CL) are computed.

Piece-wise Parameterization

In order to overcome the limitations of the parameterizations discussed in previous sections, a “piece-wise” model is introduced. Here, the energy spectrum is described as sum of power laws with a fixed spectral index (γ=2.0\gamma= 2.0) in pre-defined, fixed segments of neutrino energy. This allows for measuring the flux strength in a well-defined range of neutrino energy and enables easy comparison to predictions from the literature and to other measurements. The total flux strength is then given by (5), where the flux normalizations per bin ϕpiecei\phi_{\mathrm{piece}}^i are fit parameters and ElowiE_{\mathrm{low}}^i and EhighiE_{\mathrm{high}}^i form the bounds of bin ii.

Φastro.ν+νˉ(Eν)=∑ipiecesχ(Eν)⋅ϕpiecei⋅(Eν100 TeV)−2.0\Phi_{\mathrm{astro.}}^{\nu+\bar{\nu}}(E_\nu)=\sum_i^{\mathrm{pieces}}\chi(E_\nu)\cdot\phi_{\mathrm{piece}}^i\cdot\left(\frac{E_\nu}{100\ \mathrm{TeV}}\right)^{-2.0}
χ(Eν)={1if Elowi<Eν<Ehighi0else(5)\chi(E_\nu)= \begin{cases} 1 & \text{if } E_{\mathrm{low}}^i < E_\nu< E_{\mathrm{high}}^i \\ 0 & \text{else} \end{cases} \tag*{(5)}

Prior to performing the fit on the experimental data, the energy ranges of the segments were defined to be equally spaced in log-energy spanning the sensitive energy range of the astrophysical measurement (see Section 4) with three segments. Additionally, one segment above and below have been added, respectively, to cover the full energy range. The full parameterization of the astrophysical flux is given in (5), and the energy ranges and obtained best-fit normalizations ϕpiecei\phi_{\mathrm{piece}}^i are listed in Table 5. Figure 4 visualizes the obtained flux measurement of the piece-wise parameterization together with the results of the single power law, power law with cutoff, and log-parabola models. In all models beyond the single power law, hints for a softening of the spectral shape as a function of energy are found.

Energy Range (EνE_\nu)Norm. (ϕpiecei/Cunits\phi_{\mathrm{piece}}^i/C_{\mathrm{units}})
Piece 1a^a100 GeV–15 TeV0.0−0.0+0.310.0^{+0.31}_{-0.0}
Piece 215 TeV–104 TeV2.22−0.8+0.082.22^{+0.08}_{-0.8}
Piece 3104 TeV–721 TeV1.21−0.31+0.321.21^{+0.32}_{-0.31}
Piece 4721 TeV–5 PeV0.33−0.18+0.220.33^{+0.22}_{-0.18}
Piece 5a^a5 PeV–100 PeV0.0−0.0+0.410.0^{+0.41}_{-0.0}

Table 5. Piece-wise Parameterization: Energy Ranges and Result of the Likelihood Fit

Summary of best-fit models for the astrophysical neutrino flux

Figure 4. Summary of best-fit models for the astrophysical neutrino flux. The bins from the piece-wise unfolding are marked in green and gray wherever only upper limits are calculated. The single power-law band is drawn in the sensitive energy range as defined in Section 4. All models with more degrees of freedom than the single power law show a hard spectrum at low and medium energies to a softer spectrum at highest energies.

Flux Predictions for Specific Source Classes

Besides the wide range of generic parameterizations for the energy spectrum discussed in the sections above, it is also possible to compare the experimental data to source-class specific flux predictions directly. The total astrophysical flux may originate from multiple source classes, thus it is not expected that a single flux prediction can fully explain the observed data. Instead, we model the total astrophysical component as sum of the predicted energy spectrum model times a free normalization ϕmodel\phi_{\mathrm{model}} and a single power law to cover other potential flux contributions:

Φastro.ν+νˉ(Eν)=ϕmodel×Model⁡(Eν)+ϕSPL×(Eν100 TeV)−γSPL(6)\Phi_{\mathrm{astro.}}^{\nu+\bar{\nu}}(E_\nu)=\phi_{\mathrm{model}}\times\operatorname{Model}(E_\nu) +\phi_{\mathrm{SPL}}\times\left(\frac{E_\nu}{100\ \mathrm{TeV}}\right)^{-\gamma_{\mathrm{SPL}}} \tag*{(6)}

A representative set of different source-class specific predictions have been selected, focusing on predictions not already covered by the performed test of a single power law, and including variations of the benchmark models shown in the publications (see Table 6). All these predictions model the cumulative expected flux at Earth for the given source class. The obtained fit results using these predictions are listed in Table 6. The test statistic TSfree modelTS_{\mathrm{free\ model}} from (7) compares the best-fit result including the additional component of the source-class specific flux prediction to the hypothesis of only a single power law. That is, TSfree model=0TS_{\mathrm{free\ model}}=0 implies that the description of the experimental data can not be improved with an additional contribution from the model prediction and the fit instead prefers the single power-law model. For these cases, upper limits on the model normalization are computed at 90% CL employing Wilks’ theorem.

ModelVariationTSfree modelTS_{\mathrm{free\ model}}ϕmodel\phi_{\mathrm{model}}ϕastro.SPL\phi_{\mathrm{astro.}}^{\mathrm{SPL}}γSPL\gamma_{\mathrm{SPL}}ULmodel90%UL_{\mathrm{model}}^{90\%}
Biehl et al. ([26]; GRB)Sum model A0.000.001.442.370.19
Sum model B0.000.001.442.374.92
Senno et al. ([57]; SFG w. HNe)Diffusion ∝E1/2\propto E^{1/2}−0.14-0.140.12−0.34+0.270.12^{+0.27}_{-0.34}1.122.40-
Diffusion ∝E1/3\propto E^{1/3}−0.41-0.410.23−0.37+0.290.23^{+0.29}_{-0.37}0.872.42-
Murase et al. ([52]; AGN inner Jets)Γ=2.0\Gamma=2.0, Blazar0.000.001.442.370.48
Γ=2.0\Gamma=2.0, Torus0.000.001.442.370.58
Γ=2.3\Gamma=2.3, Blazar0.000.001.442.370.48
Γ=2.3\Gamma=2.3, Torus0.000.001.442.370.27
Liu et al. ([49]; AGN winds)CR (Γ=2.1\Gamma=2.1)−0.98-0.980.87−0.88+0.490.87^{+0.49}_{-0.88}0.472.47-
CR (Γ=2.3\Gamma=2.3)−4.12-4.1214.3−5.97+3.6114.3^{+3.61}_{-5.97}0.122.04-
Padovani et al. ([53]; BL Lac)FνFγ=0.3\frac{F_{\nu}}{F_{\gamma}}=0.30.000.001.442.370.27
FνFγ=0.8\frac{F_{\nu}}{F_{\gamma}}=0.80.000.001.442.370.1
Kimura et al. ([47]; lowL AGN)Model B1−1.69-1.690.33−0.25+0.230.33^{+0.23}_{-0.25}0.892.46-
Model B40.000.001.442.370.24
Biehl et al. ([27]; TDE)No variations0.000.001.442.370.63
Tavecchio & Ghisellini ([64]; lowL BL Lac)No variations−1.74-1.740.32−0.24+0.220.32^{+0.22}_{-0.24}0.822.47-
Senno et al. ([58]; GRB w. choked Jets)No variations−4.36-4.362.6−1.16+0.42.6^{+0.4}_{-1.16}0.00--

Table 6. Results of the Likelihood Fits with an Additional Astrophysical Component following a Source-class Specific Flux Prediction

TSfree model=−2×log⁡(L(D∣ϕ^model,ϕ^SPL,γ^SPL,ξ^)L(D∣ϕmodel≡0.0,ϕSPL,γSPL,ξ⃗))(7)TS_{\mathrm{free\ model}} = -2 \times\log\left( \frac{\mathcal{L}(D|\hat{\phi}_{\mathrm{model}},\hat{\phi}_{\mathrm{SPL}},\hat{\gamma}_{\mathrm{SPL}},\hat{\xi})} {\mathcal{L}(D|\phi_{\mathrm{model}}\equiv0.0,\phi_{\mathrm{SPL}},\gamma_{\mathrm{SPL}},\vec{\xi})} \right) \tag*{(7)}

Similar to the results using the generic parameterizations for the energy spectrum, we find that model predictions with a softening of the spectral shape in the energy range 100 TeV – 1 PeV describe the experimental data better than the single power law. For example, the addition of spectral components as predicted by the models from [58] and [49] Note. The normalization is added as an additional fit parameter. A negative log-likelihood difference indicates that the data is better described with the additional component compared to the single power-law model. The last column shows the 90% CL (Wilks’ theorem) upper limit if the model normalization is fitted to zero. Tested source classes include different GRB and AGN scenarios, star-forming galaxies (SFG) with hypernovae (HNe), and TDEs. Low luminosity models are preceded by “lowL.”

(Γ=2.3)(\Gamma=2.3) are favored compared to the pure single power-law hypothesis, with test statistics reaching −TSfree model>4-TS_{\mathrm{free\ model}}>4. For these two models, the best-fit model normalizations are multiples of the model prediction, substantially reducing the strength of the single power-law component, while neither model is claiming to account for the entire diffuse emission. When fixing the model normalization (Φmodel=1.0\Phi_{\mathrm{model}}=1.0 in (6)) and comparing to the single power-law model, only the models from [49] and [58] result in negative test statistics, indicating a small but insignificant preference of those model predictions [62]. For every other model in Table 6 which is yielding negative test statistics when allowing a free model normalization, the fitted normalization is smaller than the nominal prediction in the original publication.

Models in Table 6 predicting some spectral hardening in the considered energy range, like [27], [52] and [53] are mildly disfavored. While this does not disfavor them as source models for individual neutrino sources, they are less likely to be main contributors to the overall diffuse astrophysical neutrino flux.

Discussion and Outlook

We have presented an updated measurement of the astrophysical muon–neutrino flux from the northern celestial sky. The last measurement from [3] observed a flux compatible with a single power law described by a spectral index of γ=2.13−0.13+0.13\gamma=2.13^{+0.13}_{-0.13}. Our update consists of more than three years of additional data (roughly doubling statistics), a re-processing of the full data with latest calibration and filtering standards (Pass-2) and an updated treatment of systematic uncertainties (detector effects and atmospheric fluxes). Assuming a single power-law energy spectrum, we find a spectral index of 2.37−0.09+0.092.37^{+0.09}_{-0.09}. This is in agreement with but slightly softer than earlier iterations of this analysis. This change is partly caused by the updated atmospheric flux models and uncertainty treatment but also by updated detector simulations of photon detection as well as the added data. In addition, we tested parameterizations beyond the single power law for the first time and find hints for a softening of the spectral shape as a function of energy at the two-sigma confidence level (cutoff, log-parabola, and piece-wise models). Figure 5 shows the result of the spectral fit in comparison to other measurements of the diffuse astrophysical flux by IceCube ([12], [7]; [19]) and the mild excess over expected atmospheric backgrounds observed by ANTARES ([37]). In the energy range where all referred analyses are sensitive, the observed event rates agree well with each other. The analyses themselves are based on event samples of varying statistical sizes, covering different energies and neutrino flavors. The advantages of the different analysis methods and the relations between the spectral results shown in Figure 5 are discussed in detail in [19]. With continued data taking of IceCube, it is expected that these measurements can be further improved in the future. Furthermore, with the future IceCube-Gen2 Observatory ([2]), we expect a substantial increase of exposure by a factor of ≈6\approx6, which will improve the statistics in the here probed energy range and also allow for extensions of the energy range to higher energies. A further goal of the collaboration is combining the measured astrophysical fluxes from different detection channels into a single consistent analysis based on a global fit of the data.

Summary of astrophysical neutrino flux measurements and uncertainty contours

Figure 5. Summary of astrophysical neutrino flux measurements. Best-fit parameters and uncertainty contours for the single power-law hypothesis are drawn for studies based on high-energy starting events [19], cascade-like events [13], and an inelasticity study [12] by IceCube. ANTARES observes a mild excess of events over the expected atmospheric backgrounds in a combined study of tracks and cascades [37].

The IceCube collaboration acknowledges the significant contributions to this manuscript from Philipp Fürst, Jöran Stettner, and Christopher Wiebusch. We acknowledge the support from the following agencies: USA—U.S. National Science Foundation-Office of Polar Programs, U.S. National Science Foundation-Physics Division, U.S. National Science Foundation-EPSCoR, Wisconsin Alumni Research Foundation, Center for High Throughput Computing (CHTC) at the

Notes. The reconstructed muon energy is as used in the analysis [17]. The directional reconstruction is based on a more sophisticated reconstruction algorithm that is also applied to real-time alerts [10, 11]. The given statistical uncertainty ranges are 90% CL, derived from reconstructions performed on a sample of similar events [54]. See text for a description of the Pass-2 re-calibration campaign that leads to changes of reconstructed energy and direction for some events. Seven events that have been reported in the last publications [3, 41] are marked as dropped, either because their reconstructed energy falls below threshold or because they do not pass the event selection anymore.

University of Wisconsin-Madison, Open Science Grid (OSG), Extreme Science and Engineering Discovery Environment (XSEDE), Frontera computing project at the Texas Advanced Computing Center, U.S. Department of Energy-National Energy Research Scientific Computing Center, Particle astrophysics research computing center at the University of Maryland, Institute for Cyber-Enabled Research at Michigan State University, and Astroparticle physics computational facility at Marquette University; Belgium—Funds for Scientific Research (FRS-FNRS and FWO), FWO Odysseus and Big Science programmes, and Belgian Federal Science Policy Office (Belspo); Germany—Bundesministerium für Bildung und Forschung (BMBF), Deutsche Forschungsgemeinschaft (DFG), Helmholtz Alliance for Astroparticle Physics (HAP), Initiative and Networking Fund of the Helmholtz Association, Deutsches Elektronen Synchrotron (DESY), and High Performance Computing cluster of the RWTH Aachen; Sweden—Swedish Research Council, Swedish Polar Research Secretariat, Swedish National Infrastructure for Computing (SNIC), and Knut and Alice Wallenberg Foundation; Australia—Australian Research Council; Canada—Natural Sciences and Engineering Research Council of Canada, Calcul Québec, Compute Ontario, Canada Foundation for Innovation, West-Grid, and Compute Canada; Denmark—Villum Fonden and Carlsberg Foundation; New Zealand—Marsden Fund; Japan—Japan Society for Promotion of Science (JSPS) and Institute for Global Prominent Research (IGPR) of Chiba University; Korea—National Research Foundation of Korea (NRF); Switzerland—Swiss National Science Foundation (SNSF); United Kingdom—Department of Physics, University of Oxford.

Appendix ASupplementary Material

Barr Treatment of Atmospheric Uncertainties

The nuisance parameters are included in the fit to cover the systematic uncertainties affecting this measurement, with the goal of measuring an unbiased result. The scheme from Barr et al. [22] was adopted to cover atmospheric flux uncertainties. Previous analyses used a parameter describing the ratio between the integrated neutrino fluxes arising from kaon and pion decays, respectively [3]; [41]. In principle, the Barr scheme allows for an uncertainty in the production yield of each individual meson, for example pions and antipions. Different from a global scaling of these production yields, each parameter describes uncertainties in a specific region of meson production phase space, with the goal of having different parameters for regions dominated by different physical effects and with different experimental coverage. Since the total ν+νˉ\nu+ \bar{\nu} flux is measured in this analysis, the parameters for mesons and antimesons can be combined into single parameters here. Flux gradients are then calculated from flux predictions obtained with different parameter values, and the Barr parameters in the fit then scale this gradient to obtain a flux prediction depending on parameter value. Since they affect the neutrino production, the Barr parameters are correlated to the absolute normalization of the conventional atmospheric flux, but crucially also introduce energy-dependent flux variations [62]. The correlations between the nuisance parameters are shown in Figure 6.

Correlation Coefficients

The Pearson correlation coefficients between signal and nuisance parameters are calculated and shown in Figure 6

Pull Density Distribution

The per-bin difference between experimental data and MC expectation is shown in Figure 7, which shows the statistical pull per bin in the upper part and the respective pull distribution in the lower part.

ORCID iDs

R. Abbasi https://orcid.org/0000-0001-6141-4205

M. Ackermann https://orcid.org/0000-0001-8952-588X

J. A. Aguilar https://orcid.org/0000-0003-2252-9514

J. Ahlers https://orcid.org/0000-0003-0709-5631

J. M. Alameddine https://orcid.org/0000-0002-9534-9189

G. Anton https://orcid.org/0000-0003-2039-4724

C. Argüelles https://orcid.org/0000-0003-4186-4182

A. Balagopal V. https://orcid.org/0000-0001-5367-8876

A. Barbano https://orcid.org/0000-0002-4836-7093

S. W. Barwick https://orcid.org/0000-0003-2050-6714

V. Basu https://orcid.org/0000-0002-9528-2009

S. Baur https://orcid.org/0000-0002-3329-1276

J. J. Beatty https://orcid.org/0000-0003-0481-4952

K.-H. Becker https://orcid.org/0000-0002-1748-7367

S. BenZvi https://orcid.org/0000-0001-5537-4710

E. Bernardini https://orcid.org/0000-0003-3108-1141

E. Blaufuss https://orcid.org/0000-0001-5450-1757

S. Blot https://orcid.org/0000-0003-1089-3001

S. Böser https://orcid.org/0000-0002-5918-4890

O. Botner https://orcid.org/0000-0001-8588-7306

F. Bradascio https://orcid.org/0000-0002-7750-5256

A. Burgman https://orcid.org/0000-0003-1276-676X

M. A. Campana https://orcid.org/0000-0003-4162-5739

C. Chen https://orcid.org/0000-0002-8139-4106

D. Chirkin https://orcid.org/0000-0003-4911-1345

B. A. Clark https://orcid.org/0000-0003-4089-2245

K. Clark https://orcid.org/0000-0003-2467-6825

A. Coleman https://orcid.org/0000-0003-1510-1712

J. M. Conrad https://orcid.org/0000-0002-6393-0438

P. Coppin https://orcid.org/0000-0001-6869-1280

P. Correa https://orcid.org/0000-0002-1158-6735

R. Cross https://orcid.org/0000-0003-0081-8024

P. Dave https://orcid.org/0000-0002-3879-5115

C. De Clercq https://orcid.org/0000-0001-5266-7059

J. J. DeLaunay https://orcid.org/0000-0001-5229-1995

D. Delgado López https://orcid.org/0000-0002-4306-8828

H. Dembinski https://orcid.org/0000-0003-3337-3850

A. Desai https://orcid.org/0000-0001-7405-9994

P. Desiati https://orcid.org/0000-0001-9768-1858

K. D. de Vries https://orcid.org/0000-0002-9842-4068

G. de Wasseige https://orcid.org/0000-0002-1010-5100

T. DeYoung https://orcid.org/0000-0003-4873-3783

A. Diaz https://orcid.org/0000-0001-7206-8336

J. C. Díaz-Vélez https://orcid.org/0000-0002-0087-0693

H. Dujmovic https://orcid.org/0000-0003-1891-0718

M. A. DuVernois https://orcid.org/0000-0000-2987-9691

P. Eller https://orcid.org/0000-0001-6354-5209

A. R. Fazely https://orcid.org/0000-0002-6907-8020

C. Finley https://orcid.org/0000-0003-3509-390X

D. Fox https://orcid.org/0000-0002-3714-672X

A. Franckowiak https://orcid.org/0000-0002-5605-2219

P. Fürst https://orcid.org/0000-0002-7951-8042

T. K. Gaisser https://orcid.org/0000-0003-4717-6620

E. Ganster https://orcid.org/0000-0003-4393-6944

A. Garcia https://orcid.org/0000-0002-8186-2459

S. Garrappa https://orcid.org/0000-0003-2403-4582

A. Ghadimi https://orcid.org/0000-0002-6350-6485

T. Glauch https://orcid.org/0000-0003-1804-4055

T. Glüsenkamp https://orcid.org/0000-0002-2268-9297

S. Griswold https://orcid.org/0000-0002-7321-7513

P. Gutjahr https://orcid.org/0000-0001-7980-7285

A. Hallgren https://orcid.org/0000-0001-7751-4489

L. Halve https://orcid.org/0000-0003-2237-6714

F. Halzen https://orcid.org/0000-0001-6224-2417

A. Haungs https://orcid.org/0000-0002-9638-7574 K. Helbing https://orcid.org/0000-0003-2072-4172

F. Henningsen https://orcid.org/0000-0002-0680-6588

C. Hill https://orcid.org/0000-0003-0647-9174

F. Huang https://orcid.org/0000-0002-6014-5928

T. Huber https://orcid.org/0000-0002-6515-1673

N. Iovine https://orcid.org/0000-0001-7965-2252

G. S. Japaridze https://orcid.org/0000-0002-7000-5291

B. J. P. Jones https://orcid.org/0000-0003-3400-8986

D. Kang https://orcid.org/0000-0002-5149-9767

W. Kang https://orcid.org/0000-0003-3980-3778

A. Kappes https://orcid.org/0000-0003-1315-3711

T. Karg https://orcid.org/0000-0003-3251-2126

M. Karl https://orcid.org/0000-0003-2475-8951

A. Karle https://orcid.org/0000-0001-9889-5161

U. Katz https://orcid.org/0000-0002-7063-4418

M. Kauer https://orcid.org/0000-0003-1830-9076

J. L. Kelley https://orcid.org/0000-0002-0846-4542

A. Kheirandish https://orcid.org/0000-0001-7074-0539

S. R. Klein https://orcid.org/0000-0003-2841-6553

R. Koirala https://orcid.org/0000-0002-7735-7169

H. Kolanoski https://orcid.org/0000-0003-0435-2524

C. Kopper https://orcid.org/0000-0001-6288-7637

D. J. Koskinen https://orcid.org/0000-0002-0514-5917

P. Koundal https://orcid.org/0000-0002-5917-5230

M. Kovacevich https://orcid.org/0000-0002-5019-5745

M. Kowalski https://orcid.org/0000-0001-8594-8666

N. Kurahashi https://orcid.org/0000-0003-1047-8094

C. Lagunas Gualda https://orcid.org/0000-0002-9040-7191

M. J. Larson https://orcid.org/0000-0002-6996-1155

F. Lauber https://orcid.org/0000-0001-5648-5930

J. P. Lazar https://orcid.org/0000-0003-0928-5025

K. Leonard https://orcid.org/0000-0002-8795-0601

A. Leszczyńska https://orcid.org/0000-0003-0935-6313

Q. R. Liu https://orcid.org/0000-0003-3379-6423

L. Lu https://orcid.org/0000-0003-3175-7770

F. Lucarelli https://orcid.org/0000-0002-9558-8788

A. Ludwig https://orcid.org/0000-0001-9038-4375

W. Luszczak https://orcid.org/0000-0003-3085-0674

Y. Lyu https://orcid.org/0000-0002-2333-4383

W. Y. Ma https://orcid.org/0000-0003-1251-5493

J. Madsen https://orcid.org/0000-0003-2415-9959

I. C. Mariş https://orcid.org/0000-0002-5771-1124

R. Maruyama https://orcid.org/0000-0003-2794-512X

F. McNally https://orcid.org/0000-0002-0785-2244

K. Meagher https://orcid.org/0000-0003-3967-1533

M. Meier https://orcid.org/0000-0002-9483-9450

S. Meighen-Berger https://orcid.org/0000-0001-6579-2000

T. Montaruli https://orcid.org/0000-0001-5014-2152

R. W. Moore https://orcid.org/0000-0003-4160-4700

M. Moulai https://orcid.org/0000-0001-7909-5812

R. Naab https://orcid.org/0000-0003-2512-466X

R. Nagai https://orcid.org/0000-0001-7503-2777

J. Necker https://orcid.org/0000-0003-0280-7484

H. Niederhausen https://orcid.org/0000-0002-9566-4904

M. U. Nisa https://orcid.org/0000-0002-6839-3944

A. Obertacke Pollmann https://orcid.org/0000-0002-2492-043X

B. Oeyen https://orcid.org/0000-0003-2940-3164

E. O’Sullivan https://orcid.org/0000-0003-1882-8802

H. Pandya https://orcid.org/0000-0002-6138-4808

N. Park https://orcid.org/0000-0002-4282-736X

E. N. Paudel https://orcid.org/0000-0001-9276-7994

C. Pérez de los Heros https://orcid.org/0000-0002-2084-5866

A. Pizzuto https://orcid.org/0000-0002-8466-8168

M. Plum https://orcid.org/0000-0001-8691-242X

A. Porcelli https://orcid.org/0000-0002-3220-6295

C. Raab https://orcid.org/0000-0001-9921-2668

M. Rameez https://orcid.org/0000-0001-5023-5631

A. Rehman https://orcid.org/0000-0001-7616-5790

R. Reimann https://orcid.org/0000-0002-1983-8271

E. Resconi https://orcid.org/0000-0003-0705-2770

W. Rhode https://orcid.org/0000-0003-2636-5000

B. Riedel https://orcid.org/0000-0002-9524-8943

M. Rongen https://orcid.org/0000-0002-7057-1007

C. Rott https://orcid.org/0000-0002-6958-6033

D. Rysewyk Cantu https://orcid.org/0000-0002-3612-6129

I. Safa https://orcid.org/0000-0001-8737-6825

A. Sandrock https://orcid.org/0000-0000-0002-6779-1172

J. Sandroos https://orcid.org/0000-0002-0629-0630

M. Santander https://orcid.org/0000-0001-7297-8217

S. Sarkar https://orcid.org/0000-0002-3542-858X

S. Sarkar https://orcid.org/0000-0002-1206-4330

K. Satalecka https://orcid.org/0000-0002-7669-266X

A. Schneider https://orcid.org/0000-0002-0895-3477

J. Schneider https://orcid.org/0000-0001-7752-5700

F. G. Schröder https://orcid.org/0000-0001-8495-7210

S. Sclafani https://orcid.org/0000-0001-9446-1219

M. Silva https://orcid.org/0000-0001-6940-8184

B. Smithers https://orcid.org/0000-0003-1273-985X

G. M. Spiczak https://orcid.org/0000-0002-0030-0519

C. Spiering https://orcid.org/0000-0001-7372-0074

R. Stein https://orcid.org/0000-0003-2434-0387

J. Stettner https://orcid.org/0000-0003-1042-3675

T. Stezelberger https://orcid.org/0000-0003-2676-9574

T. Stuttard https://orcid.org/0000-0001-7944-279X

G. W. Sullivan https://orcid.org/0000-0002-2585-2352

I. Taboada https://orcid.org/0000-0003-3509-3457

S. Ter-Antonyan https://orcid.org/0000-0002-5788-1369

K. Tollefson https://orcid.org/0000-0001-9725-1479

S. Toscano https://orcid.org/0000-0002-1860-2240

C. F. Tung https://orcid.org/0000-0001-6920-7841

C. F. Turley https://orcid.org/0000-0002-9689-8075

M. A. Unland Elorrieta https://orcid.org/0000-0002-6124-3255

J. Vandenbroucke https://orcid.org/0000-0002-9867-6548

N. van Eijndhoven https://orcid.org/0000-0001-5558-3328

J. van Santen https://orcid.org/0000-0002-2412-9728

S. Verpoest https://orcid.org/0000-0002-3031-3206

T. B. Watson https://orcid.org/0000-0002-8631-2253

C. Weaver https://orcid.org/0000-0003-2385-2559

C. Wendt https://orcid.org/0000-0001-8076-8877

N. Whitehorn https://orcid.org/0000-0002-3157-0407

C. H. Wiebusch https://orcid.org/0000-0002-6418-3008

M. Wolf https://orcid.org/0000-0001-9991-3923

S. Yoshida https://orcid.org/0000-0003-2480-5105

T. Yuan https://orcid.org/0000-0001-5710-508X

References

  1. [1]Aartsen, M., Abbasi, R., Ackermann, M., et al. 2021a, Natur, 591, 220
  2. [2]Aartsen, M., Abbasi, R., Ackermann, M., et al. 2021b, JPhG, 48, 060501
  3. [3]Aartsen, M., Abraham, K., Ackermann, M., et al. 2016, ApJ, 833, 3
  4. [6]Aartsen, M., Ackermann, M., Adams, J., et al. 2018a, Sci, 361, 147
  5. [7]Aartsen, M., Ackermann, M., Adams, J., et al. 2020a, PhRvL, 124, 051103
  6. [8]Aartsen, M., Ackermann, M., Adams, J., et al. 2020b, JInst, 15, P06032
  7. [9]Aartsen, M. G., Abbasi, R., Abdou, Y., et al. 2013, PhRvL, 111, 021103
  8. [10]Aartsen, M. G., Ackermann, M., Adams, J., et al. 2017d, ApJ, 843, 112
  9. [11]Aartsen, M. G., Ackermann, M., Adams, J., et al. 2018b, Sci, 361, 6398
  10. [12]Aartsen, M. G., Ackermann, M., Adams, J., et al. 2019, PhRvD, 99, 032004
  11. [13]Aartsen, M. G., Ackermann, M., Adams, J., et al. 2020c, PhRvL, 125, 121104
  12. [14]Abbasi, R., Abdou, Y., Abu-Zayyad, T., et al. 2010, NIMPA, 618, 139
  13. [15]Abbasi, R., Abdou, Y., Abu-Zayyad, T., et al. 2012a, Natur, 484, 351
  14. [16]Abbasi, R., Abdou, Y., Abu-Zayyad, T., et al. 2012b, APh, 35, 615
  15. [17]Abbasi, R., Ackermann, M., et al. 2013, NIMPA, 703, 190
  16. [18]Abbasi, R., Ackermann, M., Adams, J., et al. 2009, Natur, 460, 294
  17. [19]Abbasi, R., Ackermann, M., Adams, J., et al. 2021, PhRvD, 104, 022002
  18. [20]Ahlers, M., & Halzen, F. 2018, PrPNP, 102, 73
  19. [21]Ahrens, J., Bai, X., Bay, R., et al. 2004, NIMPA, 524, 169
  20. [22]Barr, G. D., Robbins, S., Gaisser, T. K., & Stanev, T. 2006, PhRvD, 74, 094009
  21. [23]Becker, J. K. 2008, PhR, 458, 173
  22. [24]Becker Tjus, J., & Merten, L. 2020, PhR, 872, 1
  23. [25]Bhattacharya, A., Enberg, R., Reno, M. H., Sarcevic, I., & Stasto, A. 2015, JHEP, 2015, 110
  24. [26]Biehl, D., Boncioli, D., Fedynitch, A., & Winter, W. 2018, A&A, 611, A101
  25. [27]Biehl, D., Boncioli, D., Lunardini, C., & Winter, W. 2018, NatSR, 8, 10828
  26. [29]Chirkin, D., Díaz-Vélez, J. C., Kopper, C., et al. 2019, 2019 15th International Conference on eScience (eScience) (San Diego, CA: IEEE)
  27. [30]Cooper-Sarkar, A., Mertsch, P., & Sarkar, S. 2011, JHEP, 2011, 42
  28. [31]Dai, L., & Fang, K. 2017, MNRAS, 469, 1354
  29. [32]Eller, P. 2019, Unified Hole-ice Model: Angular-acceptance Code, GitHub, https://github.com/philippeller/angular_acceptance
  30. [33]Engberg, R., Reno, M. H., & Sarcevic, I. 2008, PhRvD, 78, 043005
  31. [34]Farrar, G. R., & Piran, T. 2014, arXiv:1411.0704arxiv.org/abs/1411.0704
  32. [35]Fedynitch, A., Engel, R., Gaisser, T. K., Riehn, F., & Stanev, T. 2015, in ISVHECRI 2014 – 18th Int. Symp. on Very High Energy Cosmic Ray Interactions, 99, ed. D. Berge et al., 08001
  33. [36]Fedynitch, A., Riehn, F., Engel, R., Gaisser, T. K., & Stanev, T. 2019, PhRvD, 100, 103018
  34. [37]Fusco, L. A., & Versari, F. 2019, ICRC (Madison, WI), 358, 891
  35. [38]Gaisser, T. K. 2012, APh, 35, 801
  36. [39]Gaisser, T. K., Stanev, T., & Tilav, S. 2013, FrPhy, 8, 748
  37. [40]Guépin, C., & Kotera, K. 2017, A&A, 603, A76
  38. [41]Haack, C., Wiebusch, C. & (IceCube Collaboration) 2018, ICRC (Busan), 301, 1005
  39. [42]Halzen, F., & Klein, S. R. 2010, RSci, 81, 081101
  40. [43]Hayasaki, K. 2021, NatAs, 5, 436
  41. [44]Heck, D., Knapp, J., Capdevielle, J., et al. 1998, Karlsruhe, Forschungszentrum, Technical Report, FZKA-6019
  42. [45]Honda, M., Kajita, T., Kasahara, K., Midorikawa, S., & Sanuki, T. 2007, PhRvD, 75, 043006
  43. [46]Köhne, J. H., Frantzen, K., Schmitz, M., et al. 2013, CoPhC, 184, 2070
  44. [47]Kimura, S. S., Murase, K., & Toma, K. 2015, ApJ, 806, 159
  45. [48]Learned, J. G., & Mannheim, K. 2000, ARNPS, 50, 679
  46. [49]Liu, R.-Y., Murase, K., Inoue, S., Ge, C., & Wang, X.-Y. 2018, ApJ, 858, 9
  47. [50]Loeb, A., & Waxman, E. 2006, JCAP, 05, 003
  48. [51]Lunardini, C., & Winter, W. 2017, PhRvD, 95, 123001
  49. [52]Murase, K., Inoue, Y., & Dermer, C. D. 2014, PhRvD, 90, 023007
  50. [53]Padovani, P., Petropoulou, M., Giommi, P., & Resconi, E. 2015, MNRAS, 452, 1877
  51. [54]Rädel, L. 2017, PhD thesis, RWTH Aachen Univ., http://publications.rwth-aachen.de/record/709576
  52. [55]Rongen, M. 2019, PhD thesis, RWTH Aachen Univ., http://publications.rwth-aachen.de/record/771097
  53. [56]Schatto, K. 2014, PhD thesis, Johannes Gutenberg-Universität Mainz, https://opscience.ub.uni-mainz.de/handle/20.500.12030/2899
  54. [57]Senno, N., Mészáros, P., Murase, K., Baerwald, P., & Rees, M. J. 2015, ApJ, 806, 24
  55. [58]Senno, N., Murase, K., & Meszaros, P. 2016, PhRvD, 93, 083003
  56. [59]Senno, N., Murase, K., & Meszaros, P. 2017, ApJ, 838, 3
  57. [60]Stein, R., Velzen, S. v., Kowalski, M., et al. 2021, NatAs, 5, 510
  58. [61]Stettner, J. & IceCube Collaboration 2019, ICRC (Madison, WI), 36, 1017
  59. [62]Stettner, J. B. 2021, PhD thesis, RWTH Aachen Univ. https://publications.rwth-aachen.de/record/811376
  60. [63]Tamborra, I., Ando, S., & Murase, K. 2014, JCAP, 09, 043
  61. [64]Tavecchio, F., & Ghisellini, G. 2015, MNRAS, 451, 1502
  62. [65]Wilks, S. S. 1938, Ann. Math. Statist., 9, 60

Paper details

Contents