An eV-scale sterile neutrino search using eight years of atmospheric muon neutrino data from the IceCube Neutrino Observatory

The results of a 3+1 sterile neutrino search using eight years of data from the IceCube Neutrino Observatory are presented. A total of 305,735 muon neutrino events are analyzed in reconstructed energy-zenith space to test for signatures of a matter-enhanced oscillation that would occur given a sterile neutrino state with a mass-squared differences between 0.01\,eV$^2$ and 100\,eV$^2$. The best-fit point is found to be at $\sin^2(2\theta_{24})=0.10$ and $\Delta m_{41}^2 = 4.5{\rm eV}^2$, which is consistent with the no sterile neutrino hypothesis with a p-value of 8.0\%.

Wulff, 11 X. W. Xu, 7 Y. Xu, 48 J. P. Yanez, 24 G. Yodh, 29 S. Yoshida, 15 T. Yuan, 36 Z. Zhang, 48 and M. Zöcklein (IceCube Collaboration) INTRODUCTION The three-flavor massive neutrino oscillation formalism has been well-established experimentally [1][2][3][4]. The standard paradigm has also been challenged, by several experiments exhibiting anomalous ν e (ν e ) appearance in ν µ (ν µ ) beams [5,6]. These anomalies can be interpreted as evidence for subleading oscillations of ν µ → ν e orν µ →ν e caused by additional neutrinos with large mass-squared differences in the range of ∆m 2 ∼ 0.1 − 10 eV 2 [7]. On the other hand, measurements of the Z-boson decay width to invisible final states demonstrate that only three light neutrinos participate in weak interactions [8], so any additional neutrino flavor states must be nonweakly-interacting, or "sterile." The simplest such model is referred to as a "3+1" model, where in addition to the three known mass states, a fourth heavier one is added.
The relationship between the flavor and mass states is described by a unitary matrix, U P N M S , which in the three-neutrino model can be parameterized in terms of three mixing angles and one oscillation-accessible CPviolating phase. Adding a sterile state expands the mixing matrix to four dimensions, in which the added degrees of freedom can be parameterized by introducing three new rotations with angles θ 14 , θ 24 , and θ 34 , and two new oscillation-accessible CP -violating phases, δ 14 and δ 24 . The oscillation phenomenology of the 3+1 model adds both shorter baseline vacuumlike oscillations, and also novel oscillation effects in the presence of matter [9][10][11][12][13]. For eV-scale sterile neutrino states, for example, a matter-enhanced resonance [14][15][16][17][18][19] would result in the near complete disappearance of TeV-scale muon antineutrinos passing through the Earth's core, as shown in Fig. 1. By measuring and characterizing the flux of atmospheric neutrinos in the GeV to PeV energy range, the IceCube Neutrino Observatory is uniquely positioned to search for such matter-enhanced oscillations, a smokinggun signature of eV-scale sterile neutrinos. Testing the 3+1 model as an explanation of shortbaseline anomalies and constraining its free parameters requires measurements in multiple oscillation channels, including ν µ → ν µ [21][22][23][24][25][26][27][28][29], ν e → ν e [30][31][32][33][34][35][36][37][38], and ν µ → ν e [5,6,[39][40][41]. Fits to global data [20,42,43] find a strong preference for models with sterile neutrinos over the standard three-neutrino paradigm. However, even at the most preferred values of ∆m 2 ∼ 1 eV 2 , the mixing angles required to viably explain anomalies in the ν µ → ν e andν µ →ν e channels are in strong tension with measurements of ν µ andν µ disappearance. Evidence for oscillation effects beyond the three-neutrino paradigm inν µ disappearance are yet to be observed [43]. One of these nonobservations was made by IceCube, using a sample of 20,145 atmospheric ν µ andν µ events collected over one year of detector livetime [26].
This paper updates IceCube's high-energy sterile neutrino search using an eight-year dataset and improved event selection. The sample includes 305,735 wellreconstructed charged-current ν µ andν µ events collected from May 13 th 2011 to May 19 th 2019. Events are binned uniformly in log(E µ reco ) spanning E µ reco ∈ Atmosphericνµ disappearance probability vs. true energy and cosine zenith at the globally preferred sterile neutrino hypothesis of Ref. [20] (∆m 2 41 = 1.3 eV 2 , sin 2 (2θ24) = 0.07, sin 2 (2θ34) = 0.0). Effects include a matter-enhanced resonance at TeV energies, neutrino absorption at high energy and small zenith, and vacuumlike oscillation at low energies. The matter-enhanced resonance appears only in the antineutrino flux for the case of small angles and ∆m 2 41 > 0. Vertical white lines indicate transitions between inner to outer core (cos(θ true ν ) = −0.98) and outer core to mantle (cos(θ true z ) = −0.83).  [500 GeV, 9976 GeV] in 13 bins and uniformly in cos θ reco z spanning the up-going region from −1.0 to 0.0 in 20 bins. The event counts in each bin are used as inputs to a likelihood-based analysis to test for evidence of eV-scale sterile neutrinos.
The increased sample size of this analysis with respect to Ref [26] has been accompanied by a commensurate improvement in the precision of treatments of systematic uncertainties and statistical methods. This paper summarizes these advances and presents the main results of this search. A companion paper, Ref. [44], contains a more detailed exposition of the technical aspects of the analysis, as well as alternate interpretations of the data in a wider space of sterile neutrino parameters.

ICECUBE UP-GOING TRACK SAMPLE
The IceCube Neutrino Observatory is a cubickilometer neutrino detector buried in the Antarctic glacier [45]. It is comprised of photomultiplier tubes enclosed in glass pressure housings called "Digital Optical Modules" (DOMs) [46]. These are arranged in vertical strings on a hexagonal lattice. The main array consists of 78 strings spaced 125 m apart, each supporting 60 downward-facing DOMs with a 17 m vertical spacing. A denser array called DeepCore [47] instruments the clearest part of the ice within the main array. The eight strings of DeepCore are arranged with lateral spacing between 42 m and 72 m and vertical DOM separation of 7 m. This analysis uses the complete set of IceCube DOMs in both the main array and DeepCore.
The majority of IceCube events are produced by highenergy muons and neutrinos from cosmic-ray air showers. Down-going (cos θ true z > 0) atmospheric muons (and antimuons) can penetrate the 1450 m vertical overburden of the detector, triggering at a rate of ∼3 kHz [48]. These events dominate the Southern-hemisphere through-going sample. Up-going atmospheric muons, on the other hand, are effectively removed by the large overburden provided by the Earth. Thus, muons originating from the Northern hemisphere are dominated by those produced in charged-current neutrino interactions. A charged-current ν µ interaction will produce a forward secondary muon, with an energy typically between 50% and 80% of that of the parent ν µ [49]. The muon travels through the ice emitting Cherenkov radiation. While photons travel tens to hundreds of meters before being absorbed by the impurities in the ice [50][51][52], muons with TeV energies are able to penetrate multiple kilometers of ice before falling below the Cherenkov threshold [53,54]. This produces an extended tracklike signature. These events originate either inside of the detector or from a target volume extending meters to kilometers outside the array, depending on energy [53,55].
Events used in this analysis first pass a filter that selects muon-like events for satellite transmission to the North, and are then subject to further data-reduction techniques to reject low-energy and poorly reconstructed tracks. Only data periods with 86 active IceCube strings and greater than 5,000 active DOMs in the detector are considered. A high-level event selection is applied, leveraging morphology, measures of track reconstruction quality, and the expected transmission of signal events through the zenith-dependent overburden, explained in detail in Ref. [44] and based on Ref. [56]. The reconstructed energy and direction of each event is calculated according to the time and geometry of light detected throughout the array [57,58]. The angular resolution σ cos θz varies between 0.005 and 0.015 and energy resolution of σ log 10 Eµ ∼ 0.5, as in the previous version of this analysis [26]. The energy distribution of selected events is shown in Fig. 2.
Cosmic-ray muon background contamination is assessed using CORSIKA [59], with primary cosmic-ray energies ranging from 600 GeV to 10 11 GeV. Approximately 10% of the dataset of neutrino events are predicted to contain a coincident cosmic-ray muon within the readout frame. The ν µ and cosmic-ray muon tracks are separated into sub-events using an event splitter, and each sub-event is treated independently in the event selection. After splitting and event selection, the sample is predicted to be > 99.9% pure in ν µ /ν µ induced events [44].

STERILE NEUTRINO ANALYSIS
In this analysis, we consider a sterile neutrino model parameterized by one mass-squared difference, ∆m 2 41 , and one mixing angle, sin 2 (θ 24 ). For each hypothesis point on a grid of ∆m 2 41 from 10 −2 eV 2 to 10 2 eV 2 and sin 2 (2θ 24 ) from 10 −3 to 1, the neutrino flux incident on the detector is calculated using the four-flavor formalism.
The neutrino flux includes contributions from both atmospheric and astrophysical neutrinos. The conventional atmospheric ν µ andν µ flux is produced by the decay of pions and kaons and is calculated using the MCEq cascade equation solver [60,61]. The hadronic interactions are modeled with Sibyll2.3c [62]. The primary cosmicray spectrum is a three-population model [63,64], in which each population contains five groups of nuclei. The zenith-dependent seasonal atmospheric density profile, through which the cascade develops, is determined using data from the Atmospheric Infrared Sounder (AIRS) satellite [65]. The prompt ν µ component from the decay of charmed mesons is implemented as in Ref. [66]. The astrophysical neutrino flux is assumed to have equal parts of each neutrino flavor and to be symmetric in neutrinos and antineutrinos [67][68][69]; be isotropically distributed; and have a single power-law energy spectrum consistent with previous IceCube measurements [70]. These fluxes are subject to a suite of systematic uncertainties, summarized in the following section.
For each sterile neutrino hypothesis, the atmospheric and astrophysical neutrino fluxes are propagated through the Earth using the nuSQuIDS neutrino evolution code [71,72]. This accounts for both coherent and non-coherent interactions [73]; namely charged-current, neutral-current, and Glashow resonance interactions [74], as well as tau-neutrino regeneration [75]. We use the CSMS [76] neutrino-nucleon cross section to describe both interactions during neutrino propagation and near the detector. This requires as an input the Earth density profile, which we parameterize via the spherically symmetric PREM model [77]. Using the above, we obtain a prediction for the up-going ν µ flux at the detector for each physics parameter point. These fluxes are used to weight detector Monte Carlo (MC) event sets, with effective livetime ≥ 50× the sample size.
We account for systematic uncertainties by means of nuisance parameters, which reweight the MC by applying continuous parameterizations of the effects discussed in the following section. We then compare the data to expectation using a modified version of the Poisson likelihood to account for MC statistical uncertainty [78]. For our frequentist analysis, the likelihood is profiled over the eighteen nuisance parameters to construct a test statistic. Frequentist contours are constructed using Wilks' theorem [79], validated at an array of parameter points using MC ensembles and the Feldman-Cousins [80] procedure.
A Bayesian hypothesis test is also performed, by means of comparing the model evidences [81] with respect to the no sterile neutrino hypothesis. The model evidences, as a function of sterile neutrino parameters, are computed by integrating the likelihood over the nuisance parameters using MultiNest [82]. These two statistical approaches are complementary: the Bayesian approach conveys the likelihood of the model given observed data and prior knowledge, whereas the frequentist approach yields intervals that are likely to contain the true model parameters for repeated experiments, enabling direct comparison with previous publications.

SYSTEMATIC UNCERTAINTIES
Dominant sources of uncertainty derive from the shape and normalization of astrophysical and atmospheric neutrino fluxes; the bulk properties of the South Pole ice; the local response of the IceCube DOMs; and neutrino interaction cross sections. Other uncertainties, such as the Earth density profile, neutrino interactions in the rock/ice transition region, prompt neutrino flux, and ν µ /ν µ astrophysical ratio were investigated but established as negligible relative to statistical uncertainty.
Atmospheric Neutrino Flux: In the relevant energy range the spectrum of cosmic-ray primaries follows approximately an E −2.65 energy (E) dependence. To account for the uncertainty in the cosmic-ray spectral index, we apply a spectral shift ∆γ with an uncertainty of 0.03 pivoting at 2.2 TeV [83][84][85][86]. The meson production uncertainty in the interaction between the primary cosmic ray and air and in subsequent hadronic interactions is described through the Barr et al. scheme [87]. In this scheme, the uncertainty in the differential cross section for meson production is quantified in regions of primary proton energy E p and meson fractional momenta x lab . The charged-kaon production yield carries the leading uncertainty. We parameterize its production over three kinematic regions: x lab < 0.1 and E p > 30 GeV; x lab ≥ 0.1 and 30 GeV < E p < 500 GeV; and x lab > 0.1 and E p ≥ 500 GeV. We include two collider-constrained nuisance parameters for each region, one for particles and one for antiparticles, which rescale the production cross section. The highest-energy uncertainties are obtained through extrapolation, and both the scale and energy dependence have ascribed uncertainties. Kaon energy losses by interaction with oxygen and nitrogen nuclei are accounted for via the total kaon-nucleus cross-sectional uncertainty [88]. The charged-pion production and interaction uncertainties were studied and found negligible. The atmospheric density profile is inferred from the zenith-dependent vertical temperature profile measured by the AIRS satellite. To incorporate its uncertainty, showers are recomputed through randomly perturbed density models within the statistical and system-atic uncertainties reported in the AIRS measurements. Finally, the total conventional atmospheric ν µ flux carries an additional 40% normalization uncertainty following Ref. [61].
Astrophysical Neutrino Flux: The central astrophysical model is a single power law with an equal normalization for all neutrino and antineutrino flavors at 100 TeV of 0.787×10 −18 GeV −1 sr −1 s −1 cm −2 and a spectral index of 2.5. The Gaussian priors on the normalization and spectral index are correlated and selected to accommodate all IceCube astrophysical neutrino flux measurements to date [70,[89][90][91][92][93], with allowed spectral indices of γ astro ∼ 2.2−2.8 at 68% confidence level (C.L.). This represents a significant contribution to the total flux in the top two energy bins, depending strongly on the value of γ astro .
Bulk Ice Model: The uncertainty associated with the measured scattering and absorption of the undisturbed glacial ice is implemented as described in Ref. [94]. This treatment expresses the depth dependence of the ice optical properties using a Fourier decomposition. The covariance of the Fourier mode coefficients are determined using LED flasher calibration data [52]. Only the six lowest modes contribute a sizeable shape difference in the reconstructed event distributions. The effect of these modes is parameterized using two empirical energy-dependent basis functions. The two associated amplitudes are incorporated as nuisance parameters with a correlated bivariate Gaussian prior.
DOM Response and Local Ice Effects: The ice in the immediate vicinity of the DOMs has optical properties distinct from the bulk ice between strings [95], caused by bubble formation during the refreezing process after their deployment. This introduces uncertainties via two effects. First, the global photon detection efficiency is impacted. This is modeled by an efficiency correction with an effectively flat prior, ultimately constrained with a tight posterior through its effect on the overall energy scale. Second, the bubble column influences the angular dependence of photon detection. This is encoded in two parameters tuned to detailed optical simulations of bubble scattering near the DOM [96], with only one having a substantial impact.
Neutrino Cross Section: The neutrino-nucleon cross section enters the analysis in two ways, influencing: 1) absorption during the neutrino propagation through the Earth [49,97] and 2) the rates and inelasticities of interactions near the detector [49,76,98]. The latter source of uncertainty was previously investigated in Refs. [99,100] and found to be negligible. The former is found to be non-negligible and is taken into account by separately parameterizing the change in neutrino absorption when the ν µ andν µ cross sections are scaled. The priors on these parameters are fixed at the largest uncertainties in our energy range from Ref. [76], which are 3% for ν µ and 7% forν µ . The signal shape at the best-fit point compared to the null hypothesis with the same nuisance parameters. Bottom: data compared to the null hypothesis with the nuisance parameters held at the same values.
The frequentist analysis best-fit point is ∆m 2 41 = 4.5 eV 2 and sin 2 (2θ 24 ) = 0.10. At this point, the largest nuisance parameter pull was observed in the cosmic-ray spectral index, which shifted the cosmic-ray spectrum by 0.068 (2.3σ); the other nuisance parameter best-fit values are within one sigma of their respective central values and can be found in the accompanying Ref. [44]. Fig. 3 shows the signal shape at the best-fit point, given the best-fit nuisance parameters, as well as the pull between data and no sterile neutrino hypothesis, evaluated at those same nuisance parameters. Fig. 4 shows the 90% and 99% C.L. contours calculated according to Wilks' theorem with two degrees of freedom. Sensitivity envelopes, illustrating symmetrically counted ensembles of 68% and 95% non-closed contours derived from 2,000 pseudoexperiments, are shown overlaid for the 99% contour. The IceCube 90% C.L. preferred region is consistent with constraints from previous ν µ disappearance  [22][23][24][25][26][27]101]; where results were not available at 99% C.L., methods of Ref. [20] were applied using public data releases. Finally, the star marks the analysis best-fit point location.
experiments, and the 99% contour is stronger than other exclusion limits at values of ∆m 2 up to 1 eV 2 . Fig. 5 shows the corresponding Bayesian result, where the point-wise Bayes factor is calculated relative to the no sterile neutrino hypothesis. The best-model location is at ∆m 2 ∼ 4.5 eV 2 and sin 2 (2θ 24 ) ∼ 0.9 and is strongly preferred, by a factor of 10.7, to the no sterile neutrino hypothesis. Contours are drawn in logarithmic Bayes Factor steps of 0.5, quantifying strength of evidence [102].
The best-fit point and inferred confidence regions are found to be robust under the removal of any one of the eight years of data, showing only minor changes in the contour position. This is also the case for removal of any of the following group of uncertainties: neutrino cross sections, detector effects, atmospheric flux, and astrophysical flux. Details can be found in Ref. [44]. Furthermore, a similar best-fit point is obtained when fitting any one year of data independently, suggesting a small effect of physical or systematic rather than statistical origin.
The difference in likelihood to the null hypothesis is 4.94, corresponding to a p-value of 8% against the null hypothesis. The location of this point was found to be compatible with expectations based on simulated no sterile neutrino pseudoexperiments, which by definition pro- duce closed contours at 90% C.L. in 10% of trials.
In summary, we have studied 305,735 up-going atmospheric and astrophysical muon-neutrinos to search for evidence of eV-sterile neutrino signatures. The best-fit point is consistent with the no sterile neutrino hypothesis at a p-value of 8%. Because of its unique statistical strength this result is expected to have a substantial impact on the global sterile neutrino landscape.