Measurement of the high-energy all-flavor neutrino-nucleon cross section with IceCube

Abbasi, R.; Ackermann, M.; Adams, J.; Aguilar, J. A.; Ahlers, M.; Ahrens, M.; Alispach, C.; Alves, A. A.; Amin, N. M.; Andeen, K.; Anderson, T.; Ansseau, I; Anton, G.; Arguelles, C.; Axani, S.; Bai, X.; Balagopal, A.; Barbano, A.; Barwick, S. W.; Bastian, B.; Basu, Alakananda; Baum, Oliver; Baur, S.; Bay, R.; Beatty, J. J.; Becker, K-H; Tjus, J. Becker; Bellenghi, C.; BenZvi, S.; Berley, D.; Bernardini, E.; Besson, D. Z.; Binder, G.; Bindig, D.; Blaufuss, E.; Blot, S.; Boeser, S.; Botner, O.; Boettcher, J.; Bourbeau, E.; Bourbeau, J.; Bradascio, F.; Braun, J.; Bron, S.; Brostean-Kaiser, J.; Burgman, A.; Busse, R. S.; Koskinen, D. J.; Rameez, M.; Stuttard, T.; Icecube Collaboration

The flux of high-energy neutrinos passing through the Earth is attenuated due to their interactions with matter. The interaction rate is determined by the neutrino interaction cross section and affects the flux arriving at the IceCube Neutrino Observatory, a cubic-kilometer neutrino detector embedded in the Antarctic ice sheet. We present a measurement of the neutrino cross section between 60 TeV and 10 PeV using the high-energy starting event (HESE) sample from IceCube with 7.5 years of data. The result is binned in neutrino energy and obtained using both Bayesian and frequentist statistics. We find it compatible with predictions from the Standard Model. While the cross section is expected to be flavor independent above 1 TeV, additional constraints on the measurement are included through updated experimental particle identification (PID) classifiers, proxies for the three neutrino flavors. This is the first such measurement to use a ternary PID observable and the first to account for neutrinos from tau decay.

I. INTRODUCTION
In the Standard Model (SM), neutrino interactions are mediated by W AE and Z 0 bosons for charged-current (CC) and neutral-current (NC) channels, respectively. At energies above a few GeV, the dominant process is deep inelastic scattering (DIS) off of individual partons within the nucleon. Calculations in the perturbative QCD (pQCD) formalism rely on parton distribution functions (PDFs) obtained mostly from DIS experiments [1][2][3]. Uncertainties on the PDFs lead to uncertainties on the cross section. An alternative approach [4] based on an empirical color dipole model of the nucleon along with the assumption that all cross sections increase at high energies as ln 2 s results in good agreement with the latest pQCD calculations. Proposed extensions of the SM based on large extra dimensions opening up above the Fermi scale predict a sharp rise in the neutrino-nucleon cross section above the SM value. One such model [5], which was motivated by the claimed detection of cosmic rays above the Greisen-Zatsepin-Kuzmin (GZK) bound, assumes that neutrinonucleon interaction is mediated by a massive spin-2 boson. This allows the neutrino-nucleon cross section to climb above 10 −27 cm 2 at E ν > 10 19 eV. Another possibility if spacetime has greater than four dimensions allows for the production of microscopic black holes in high-energy particle interactions and also leads to an increased neutrino-nucleon cross section above approximately 1 PeV [6]. Such scenarios where the cross section increases steeply with energy could also be due to the existence of exotic particles such as leptoquarks [7] or sphalerons [8], both of which have been discussed in the context of neutrino telescopes and could be probed via measurements of the high-energy neutrino cross section.
At energies above 40 TeV, the Earth becomes opaque to neutrinos. For a power-law spectrum proportional to E −γ at Earth's surface, the ratio of the flux arriving at IceCube to that at Earth's surface, ΦðE ν Þ=Φ 0 ðE ν Þ, depends on the Earth column density, neutrino energy, E ν , spectral index γ (through secondaries), and neutrino cross section. The Earth column density is defined as tðθÞ ¼ R y max 0 ρðy;θÞdy, where θ is the arrival direction of the neutrino, y max is its path length through the Earth, and ρðy; θÞ is the density at a point y along the path. Figure 1 shows the electron neutrino and antineutrino ΦðE ν Þ=Φ 0 ðE ν Þ assuming a surface flux with γ ¼ 2, arising from Fermi acceleration at shocks [9,10]. The spectral index affects the arrival flux through secondaries produced by tau decay in CC interactions or NC interactions. In ν e and ν μ CC interactions, the neutrino is effectively destroyed, whereas in ν τ CC interactions, the outgoing tau-lepton may decay into lower-energy neutrinos [11]. In NC interactions, the incoming neutrino is not destroyed but cascades down in energy [12]. These flavordependent processes alter the neutrino flux as a function of the traversed path length [13], which allows for probing the neutrino cross section at high energies. The dip in theν e flux ratio due to the Glashow resonance [14] is visible in the right panel of Fig. 1 near E ν ¼ 6.3 PeV. The Glashow resonance occurs from the interaction of an electron antineutrino with a bound atomic electron and is independent of the CC and NC interactions of nucleons. A related, but subdominant, effect that has found renewed interest is the production of an on-shell W-boson off of the nucleus [15][16][17][18]. While this measurement is insensitive to the process, prospects for detection seem favorable with IceCube-Gen2 [19,20].
The IceCube Neutrino Observatory, an in-ice neutrino detector situated at the South Pole, is capable of detecting high-energy neutrinos originating from both Northern and Southern hemispheres [21][22][23][24]. IceCube comprises over 5000 digital optical modules (DOMs) encompassing approximately a cubic kilometer of ice [25][26][27]. The ice acts as a detection medium by which Cherenkov radiation from charged particles produced in neutrino interactions can be observed. The high-energy starting event (HESE) sample selects events that interact within a fiducial region of the detector across a 4π solid angle [24,28]. Here, we report a new cross section measurement using information from all three neutrino flavors with 7.5 years of data.
While SM calculations are generally consistent in the TeV-PeV energy range, few experimental measurements exist, and none have been performed with all three neutrino flavors [29,30]. Recently, an IceCube measurement of the neutrino DIS cross section using up-going, muon neutrinos gave a result consistent with the Standard Model [29]. The measurement in Ref. [30] used showers in publicly available HESE data with six years of data taking. This result, using the latest HESE sample with 7.5 years of data, includes classifiers for all three neutrino flavors and accounts for neutrinos from NC interactions and tau regeneration. Out of a total of 60 events above 60 TeV, 33 are also used in Ref. [30]. However, updates described in Ref. [28] affect the measurement and are incorporated in this work. These include updated parametrizations of absorption and scattering of light in the ice [31,32], more accurate atmospheric neutrino passing fractions [33], a likelihood construction that properly accounts for the stochastic nature of IceCube simulations [34], and improved systematics treatment [35].
As the sample updates are detailed in Ref. [28], this paper focuses on the results of the neutrino-nucleon cross section measurement. A brief description of the event selection is given in Sec. II. Section III details the analysis procedure. Section IV presents our Bayesian and frequentist results and compares them to existing measurements. We conclude in Sec. V.

II. EVENT SELECTION
The measurements presented here rely on a sample of high-energy events that start within a fiducial region of the IceCube detector [24,28]. In this context, events are taken to be the interaction byproducts of neutrino interactions, or background muons from cosmic-ray interactions in the atmosphere. The 90 m of the top and outer side layers of the detector, 10 m of the bottom of the detector, and a 60 m horizontal region near the highest concentration of dust in the ice are used as an active veto. Only events with fewer FIG. 1. Ratio of the arrival flux to surface flux for both electron neutrinos and antineutrinos as a function of E ν and zenith angle in IceCube detector coordinates. The flux at the surface is assumed to have a spectral index of γ ¼ 2. The core-mantle boundary is visible as a discontinuity near a zenith angle of 147°, and the enhanced suppression due to Glashow resonance is visible near 6.3 PeV in the electron antineutrino channel. Flux ratios for the other flavors are similar to that of electron neutrinos. than 3 photoelectrons (PEs) and fewer than 3 hit DOMs in the veto region within a predefined time window are kept. In addition, the total charge must exceed 6000 PE [28]. This removes almost all of the background due to atmospheric muons from the Southern sky. Neutrinos arriving from above and below the detector are included in the sample, thus allowing for constraints across the full allowed region in zenith.
Events are grouped into three experimental particle identification (PID) classifiers: cascades, tracks, and double cascades. These PIDs are related to the true interaction channel of the neutrino. Electromagnetic and hadronic showers appear cascadelike, stochastic energy losses from high-energy muons appear tracklike, and the production and subsequent decay of a tau can appear as double cascades (in addition to the other two PIDs) [28,36]. Since a NC interaction produces a hadronic shower, it is not directly distinguishable from a CC interaction. In addition, misclassifications can occur, and as such, the mapping from true to reconstructed observables is imperfect. To model such effects, detailed Monte Carlo (MC) simulations are performed, taking into account systematic variations in the ice model. The MC is then processed in an identical manner as the data. It thus provides the connection from the physics parameters of interest to the observed data events.

III. ANALYSIS METHOD
The cross section is probed by measuring neutrino absorption in the Earth, which effects an angular and energy dependence in the neutrino flux arriving at IceCube. In essence, we perform fits to the measured energy-zenith distribution for different hypotheses of the cross section, using a MC to model the expected signal under these hypotheses. A proper modeling accounts for neutrino interactions in the Earth as well as in the detector, with the event rate being proportional to the product of the cross section and the arrival flux. Knowledge of the number of events detected in data and the expected arrival flux can thus be used to measure the cross section. This is accomplished by binning the data and the expectation from MC (under a particular assumption of the cross section and flux) in observable space and then constructing a representative likelihood. The MC expectation varies as the cross section, as well as other systematic parameters, is modified. These modifications on theoretical parameters are captured via event-by-event reweighting of the MC [37]. For this analysis, the neutrino DIS cross section is parametrized as a function of four independent energy slices. Figure 2 illustrates the effect of scaling the DIS cross section up or down on the survival probability as a function of energy for a neutrino traveling through the full diameter of the Earth. It is plotted for each flavor individually as a function of the neutrino energy, E ν ; at a zenith angle of 180°; and for a surface flux with spectral index of γ ¼ 2.
The dependence on the spectrum arises from secondary neutrinos, which cascade down in energy and are produced in NC interactions and tau decay [11]. As the cross section increases, ΦðE ν Þ=Φ 0 ðE ν Þ decreases since the neutrinos are more likely to interact on their way through the Earth. The reason there is a slight flavor dependence is due to the fact that CC ν muons that lose most of their energy in the Earth before decaying due to their much longer lifetimes, can quickly decay to a lower-energy ν τ . Neutral-current interactions have a similar effect, and these effects are taken into account [13,29]. Furthermore, the dip in theν e flux ratio due to the Glashow resonance is again visible. This effect, taken over the full two-dimensional energy-zenith distribution, allows us to place constraints on the cross section itself.
In this paper, we report the neutrino DIS cross section as a function of energy under a single-power-law astrophysical flux assumption. Four scaling parameters, x ¼ ðx 0 ; x 1 ; x 2 ; x 3 Þ, are applied to the cross section given by CSMS [2] across four energy bins with edges fixed at 60 TeV, 100 TeV, 200 TeV, 500 TeV, and 10 PeV, where the indices correspond to the ordering of the energy bins from lowest to highest energies. Each parameter linearly scales the neutrino and antineutrino DIS cross section in each bin, while keeping the ratio of CC-to-NC contributions fixed. More specifically, for each x, neutrino events in MC are reweighted by x i ΦðE ν ; θ ν ; xÞ=ΦðE ν ; θ ν ; 1Þ, where Φ is the arrival flux as calculated by nuSQuIDS, E ν is the true neutrino energy, θ ν the true neutrino zenith angle, and x i the cross section scaling factor at E ν . The fixed CC-to-NC ratio implies that this analysis should not be interpreted as a direct test of the large extra dimensions model [5], which only applies to NC interactions. At energies above 1 TeV, the neutrino-nucleon cross sections for all three neutrino flavors converge. The cross section is therefore assumed to not depend on flavor in this measurement, but any differences in the arrival flux of ν e , ν μ , ν τ are taken into account. As the cross section is not flat in each bin, the effect of these four parameters is to convert it into a piecewise function where each piece is independently rescaled. Such an approach introduces discontinuities due to binning but allows for a measurement of the total neutrino-nucleon cross section as a function of energy. It also relaxes constraints based on the overall shape of the CSMS cross section and results in a more model-independent measurement. As the fit proceeds over all four bins simultaneously, bin-tobin correlations can be examined, though no regularization is applied.
The CSMS cross section is computed for free nucleon targets and does not correct for nuclear shadowing. The shadowing effect modifies nuclear parton densities and is stronger for heavier nuclei. At energies below 100 TeV, antishadowing can increase the cross section by 1%-2%, while above 100 TeV, shadowing can decrease the cross section by 3%-4% [38]. As this is a subdominant effect, we do not include it in this analysis. We do, however, consider the Glashow resonance in which an incidentν e creates an on-shell W − by scattering off an electron in the detector.
The effect on the expected arrival flux at the detector due to a modified cross section is calculated with nuSQuIDS, a neutrino propagation framework that properly takes into account destructive CC interactions, cascading NC interactions, and tau-regeneration effects [39]. A forward-folded fit is then performed, relying on MC to map each neutrino flavor to the experimental PID of tracks, cascades and double cascades [36], in the reconstructed zenith vs reconstructed energy distribution for tracks and cascades and in the reconstructed energy vs cascade length separation distribution for double cascades [28]. The fit uses the Poisson-like likelihood, L Eff , which accounts for statistical uncertainties in the MC and is constructed by comparing the binned MC to data [34]. The ternary PID of tracks, cascades, and double cascades is an additional constraint to the fit, which allows this measurement to incorporate interaction characteristics of all three neutrino flavors [36]. Using MC simulations, we can account for deviations between the true flavor and the PID and also estimate its accuracy. Under best-fit expectations, true ν e are classified as cascades approximately 57% of the time, true ν μ as tracks approximately 73% of the time, and true ν τ as double cascades approximately 65% of the time [28].
Systematic uncertainties are incorporated via the parameters listed in Table I. The parameter Φ conv (Φ prompt ) linearly scales the atmospheric neutrino flux normalization for neutrinos produced by π or K (charm meson) decay [40,41]. The parameters γ and Φ astro modify the spectral index and normalization of the astrophysical neutrino flux, respectively. The parameter Φ μ linearly scales the atmospheric muon flux normalization. The ratio of atmospheric neutrinos produced in pion vs kaon decay is governed by the π=K parameter. The parameter ν=ν governs the ratio of atmospheric neutrinos to antineutrinos. Finally, the parameter Δγ CR modifies the spectral index of the cosmic-ray flux [42]. Detector systematic studies were performed using Asimov data but had a negligible impact on the result. Priors on the nuisance parameters are given in Table I. The   TABLE I. Central values and uncertainties on the nuisance parameters included in the fit. Truncated Gaussians are set to zero outside the range. These modify the likelihood used in both the Bayesian and frequentist constructions. Their best-fit values over the likelihood space are also given. prior on γ astro is driven by the usual Fermi acceleration mechanism, allowing for a large uncertainty that covers those reported in a previous and independent IceCube measurement of the diffuse neutrino flux [43]. Such a large uncertainty minimizes the impact of changing the central value on the measured cross section. None of the x i parameters shifted by more than 1% in postunblinding checks where γ astro ¼ 3.0 AE 1.0. Out of all the nuisance parameters, γ astro and Φ astro exhibited the largest correlation with the cross section parameters. They are most strongly correlated with x 0 , the cross section in the lowest-energy bin. This is believed to be related to the fact that lower-energy neutrinos are subject to less Earth absorption so the main effect of varying the lowenergy cross section is a near-linear scaling at the detector. This makes x 0 essentially inversely proportional to the astrophysical flux. By allowing the cross section to float, the data seem to prefer the softer index, as given in Table I.
The interaction rate of high-energy neutrinos traveling through the Earth is also dependent on the Earth density. Here, we fix the density to the preliminary reference Earth model [44]. This is a parametric description of the density as a function of radial distance from the center of the Earth, evaluated using several sources of surface and body seismic wave data. Since the density uncertainty is at the few percent level, it is negligible in comparison to the flux uncertainty and is fixed for the purposes of this measurement [45].
Note that the Glashow resonance occurs for an incident ν e with an energy around 6.3 PeV and is not varied in the fit as it is calculable from first principles, using the known decay width of the W boson. However, unlike high-energy neutrino-nucleon scattering, the expected number of events due to the Glashow resonance is strongly dependent on the ratio of neutrinos and antineutrinos in the incident flux. We therefore performed a test that varied the astrophysical flux from a pure neutrino flux to a pure antineutrino flux. It was found only to have a minimal effect in the highest-energy bin, where the measurement uncertainty is largest. This is due in part to the steeply falling spectrum, which causes the flux at 6.3 PeV to be much smaller than that at lower energies. As the effect on the cross section is minimal, we keep the ratio of the flux of astrophysical neutrinos and antineutrinos fixed to unity.
We report both Bayesian highest posterior density (HPD) credible intervals and frequentist confidence intervals (CI). In the Bayesian construction, the posterior on the four scaling parameters is obtained with a Markov chain Monte Carlo (MCMC) sampler, emcee, marginalizing over nuisance parameters [46]. A uniform prior from 0 to 50 is assumed for all four cross section scaling parameters. Such a prior gives more weight to parameter values greater than 1. To test its effect, the MCMC was also run assuming a log-uniform prior, which gives results consistent with those assuming a uniform prior. The MCMC is sampled with 60 walkers over 5000 total steps, the first 1000 of which are treated as part of the initialization stage and discarded.
The frequentist confidence regions are obtained from a grid scan of the likelihood across four dimensions, profiling over the nuisance parameters and assuming Wilks theorem. For x 0 , x 1 , and x 2 , 15 equal-distant points are used from 0.1 to 5. For x 3 , 29 equal-distant points are used from 0.1 to 9.9. For each x on the mesh of these points, the likelihood is minimized over all other nuisance parameters. Confidence regions in two or one dimension are then evaluated by profiling across the other cross section parameters followed by application of Wilks theorem. Though the best-fit Φ prompt ¼ 0, the prompt component is expected to be a small contribution to the overall distribution. Thus, we expect Wilks theorem to hold asymptotically in the high statistics limit.
The zenith-dependent effect of the cross section on the event rate is shown in Fig. 3, assuming the best-fit, singlepower-law flux reported in Ref. [28], which is obtained using the CSMS cross section σ ¼ σ CSMS [2]. The degeneracy in the measurements of flux and cross section is broken by the different amounts of matter traversed by neutrinos arriving from different directions. To illustrate the effect of a modified cross section, two alternative expectations are shown for σ ¼ 0.2σ CSMS and σ ¼ 5σ CSMS under FIG. 3. The zenith distribution of data and the best-fit, singlepower-law flux expectation assuming σ CSMS (orange) [2]. Predictions from two alternative cross sections are shown as well, assuming the same flux. In the Southern sky, cos θ > 0, the Earth absorption is negligible, so the effect of rescaling the cross section is linear. In the Northern sky, cos θ < 0, the strength of Earth absorption is dependent on the cross section as well as the neutrino energy and zenith angle. the same best-fit flux assumption. In the Southern sky (cos θ > 0), the Earth absorption is negligible, and the event rate is simply proportional to the cross section. In the Northern sky (cos θ < 0), the strength of Earth absorption is dependent on the zenith angle, E ν , and the cross section. Absorption alters the shape of the event-rate zenith distribution in the Northern sky. For example, with σ ¼ 5σ CSMS and near cos θ ¼ −0.5, the attenuation of the arriving flux counteracts the increased neutrino interaction probability so that the event rate falls back to that expected from the CSMS cross section. Modifications of the neutrino cross section are thus constrained by the nonobservation of energy-dependent distortions in the zenith angle distribution.

IV. RESULTS
The CC cross section, averaged over ν andν, is shown in black in Figs. 4 and 5 for the Bayesian 68.3% HPD and frequentist one sigma intervals assuming Wilks theorem, respectively. As the scale factor is applied across the entire interval within an energy bin on the CSMS calculations, the shape is preserved within each bin. The central point in each energy bin corresponds to the expected, most probable energy in dN MC =d log E, the distribution of events in the MC along the x axis. This is chosen in lieu of the linear or logarithmic bin center to better represent where most of the statistical power lies in each bin. Since we assume a fixed CC-NC cross section ratio, the NC cross section is the same result relative to the CSMS prediction and so is not shown here.
In addition, the measurement based on HESE showers with 6 years of data is shown as orange crosses [30] in FIG. 4. The charged-current, high-energy neutrino cross section as a function of energy, averaged over ν andν. The Bayesian 68.3% HPD credible interval is shown along with two cross section calculations [2,4]. The credible intervals from a previous analysis [30] are also shown for comparison.
FIG. 5. The charged-current, high-energy neutrino cross section as a function of energy, averaged over ν andν. The Wilks 1-sigma CI is shown along with two cross section calculations [2,4]. The confidence intervals from Ref. [29] are also shown for comparison.
FIG. 6. The full posterior distribution of x as evaluated with emcee [46]. In the two-dimensional distributions, the 68.3% and 95.4% HPD regions are shown. In the one-dimensional distribution, the 68.3% HPD interval is indicated by the dashed lines. Fig. 4, and the previously published IceCube measurement, using up-going muon-neutrinos, is shown as the shaded gray region [29] in Fig. 5. Since credible intervals and confidence intervals have different interpretations, we do not plot them on the same figure. Note that both previous measurements extend below 60 TeV and are truncated in this comparison. Predictions from Refs. [2] and [4] are shown as the dashed and solid lines, respectively.
A corner plot of the posterior density, marginalized over all except two or one of the cross section parameters, is shown in Fig. 6. Similarly, two-dimensional profile likelihoods are shown in Fig. 7. Both exhibit little correlation between the various cross section parameters. The largest uncertainty arises for x 3 , which has the widest posterior distribution and flattest profile likelihood.
The Bayesian and frequentist results are consistent with each other, though again we caution that their intervals cannot be interpreted in the same manner. The results are compatible with the Standard Model and are summarized in Table II.

V. CONCLUSIONS
We have described a measurement of the neutrino DIS cross section using the IceCube detector. Variations in the neutrino cross section from Standard Model predictions modify the expected flux and event rate at our detector, and a sample of high-energy events starting within the fiducial volume of IceCube has been utilized to thus measure the neutrino cross section. Previous TeV-PeV scale neutrino cross sections have been measured by IceCube [29] using a sample of through-going muons and with cascades in the HESE sample [30]. This result, however, is the first measurement of the neutrino DIS cross section to combine information from all three neutrino flavors.
Our results are compatible with Standard Model predictions, though the data seem to prefer smaller values at the lowest-energy bin and higher values at the highestenergy bin. There do not seem to be strong correlations between the cross section bins, though large uncertainties due to a dearth of data statistics make it difficult to draw strong conclusions. With additional data, or with a combined fit across multiple samples, more precise measurements are foreseen in the near future [47].