Combining sterile neutrino fits to short-baseline data with IceCube data

Recent global fits to short-baseline neutrino oscillation data have been performed finding preference for a sterile neutrino solution ( 3 þ 1 ) over null. In the most recent iteration, it was pointed out that an unstable sterile neutrino ( 3 þ 1 þ decay) may be a better description of the data. This is due to the fact that this model significantly reduces the tension between appearance and disappearance datasets. In this work, we add a 1-year IceCube dataset to the global fit obtaining new results for the standard 3 þ 1 and 3 þ 1 þ decay sterile neutrino scenarios. We find that the 3 þ 1 þ decay model provides a better fit than the 3 þ 1 , even in the presence of IceCube, with reduced appearance to disappearance tension. The 3 þ 1 þ decay model is a 5 . 4 σ improvement over the null hypothesis and a 2 . 8 σ improvement over the standard 3 þ 1 model. DOI: 10.1103/PhysRevD.101.055020


I. INTRODUCTION
Over the last decades, anomalies observed in shortbaseline neutrino oscillation experiments [1,2] have motivated global fits that expand from a model containing three neutrinos to one that also includes a fourth mass state and a nonweakly interacting flavor state [3][4][5][6]. Hence, these socalled "3 þ 1 models" extend the neutral-lepton mixing matrix to be 4 × 4, introducing three new independent mixing elements: U e4 , U μ4 , and U τ4 . The preferred regions of jU e4 j, jU μ4 j, and the new squared-mass difference, Δm 2 41 , have been reported in a recent global fit of the short-baseline data to a 3 þ 1 model [5].
A 3 þ 1 model can simultaneously explain anomalies observed in "appearance experiments" ( ν  [7][8][9]. However, in a 3 þ 1 model, the combined anomalies also predict a signal in ν ð−Þ μ disappearance experiments. For oscillations in vacuum, observed through charged-current scattering, the elements jU e4 j 2 and jU μ4 j 2 define three mixing angles that characterize the amplitude of ν e disappearance, ν μ disappearance, and ν μ → ν e appearance, respectively, sin 2 2θ ee ¼ 4ð1 − jU e4 j 2 ÞjU e4 j 2 ; sin 2 2θ μμ ¼ 4ð1 − jU μ4 j 2 ÞjU μ4 j 2 ; Thus, the mixing angles measured in the three types of searches are not independent. In fact, for small mixing elements, they are approximately related in the following way: sin 2 2θ ee sin 2 2θ μμ ≈ 4 sin 2 2θ μe : The squared-mass splitting must also be consistent for all three types of oscillations. The probability that all shortbaseline data, observed anomalies, and constraints from null observations is a realization of the 3 þ 1 model is found to be extremely small [4]. This statement follows from a recent parameter goodness of fit test [10] ("PG test") performed on the global data, which quantifies the tension as a disagreement at the 4.5σ [4] level.
If the anomalies are due to new physics, the underlying model may be more complicated than 3 þ 1. We have shown that adding another additional mass state, a model known as "3 þ 2," does not reduce the tension [5]. Other modifications are also possible [11][12][13][14][15][16][17][18][19][20][21][22], and we have developed a 3 þ 1 model that incorporates decay of the largest mass state, m 4 [5,23]. In fact, this category of 3 þ 1 models was first considered as an explanation of the LSND observation in [24] and was later realized to weaken the muon-neutrino disappearance constraints in [23]. More recently, it has been suggested that if one considers visible decay with one or more decay daughters being detectable, the low-energy excess of MiniBooNE can be well described by this category of models [25][26][27]. In the case of invisible decay, where the heavy neutrino has two decay daughters that are undetectable, this introduces only one additional parameter beyond the 3 þ 1 model, the lifetime of the heavy neutrino, τ. However, it also introduces two new, negligible-mass particles, the decay daughters, into the phenomenology. This model reduces the tension, measured by the parameter goodness of fit test, significantly to 3.2σ [5]. While this still suggests poor agreement, it points to the possibility of additional effects in the shortbaseline sample.
The global fits described above were limited to experiments studying vacuum oscillations at L=E ∼ 1 to 10 eV 2 . However, the IceCube experiment offers a relevant dataset that is very different from these short-baseline experiments. This experiment searches for a resonance signature from sterile-induced matter effects in upward-going antineutrinos, which have traveled through the Earth [28,29]. The oscillation amplitude is no longer given by the vacuumoscillation mixing angle relations, Eq. (1), and instead takes on a complicated form, where sin 2 2θ μμ , the amplitude of the disappearance, depends on jU μ4 j and also jU τ4 j; see Refs. [30][31][32][33] for a detailed discussion. Thus, IceCube brings additional information to the global fits, since the vacuum-oscillation results are insensitive to jU τ4 j [3].
A 3 þ 1 model has been explored for a combined shortbaseline and IceCube global fit for the first time in [3] and later in [4]. In this paper, we will follow the same procedure of [3] to include IceCube in the latest global fits given in [5], and we will expand the result to include 3 þ 1 þ decay.

II. ICECUBE
The IceCube Neutrino Observatory is located in the Antarctic continent close to the geographic South Pole [34]. IceCube is a gigaton-scale ice-Cherenkov detector made out of arrays of photomultiplier tubes encapsulated in pressure-resisting vessels buried in the clear Antarctic glacier ice [35]. IceCube has measured the atmospheric neutrino spectrum from 10 GeV [36], using a denser inner part of the detector called DeepCore [37], to 100 TeV [38]. The atmospheric neutrino flux is dominated by the socalled conventional component, which is due to the decay of kaons, pions, and muons [39]. At the lower energies, below 100 GeV, neutrinos from pion and muon decay dominate the neutrino flux and IceCube has used them to measure the atmospheric oscillation parameters [36]. In this energy range, IceCube is also sensitive to eV-scale sterile neutrinos by looking for modifications on top of the predicted standard three-neutrino oscillation pattern [40], as similarly done by SuperKamiokande [41]. At TeV energies, neutrino oscillations driven by the known squared-mass differences turn off, as their oscillation length becomes much larger than the Earth's diameter. It was pointed out in [28], that at TeV energies matter effects will induce a large disappearance of muon-antineutrinos for squared-mass differences compatible with the LSND anomaly, due to effects previously studied in broader contexts in [42][43][44][45][46]. IceCube has performed a search for sterile neutrinos using 1 year of data [29]. This analysis used a high-purity muon-neutrino event selection designed to search for an astrophysical component in the northern sky [47]. The analysis was performed with neutrinos between 400 GeV and 20 TeV in energy; this range was chosen to avoid contamination from high-energy astrophysical neutrinos and to minimize uncertainties in local ice properties at lower energies.
The IceCube Collaboration provided a data release associated with this analysis [29]. This data release includes over 20,000 atmospheric ν μ andν μ events, as well as the Monte Carlo that was used. In this work, we use the analysis tools we developed in [3,23], which use the nuSQuIDS package to calculate neutrino oscillation probabilities [48,49]. In our analysis of the IceCube data, we consider two sources of uncertainties. One source of uncertainty is the atmospheric neutrino flux, which we parametrize by means of nuisance parameters. We consider the overall normalization of the atmospheric flux, the cosmic-ray slope, the uncertainty in the atmospheric density, the ratio of kaon-topion production yields, and the ratio of neutrino-toantineutrino production. The second source of uncertainty is the detector. As discussed in [50,51], in this analysis we restrict ourselves to the leading detector systematic, namely the overall efficiency of the IceCube modules. Following [29], the IceCube log-likelihood can be written as where x i is the number of data events in the ith bin; μ i is the Monte Carlo expectation for the number of events in the ith bin, assuming sterile neutrino parameters ⃗ θ and nuisance parameters ⃗ η; and each nuisance parameter, η, has a Gaussian constraints of mean,η, and standard deviation, σ η . We use the systematic treatment and analysis framework from [23] to calculate the IceCube likelihood.
For both the 3 þ 1 and 3 þ 1 þ decay global fits, the IceCube likelihood was calculated for a randomly selected subset of parameter-set points from the corresponding Markov Chain Monte Carlo from [5]. This global analysis used experiments where vacuum-oscillation probabilities are valid and, in the stable 3 þ 1 model, those were parametrized in terms of jU e4 j, jU μ4 j, and Δm 2 41 . The datasets used include all relevant muon-(anti)neutrino disappearance, electron-(anti)neutrino disappearance, and muon-to-electron (anti)neutrino appearance measurements; a list of the experiments used can be found in the Supplemental Material [52]. There are four notable exceptions: (i) the Daya Bay result [53], as it does not significantly impact the preferred best-fit region (Δm 2 41 > 1 eV 2 ) [3]; (ii) the reactor antineutrino anomaly, due to the large uncertainties in the reactor flux modeling made manifest by the so-called 5 MeV reactor bump [54,55]; (iii) the low-energy atmospheric neutrino results from SuperKamiokande [41] and IceCube/DeepCore, [40] as these constraints only reach jU μ4 j > 0.2 which are covered by other experiments in the squared-mass splitting region of interest; and (iv) the recent MINOS+ two-detector fit [56,57], due to the lack of clear explanation of the origin of the sensitivity in the high-mass region and increased systematic model dependence, so instead we use the traditional far-to-near ratio measurement by MINOS [58].
Each point is a unique combination of jU e4 j, jU μ4 j, and Δm 2 41 ; the lifetime, τ, is an additional parameter for the 3 þ 1 þ decay scenario. In this work, we set jU τ4 j to zero, as the short-baseline experiments are insensitive to it and this is a conservative choice in the case of IceCube [3]. Downsampling points from the global analysis [5] is necessary because the IceCube likelihood calculation is computationally time intensive. We include the best-fit parameter-set points corresponding to 3 þ 1 and 3 þ 1 þ decay found in [5]. For the 3 þ 1 analysis, 49,000 points were used, while for the 3 þ 1 þ decay analysis 82,000 points were used.
To convert the IceCube likelihood into a χ 2 , we calculate the log-likelihood ratio [3], where SP fðx i Þg is the saturated Poisson likelihood, and then assume Wilks' theorem, namely χ 2 ¼ −2 ln LR [59,60]. To incorporate IceCube into the global fits, we add the IceCube χ 2 to the χ 2 from the short-baseline-only fits [5].
To determine the effect that IceCube data has on the tension for both the 3 þ 1 and 3 þ 1 þ decay models, we calculate the IceCube likelihood for a random downsampling of parameter points from the recent fits to only disappearance short-baseline experiments. We convert the IceCube likelihood to a χ 2 and add it to the short-baseline disappearance χ 2 . The fit to the appearance experiments remains unchanged.

III. RESULTS
In this section, we summarize the result of incorporating IceCube into our recent short-baseline-only global fit. Results for fitting a standard 3 þ 1 model are discussed in Sec. III A, while results for fitting a 3 þ 1 þ decay model are shown in Sec. III B. A summary of the χ 2 values for various fits and additional figures for all analyses performed in this work are given in the Supplemental Material [52].

A. 3 + 1
In this section, we report the impact of adding IceCube to the short-baseline-only 3 þ 1 fit. Figure 1 shows the frequentist allowed regions without IceCube data (top) and with IceCube data incorporated (bottom), for a random downsampling of parameter points used in the recent global fit [5]. The allowed regions are shown in terms of Δm 2 41 and sin 2 2θ μe , given by Eq. (1), at 90%, 95%, and 99% confidence levels. The best-fit points are indicated in each panel with a star. Before including IceCube data, the downsampled points produce allowed regions that are consistent with our previous result and are shown here for completeness. Incorporating IceCube data pushes the island at Δm 2 41 ≈ 1.3 eV 2 to slightly smaller sin 2 2θ μe and creates new islands at higher masses at the 99% confidence level. A table of the best-fit parameters is given in Table I and the best-fit Δχ 2 is given in Table III. The obtained bestfit values on the 3 þ 1 model are in tension with the MINOS+ results not included in this work for the reasons mentioned in Sec. II. The significance of the best-fit point of the 3 þ 1 model compared to the null hypothesis, i.e., three neutrinos, is 4.9σ, slightly lower than it is without IceCube, 5.1σ. We have performed the parameter goodness of fit test to quantify the tension in the dataset. The significance of this tension with IceCube included is 4.8σ, which is slightly higher than the value without including IceCube, which is 4.5σ.

B. 3 + 1 + decay
Global fit results for the 3 þ 1 þ decay model are shown in Fig. 2. The top row of this figure shows the allowed regions without including IceCube, which are consistent with previous results, while the bottom row shows the allowed regions with IceCube incorporated. The allowed regions without considering IceCube span 3 orders of magnitude in τ, and are shown in three panels, each corresponding to the indicated range in τ. The allowed regions without IceCube occur in two distinct regions in Δm 2 41 . IceCube data eliminate the lower squared-mass difference island and collapses the allowed values to a narrow range around Δm 2 41 ≈ 1.4 eV 2 . A table of the best-fit parameters is given in Table II and the best-fit Δχ 2 is given in Table III. The significance of the best-fit point of the 3 þ 1 þ decay model compared to the three-neutrino hypothesis is 5.4σ, slightly lower than it is without IceCube, which is 5.6σ. Notably, the 3 þ 1 þ decay model is preferred to the standard 3 þ 1 model by 2.8σ. Our global-fit combines three different oscillation channels and the 3 þ 1 þ decay model introduces features in all three; see Supplemental Material [52] for details. In the reactor electron-antineutrino disappearance measurements, the best-fit decay reduces the oscillation features at low energies where the statistics are larger. Similarly, in muon-neutrino disappearance measurements, it reduces the oscillation amplitude for baselines greater than the decay length. In the appearance channel, for scales relevant for MiniBooNE, the appearance probability maintains the spectral shape but has increased normalization. Tension in the fits remains, at 3.5σ with IceCube, yet it is reduced by about 1.3σ compared to the standard 3 þ 1 model. Finally, note that the inclusion of absolute reactor flux normalization could yield important information to distinguish between models. The preferred  parameters of the 3 þ 1 þ decay increase the electronantineutrino disappearance from approximately 3% in the 3 þ 1 best-fit point to 10% for reactor neutrino scales of interest. However, the estimated neutrino flux uncertainties range from 2% to 6% in 1-6 MeV [61]; note also that the baseline flux model predictions did not foresee the appearance of a deviation at the 10% level known as the 5 MeV bump. Given these uncertainties, we have performed a normalization independent analysis. Including a reliable reactor flux uncertainty estimation would impact our obtained best-fit point.

IV. CONCLUSION
We have studied the impact of adding 1 year of IceCube high-energy atmospheric data to our short-baseline light sterile neutrino global fits. We considered two models in this work: one where the heavy neutrino mass state is stable and one where it decays to invisible particles. We summarize our findings here: (i) For the case of stable neutrino states, we find that the best-fit solution to the short-baseline data is only slightly changed; the largest change is in jU μ4 j which goes from 0.135 to 0.131. Adding IceCube data makes two new solutions appear at the 99% C.L. at larger squared-mass differences; these are at approximately 10 eV 2 and 5 eV 2 , with sin 2 2θ μe ∼ 10 −3 . Even though the IceCube analysis is a null-like result, the best-fit 3 þ 1 point is still significantly preferred over the null hypothesis by a Δχ 2 ¼ 32 with 3 degrees of freedom. This corresponds to a 4.9σ rejection of the null model. As expected, adding the IceCube dataset increases the tension between appearance and disappearance datasets which, measured using the parameter goodness of fit test, worsens from a p-value of 3.4 × 10 −6 without IceCube to 8.1 × 10 −7 with it. (ii) For the model where the heaviest mass state is allowed to decay into invisible particles, the preferred values of the parameters have changed significantly. The new best-fit point results in a squared-mass difference of 1.35 eV 2 , a lifetime of 4.5 eV −1 , and mixing elements jU e4 j ¼ 0.238 and jU μ4 j ¼ 0.105. As in the case of a stable neutrino mass states, this scenario is preferred over the null three-neutrino hypothesis at high significance, 5.6σ without IceCube and 5.4σ with it. This model was already preferred with respect to the 3 þ 1 scenario at the 2.7σ level [5], but with the addition of the IceCube data this preference increases to the 2.8σ level. Finally, the tension is slightly worse when adding the IceCube data in this model; the p-value for the parameter goodness of fit test decreases from 8.0 × 10 −4 to 2.7 × 10 −4 . Nevertheless, this model is still a factor of ∼300 improvement in p-value with respect to the stable sterile neutrino scenario. In conclusion, the 3 þ 1 þ decay scenario best-fit parameters have been significantly changed with the addition of 1 year of IceCube data. The new best-fit solution is still an improvement over the 3 þ 1 scenario since it improves the overall fit and reduces the tension between appearance and disappearance experiments. We expect that the upcoming 8-year IceCube high-energy sterile analysis may significantly impact the light sterile neutrino interpretation of the short-baseline anomalies.   Comparison of best-fit χ 2 values. The first row gives the difference in χ 2 and degrees of freedom between the null hypothesis (no sterile neutrino) and the best-fit values from global fits to a 3 þ 1 model and 3 þ 1 þ decay model. The second row gives the difference in best-fit χ 2 and degrees of freedom between a 3 þ 1 model and a 3 þ 1 þ decay model. The results in this