Updated Global 3+1 Analysis of Short-BaseLine Neutrino Oscillations

We present the results of an updated fit of short-baseline neutrino oscillation data in the framework of 3+1 active-sterile neutrino mixing. We first consider $\nu_e$ and $\bar\nu_e$ disappearance in the light of the Gallium and reactor anomalies. We discuss the implications of the recent measurement of the reactor $\bar\nu_e$ spectrum in the NEOS experiment, which shifts the allowed regions of the parameter space towards smaller values of $|U_{e4}|^2$. The beta-decay constraints allow us to limit the oscillation length between about 2 cm and 7 m at $3\sigma$ for neutrinos with an energy of 1 MeV. We then consider the global fit of the data in the light of the LSND anomaly, taking into account the constraints from $\nu_e$ and $\nu_\mu$ disappearance experiments, including the recent data of the MINOS and IceCube experiments. The combination of the NEOS constraints on $|U_{e4}|^2$ and the MINOS and IceCube constraints on $|U_{\mu4}|^2$ lead to an unacceptable appearance-disappearance tension which becomes tolerable only in a pragmatic fit which neglects the MiniBooNE low-energy anomaly. The minimization of the global $\chi^2$ in the space of the four mixing parameters $\Delta{m}^2_{41}$, $|U_{e4}|^2$, $|U_{\mu4}|^2$, and $|U_{\tau4}|^2$ leads to three allowed regions with narrow $\Delta{m}^{2}_{41}$ widths at $ \Delta m^2_{41} \approx 1.7 $ (best-fit), 1.3 (at $2\sigma$), 2.4 (at $3\sigma$) eV$^2$. The restrictions of the allowed regions of the mixing parameters with respect to our previous global fits are mainly due to the NEOS constraints. We present a comparison of the allowed regions of the mixing parameters with the sensitivities of ongoing experiments, which show that it is likely that these experiments will determine in a definitive way if the reactor, Gallium and LSND anomalies are due to active-sterile neutrino oscillations or not.


Introduction
Neutrino physics is a powerful probe of new physics beyond the Standard Model. The LSND [1,2], Gallium [3][4][5][6][7] and reactor [8] anomalies are intriguing indications in favor of the existence of light sterile neutrinos connected with low-energy new physics. In order to assess the viability of the light sterile neutrino hypothesis, it is necessary to perform a global fit of neutrino oscillation data which takes into account not only the LSND, Gallium and reactor anomalies, but also the data of many other experiments which constrain activesterile neutrino mixing (see the reviews in Refs. [9][10][11][12]).
In this paper we consider 3+1 active-sterile neutrino mixing, in which the three standard active neutrinos ν e , ν µ , ν τ are mainly composed of three sub-eV massive neutrinos ν 1 , ν 2 , ν 3 and there is a sterile neutrino ν s which is mainly composed of a fourth massive neutrino ν 4 at the eV scale. This is the only allowed four-neutrino mixing scheme after the demise of the 2+2 scheme [13] and the fact that the 1+3 scheme with three massive neutrinos at the eV scale is strongly disfavored by cosmological measurements [14] and by the experimental bounds on neutrinoless double-β decay if massive neutrinos are Majorana particles (see Refs. [15,16]). We do not consider neutrino mixing schemes with more than one sterile neutrino, which are not necessary to explain the current data (see the discussions in Refs. [12,17]).
In the framework of 3+1 active-sterile mixing, short-baseline (SBL) experiments are sensitive only to the oscillations generated by the squared-mass difference ∆m 2 41 ∆m 2 42 ∆m 2 43 1 eV 2 , with ∆m 2 jk ≡ m 2 j − m 2 k , that is much larger than the the solar squaredmass difference ∆m 2 SOL = ∆m 2 21 ≈ 7.4 × 10 −5 eV 2 and the atmospheric squared-mass difference ∆m 2 ATM = |∆m 2 31 | |∆m 2 32 | ≈ 2.5 × 10 −3 eV 2 , which generate the observed solar, atmospheric and long-baseline neutrino oscillations explained by the standard threeneutrino mixing (see Refs. [18,19]). The 3+1 active-sterile mixing scheme is a perturbation of the standard three-neutrino mixing in which the 3×3 unitary mixing matrix U is extended to a 4 × 4 unitary mixing matrix with |U e4 | 2 , |U µ4 | 2 , |U τ 4 | 2 1. The effective oscillation probabilities of the flavor neutrinos in short-baseline experiments are given by [20] P (SBL) αβ δ αβ − sin 2 2ϑ αβ sin 2 ∆m 2 41 L 4E , where α, β = e, µ, τ, s, L is the source-detector distance and E is the neutrino energy. The short-baseline oscillation amplitudes depend only on the absolute values of the elements in the fourth column of the mixing matrix: Hence, the transition probabilities of neutrinos and antineutrinos are equal and it is not possible to measure in short-baseline experiments any CP-violating effect generated by the complex phases in the mixing matrix 1 .
In this paper we update the analysis of short-baseline neutrino oscillation data [7,12,32,33] revising the analysis of the rates measured in reactor neutrino experiments according to Ref. [34] and taking into account the recent measurements of the MINOS [35], IceCube [36], and NEOS [37] experiments. The MINOS and IceCube constraints on ν µ andν µ disappearance are expected [38] to disfavor the low-∆m 2 41 -high-sin 2 2ϑ µµ and the low-∆m 2 41 -high-sin 2 2ϑ eµ parts of the allowed region that we found in our previous analyses [7,12,32,33], as was found in the 3+1 global fit presented in Ref. [39], which updated Ref. [40] with the addition of the IceCube data. The NEOS [37] collaboration measured the spectrum of reactorν e 's at a distance of 24 m and normalized their data to the Daya Bay spectrum [41] measured at the large distance of about 550 m, where shortbaseline oscillations are averaged out. Analyzing this normalized spectrum with shortbaseline oscillations they found the best fit at ∆m 2 41 = 1.73 eV 2 and sin 2 2ϑ ee = 0.05, with a χ 2 which is lower by 6.5 with respect to the standard case of three-neutrino mixing without short-baseline oscillations. This is a 2.1σ indication in favor of short-baseline oscillations and it is intriguing to note that the best-fit value of ∆m 2 41 is close to the best-fit value found in our previous global analysis of short-baseline data [12], ∆m 2 41 = 1.6 eV 2 , albeit with a larger sin 2 2ϑ ee = 0.11. However, as one can see from Table 5 of Ref. [12], the lower bound of the 3σ allowed range of sin 2 2ϑ ee was 0.046, which is below the NEOS best-fit value. Hence, the NEOS data are not incompatible with the global indications of short-baseline oscillations and we expect that their inclusion in the fit will lead to a shift of the allowed region towards smaller values of sin 2 2ϑ ee and, consequently, of |U e4 | 2 .
It is well known [11,12,20,33,[42][43][44][45][46][47][48][49] that the global fits of short-baseline data are affected by the so-called "appearance-disappearance" tension which is present [17] for any number N s of sterile neutrinos in 3+N s mixing schemes which are perturbations of the standard three-neutrino mixing required for the explanation of the observation of solar, atmospheric and long-baseline neutrino oscillations (see Refs. [18,19]). We expect that the inclusion in the global fit of the recent measurements of the MINOS [35], IceCube [36], and NEOS [37] experiment will increase somewhat the appearance-disappearance tension. In Ref. [33] we proposed a "pragmatic approach" in which the appearance-disappearance tension is alleviated by excluding from the global fit the low-energy bins of the MiniBooNE experiment [50,51] which have an anomalous excess of ν e -like events that is under investigation in the MicroBooNE experiment at Fermilab [52]. In this paper we will discuss the effect of MINOS, IceCube and NEOS data on the appearance-disappearance tension and how much it is alleviated in the pragmatic approach.
The plan of the paper is as follows. In section 2 we consider the experimental data on short-baseline ν e andν e disappearance, motivated by the Gallium and reactor anomalies. In section 3 we consider the global fit of appearance and disappearance data, which is motivated by the addition of the LSND anomaly to the Gallium and reactor anomalies. We draw our conclusions in section 4.
2 ν e andν e disappearance In this section we consider only the experimental data on short-baseline ν e andν e disappearance, which include the Gallium neutrino anomaly [3][4][5][6][7] and the reactor antineutrino anomaly [8]. First, we discuss in subsection 2.1 our evaluation of the reactor antineutrino anomaly by considering the updated results of the reactorν e rates measured in several reactor neutrino experiments. In subsection 2.2 we add the constraints of the spectra measured in the old Bugey-3 experiment [53] and in the recent NEOS experiment [37]. Finally, in subsection 2.3 we present our results for the combined fit of reactor and Gallium data and for the global fit of all the ν e andν e disappearance data.

Reactor rates
The reactor neutrino experiments which measured the absolute antineutrino flux that are considered in our analysis 2 are listed in Table 1. For each experiment labeled with the index a, we listed the corresponding four fission fractions f a k , the ratio of measured and predicted rates R exp a , the corresponding relative experimental uncertainty σ exp a , the relative uncertainty σ cor a which is correlated in each group of experiments indicated by the braces, and the relative theoretical uncertainty σ the a which is correlated among all the experiments. The ratios R exp a and the uncertainties σ exp a and σ cor a are the same as those in Ref. [34]. In the following we repeat for convenience 3 their derivation and we explain the derivation of the relative theoretical uncertainty σ the a .
We considered the Krasnoyarsk99-34 experiment [66] that was not considered in Refs. [8,67], by rescaling the value of the corresponding experimental cross section per fission in comparison with the Krasnoyarsk94-57 result. For the long-baseline experiments Chooz [68] and Palo Verde [69], we applied the rescaling in Eq. (2.1) with the ratios R exp a,S given in Ref. [67], divided by the corresponding survival probability P sur caused by ϑ 13 . For Nucifer [70], Daya Bay [41], RENO [71], and Double Chooz 5 we use the ratios provided by the respective experimental collaborations.
The experimental uncertainties and their correlations listed in Table 1 have been obtained from the corresponding experimental papers. In particular: • The Bugey-4 and Rovno91 experiments have a correlated 1.4% uncertainty, because they used the same detector [55].
• The Rovno88 experiments have a correlated 2.2% reactor-related uncertainty [63]. In addition, each of the two groups of integral (Rovno88-1I and Rovno88-2I) and spectral (Rovno88-1S, Rovno88-2S, and Rovno88-3S) measurements have a correlated 3.1% detector-related uncertainty [63]. m from the same reactor of the Krasnoyarsk87-33 experiment [72]. There may be reactor-related uncertainties correlated among the four Krasnoyarsk experiments, but, taking into account the time separations and the absence of any information, we conservatively neglected them.
• Following Ref. [67], we considered the two SRP measurements as uncorrelated, because the two measurements would be incompatible with the correlated uncertainty estimated in Ref. [64].
For each experiment labeled with the index a, the relative theoretical uncertainty σ the a in Table 1 is given by where ρ SH is the covariance matrix of the cross sections per fission of the four fissionable isotopes given in Tab. 3. In this covariance matrix, σ SH f,238 is uncorrelated from the other cross sections per fission and the corresponding uncertainty is that given in Ref. [8] (we neglected the correlation due to the cross section uncertainty, which is of the order of 0.1%). The other three cross sections per fission have been calculated using the Huber [65] antineutrino fluxes which have been obtained by inverting the spectra of the electrons emitted by the β decays of the products of the thermal fission of 235 U, 239 Pu, and 241 Pu which have been measured at ILL in the 80's [73][74][75]. As explained in Ref. [65], the values of the three antineutrino fluxes are correlated. We calculated the uncertainties and correlations of σ SH f,235 , σ SH f,239 , and σ SH f,241 using the information given in Ref. [65]. The square roots of the diagonal elements of the covariance matrix ρ SH in Tab. 3 give the uncertainties of the cross sections per fission reported in Tab. 2. One can see that the uncertainties of σ SH f,235 , σ SH f,239 , and σ SH f,241 are slightly larger than those calculated by Saclay group in Ref. [8].
Let us note that after the discovery of the unexpected "5 MeV bump" in the spectrum of the RENO [76], Double Chooz [77], and Daya Bay [78] experiments it is believed [79,80] that the theoretical uncertainties of the reactor antineutrino fluxes may be larger than those calculated in Refs. [65,81]. However, since there is no well-motivated quantitative estimation of how much the theoretical uncertainties should be increased, we are compelled to use the uncertainties calculated in Refs. [65,81]. Figure 1 shows the experimental ratios as functions of the reactor-detector distance L. The horizontal band shows the average ratio R and its uncertainty, which has been obtained by summing in quadrature the experimental and theoretical uncertainties. Hence, the reactor antineutrino anomaly is at the level of about 2.8σ. The slight difference of the value of R in Eq. (2.3) with respect to our previous estimates in Refs. [7,12] is due to the following three changes in our analysis: 1. The revaluation [34] of the experimental ratios R exp a listed in Table 1. 2. The new treatment of the theoretical uncertainties σ the a according to Eq. (2.2) instead of considering an unrealistic common 2.0% [8].
3. The new data of the Nucifer, Daya Bay, RENO and Double Chooz experiments and the consideration for the first time of the Krasnoyarsk99-34 experiment.
The reactor antineutrino anomaly can be explained in the framework of 3+1 neutrino mixing through neutrino oscillations generated by the effective mixing angle sin 2 2ϑ ee = 4|U e4 | 2 1 − |U e4 | 2 , which determines the survival probability of ν e 's andν e 's according to Eq. (1.1). The result of the fit of the reactor rates are given in the first column of Table 4 and in Fig. 2(a), where we have drawn the allowed regions in the sin 2 2ϑ ee -∆m 2 41 plane. From Fig. 2(a) one can see that the allowed 1σ region 6 in the sin 2 2ϑ ee -∆m 2 41 plane is at a rather low value of ∆m 2 41 , but the allowed regions at 2σ and 3σ extend to higher values of ∆m 2 41 , without an upper bound. The favorite values of the amplitude sin 2 2ϑ ee of ν edisappearance oscillations are around 0.1, but the allowed 3σ region in the sin 2 2ϑ ee -∆m 2 41 plane covers the range 0.0066 sin 2 2ϑ ee 0.28, which corresponds to 0.0017 |U e4 | 2 0.076. Table 4 gives the χ 2 difference ∆χ 2 NO between the χ 2 of no oscillations and χ 2 min , and the corresponding number of σ's for the two degrees of freedom corresponding to the two fitted parameters ∆m 2 41 and sin 2 2ϑ ee . The case of no oscillations turns out to be disfavored at the level of 3.2σ.

Reactor spectra
In our previous analyses [7,12,32,33] we considered the ratio of the spectra measured at 40 m and 15 m from the source in the Bugey-3 experiment [53]. These data provide robust information on short-baselineν e disappearance, which is independent of any theoretical calculation of the spectrum and of the solution of the "5 MeV bump" problem mentioned in subsection 2.1.
In this paper we add the constraints obtained in the recent NEOS experiment by taking into account the χ 2 corresponding to Fig. 4 of Ref. [37], which has been kindly provided to us by the NEOS Collaboration 7 . The NEOS constraints are mostly independent of theoretical calculations of the spectrum and of the solution of the "5 MeV bump" problem, because the NEOS χ 2 has been obtained by fitting the NEOS spectrum normalized to the Daya Bay spectrum [41] measured at the large distance of about 550 m, where the shortbaseline oscillations due to ∆m 2 41 are averaged out. A small dependence on the theoretical calculation of the spectrum [65,81] comes from the corrections due to the differences of the fission fractions of the NEOS and Daya Bay reactors [37,41] and a small dependence on the "5 MeV bump" problem comes from a possible dependence of the "5 MeV bump" on 6 In all the paper we consider allowed regions at 1σ, 2σ, and 3σ, which correspond, respectively, to 68.27%, 95.45%, and 99.73% confidence level. The allowed regions in two-parameter planes are drawn considering two degrees of freedom, which correspond, respectively, to ∆χ 2 = 2.30, 6.18, and 11.83 from the minimum χ 2 min . 7 NEOS Collaboration, Private Communication.
the different fission fractions of NEOS and Daya Bay [82]. We neglect these possible small effects.
The results of the fit of the Bugey-3 and NEOS spectra are given in the second column of Table 4 and in Fig. 2(b), where one can see that the NEOS constraints are dominating. There are closed allowed islands at 2σ which are determined mainly by the NEOS data and the best-fit values of the oscillation parameters in Table 4 correspond to the best fit reported in Ref. [37]. Hence, the NEOS constraints can be interpreted as a weak indication in favor of short-baseline oscillations which may be compatible with the reactor antineutrino anomaly based on the reactor rates discussed in subsection 2.1. This is confirmed by the disfavoring of the case of no oscillations at the level of 2.1σ, as shown in Table 4.
The third column of Table 4 and Fig. 3(a) show the results of the combined fit of the rate and spectral data of reactor antineutrino experiments. As reported in Table 4, the combined fit disfavors the case of no oscillations at the level of 2.9σ, which is about the same level obtained from the analysis of the reactor rates alone. Hence, the NEOS constraints do not exclude the reactor antineutrino anomaly. However, in spite of the weak NEOS indication in favor of short-baseline oscillations discussed above, the statistical significance of the anomaly does not increase by including the NEOS data because there is a mild tension with the reactor rates which is illustrated by the 2σ contours in Fig. 3(a). Indeed, the rates-spectra parameter goodness-of-fit is only 2% (∆χ 2 /NDF = 8.3/2).

Global ν e andν e disappearance
In this subsection we discuss the combination of the reactor data with the data of the Gallium neutrino anomaly, other ν e andν e disappearance data and the β-decay constraints of the Mainz [83] and Troitsk [84,85] experiments.
The fourth column of Table 4 and Fig. 3(b) show the results of the combined fit of reactor and Gallium data. Since both sets of data indicate short-baseline ν e andν e disappearance, the statistical significance of active-sterile neutrino oscillations increases to 3.6σ and the 3σ allowed regions in the sin 2 2ϑ ee -∆m 2 41 plane are confined to 0.010 sin 2 2ϑ ee 0.30 and ∆m 2 41 0.35 eV 2 . Besides the reactor and Gallium data, short-baseline ν e disappearance 8 is constrained by solar and KamLAND neutrino data [7,[87][88][89][90], by the KARMEN [91] and LSND [92] ν e + 12 C → 12 N g.s. + e − scattering data [46,93] and by the T2K near detector constraints [94].
We updated our 2012 solar+KamLAND constraint [7] by including the latest solar data: the new results from the fourth phase of the Super-Kamiokande experiment [95] and the final results of Borexino phase-I [96]. We also updated the KamLAND data analysis by using the Saclay+Huber cross sections per fission [34]. Finally, we used the updated value of ϑ 13 in the 2016 Review of Particle Physics [97]. We obtained the new marginal ∆χ 2 shown in Fig. 4, where it is confronted with the old one obtained in Ref. [7]. The new solar+KamLAND constraint is weaker than the 2012 one because of the larger Saclay+Huber reactor rate prediction used in the analysis of KamLAND data and because the new value of ϑ 13 is smaller than that in 2012.
The results of the combined analysis of all ν e andν e disappearance data are shown in the fifth column of Table 4 and Fig. 5(a). Since the analysis of solar+KamLAND, ν e -12 C, and T2K data do not show any indication of short-baseline ν e disappearance, the combination with the reactor and Gallium data shifts the allowed regions in the sin 2 2ϑ ee -∆m 2 41 plane in Fig. 5(a) to smaller values of sin 2 2ϑ ee with respect to Fig. 3 can be constrained with the data of β-decay experiments (see Ref. [12]). As in Ref. [32], we use the β-decay constraints of the Mainz [83] and Troitsk [84,85] experiments, which give the allowed regions in the sin 2 2ϑ ee -∆m 2 41 plane shown in Fig. 5(b). One can see that the allowed regions are confined to the range For the oscillation length L osc 41 = 4πE/∆m 2 41 we have This is a range of oscillation lengths which can be explored in a model independent way in the new short-baseline reactor neutrino experiments (DANSS [98], Neutrino-4 [54], PROSPECT [99], SoLid [100], STEREO [101]) and in the SOX [102] and BEST [103] radioactive source experiments by measuring the reactor antineutrino rate as a function of distance. However, they will need to be sensitive to small oscillations with the amplitude 0.022 (0.0050) sin 2 2ϑ ee 0.19 (0.23) at 2σ (3σ).
(2.6) Figure 6(a) shows the sensitivities of the short-baseline reactor antineutrino experiments DANSS [98], Neutrino-4 [104], PROSPECT [99], SoLid [105], and STEREO [106] in comparison with the allowed regions in the sin 2 2ϑ ee -∆m 2 41 plane in Fig. 5(b). One can see that they will cover most of the allowed regions for ∆m 2 41 10 eV 2 and not too small sin 2 2ϑ ee . Figure 6(b) shows the sensitivities of the CeSOX [102] and BEST [103] source experiments, of IsoDAR@KamLAND [107] and C-ADS [108], and of the KATRIN [109]) electron neutrino mass experiment 9 . The source experiments will cover the large-sin 2 2ϑ ee parts of the allowed regions, the IsoDAR@KamLAND and C-ADS experiments can cover almost all the allowed regions, except the large-∆m 2 41 part and the small-sin 2 2ϑ ee -small-∆m 2 41 parts, and KATRIN will cover the large-∆m 2 41 part. Hence, there are favorable perspectives for a definitive solution of the short-baseline (−) ν e disappearance problem in the near future.

Fits of appearance and disappearance data
In this section we present the results of 3+1 fits of short-baseline neutrino oscillation data which include ν µ → ν e andν µ →ν e appearance data and ν µ andν µ disappearance data, in addition to the ν e andν e disappearance data considered in section 2. Our fits are based on a χ 2 analysis in the four-dimensional space of the mixing parameters ∆m 2 41 , |U e4 | 2 , |U µ4 | 2 , and |U τ 4 | 2 .
There is no indication in favor of short-baseline ν µ andν µ disappearance from any experiment. Therefore, the current ν µ andν µ disappearance data lead to constraints on |U µ4 | 2 . We consider the constraints obtained from the CDHSW experiment [120], from the analysis [121] of the data of atmospheric neutrino oscillation experiments, from the analysis of the SciBooNE-MiniBooNE neutrino [122] and antineutrino [123] data, which were included in our previous fits [12,33,38]. In addition, we take into account the recent data of the MINOS [35] and IceCube [36] experiments. The MINOS constraint was easily included in our numerical computation by using the ROOT program in the MINOS data release 10 , which computes the χ 2 for input values of the 3+1 mixing parameters. On the other hand, we had to calculate the IceCube χ 2 , as described in the following subsection 3.1.

Analysis of IceCube data
The IceCube detector measures the incoming (anti)muons generated by the interaction of atmospheric muon (anti)neutrinos with the surrounding earth and ice, as a function of the neutrino energy and of the zenith angle. For high-energy, up-going atmospheric neutrinos that reach the detector after having crossed the Earth, the ratio L/E is of the same order of that in SBL experiments. In this case, the oscillations arising from the usual atmospheric and solar squared mass differences have a very long wavelength and can be neglected, but the SBL squared mass difference ∆m 2 41 plays an active role. The sterile neutrino influence on the observed flux is given by the matter effects that modify the oscillation patterns inside the Earth. This happens because the matter potential is different for the different active neutrino flavors, for which the charged and neutral current interactions are not the same [124], and there is no potential for the sterile neutrinos.
We use the 20,145 released IceCube events in the approximate energy range between 320 GeV and 20 TeV, detected over 343.7 days in the 86-strings configuration [125] for constraining the active-sterile mixing parameters. The 99.9% of the IceCube events is expected to come from neutrino-induced muon events, where the neutrinos originate from the decays of atmospheric pions and kaons. The contribution from charmed meson decays is negligible [125,126].
The calculation of the χ 2 contribution from IceCube is divided into three parts: the calculation of the theoretical flux for each set of mixing parameters, for which one needs to propagate the atmospheric neutrinos through the Earth, the estimate of the expected number of events in the detector, for which we use the IceCube Monte Carlo data, and finally the computation of the χ 2 , obtained comparing theoretical and observed events. For all these parts we use the data 11 and we follow the prescriptions presented in Ref. [36].
To obtain the predicted neutrino flux at the detector, we use the ν-SQuIDS code 12 , a C++ package based on the Simple Quantum Integro-Differential Solver (SQuIDS) 13 [127], that contains all the necessary tools to numerically solve the master equation that rules the neutrino evolution in the Earth [124].
The initial flux we consider is the unoscillated HKKM flux [128][129][130][131] with the H3a knee correction [132], that we use for obtaining the initial spectrum of neutrinos from pion and kaon decays. This model is usually referred to as the "Honda-Gaisser" model. We do not employ here the other six atmospheric flux variants that are considered in Ref. [36], but we tested that our results do not change significantly if another model is used instead of the Honda-Gaisser one. Since our analysis is based not only on the IceCube data, our final result would be almost unaffected.
The unoscillated flux is propagated inside the Earth with the ν-SQuIDS code, which uses the Preliminary Reference Earth Model [133] for the inner density profile of the Earth. For the neutrino-matter interactions, the charged current cross section is dominated by deep inelastic scattering, which involves neutrino-nucleon scattering. The main uncertainty in this case is in the parton distribution functions. In the ν-SQuIDS code, the perturbative QCD calculation in Refs. [134] are used for the neutrino-nucleon cross-section calculation. We do not treat the uncertainties on the Earth density profile and on the deep inelastic scattering cross section.
The full expression for the (anti)neutrino flux at the detector is given by [36] (see also Refs. [135,136] (3.1) Here, θ is the zenith angle and E ν(ν) the energy of the incoming (anti)neutrino, while φ is the oscillated (anti)neutrino flux from pion (kaon) decays. The free parameters in the above equation are: the neutrino and antineutrino flux normalizations, N ν 0 and Nν 0 ; the pion-to-kaon ratio, R π/K ; the spectral index correction, ∆γ. These are treated as continuous nuisance parameters in our analysis, as explained in Ref. [36]. The pivot energy E m is fixed to be approximately near the median of the energy distribution of the measured events, being E m = 2 TeV.
The function F(δ) parameterizes the atmospheric density uncertainties. This function is assumed to be linear and it is obtained by imposing the AIRS constraints on the 11 http://icecube.wisc.edu/science/data/IC86-sterile-neutrino/ 12 http://github.com/Arguelles/nuSQuIDS 13 http://github.com/jsalvado/SQuIDS atmospheric temperature 14 . The expression reads as [135]: where E 0 = 360 GeV, E 1 = 11.279 TeV, κ = 200 and cos θ 0 = 0.4. The parameter δ represents the last one of our nuisance parameters. The theoretical flux is converted into a number of expected events using the Monte Carlo (MC) data released by the IceCube collaboration [36]. The MC data are needed to model the detector capabilities to measure the incoming events as a function of the real energy and zenith angle of the muon (anti)neutrino, and of the corresponding quantities for the reconstructed (anti)muon event. For each combination of mixing and nuisance parameters, we use the MC data to convert the obtained theoretical flux into the expected number of events that we compare with the data as explained below. Since IceCube cannot distinguish a muon from an antimuon, neutrinos and antineutrinos events are summed up together. It is however important to treat properly both the components, since the matter oscillation patterns for neutrinos and antineutrinos are different, with the consequence that the disappearance of neutrinos and antineutrinos is not the same.
We build our χ 2 using a binned Poissonian likelihood, written as where n i represents the number of observed events in the bin i and µ i (θ) the corresponding number of expected events as a function of the model parameters θ, that includes both mixing and nuisance parameters. Following Ref. [36], we consider a grid with 10 logarithmic bins in the reconstructed energy, with 400 GeV ≤ E reco µ(μ) ≤ 20 TeV, and 21 linear bins for the cosine of the reconstructed zenith angle, in the range −1 ≤ cos θ reco µ(μ) ≤ 0.2. For each combination of the mixing parameters, we minimize the χ 2 over the five nuisance parameters described above (N ν 0 , Nν 0 , R π/K , ∆γ, δ). We adopt a standard Nelder-Mead algorithm for the minimization [137]. It is important to note that for each point in the mixing parameter space we needed to minimize independently over the nuisance parameters. We checked that the preferred values of the nuisance parameters do not vary significantly outside the adopted Gaussian priors [36].
We show in Fig. 7 the comparison of the official IceCube 90% and 99% CL exclusion curves in the sin 2 2ϑ µµ -∆m 2 41 plane for |U e4 | 2 = |U τ 4 | 2 = 0 [36] with our results. In our analysis of IceCube data we do not vary the efficiency of the Digital Optical Modules because we do not have sufficient information. Despite this fact, one can see from Fig. 7 that the results of our analysis are in good agreement with those of the IceCube collaboration. Moreover, we emphasize that the IceCube data are just one of the datasets in our global analyses, and small differences in the IceCube analysis do not play a significant role when computing the global fit.
Since the calculation of the χ 2 given a set of mixing parameters is a highly demanding computational task, it is impossible to directly include the code that calculates the χ 2 of the IceCube data in our complete fitting routine without slowing it down in an unacceptable way. Therefore, we adopted the following method. Since we are more interested in scanning the region near the expected 3+1 mixing best-fit, we employed the results of the 3+1 fit of SBL data without IceCube in order to generate with a Markov Chain Monte Carlo (MCMC) 3,000 random points whose distribution covers an area of the parameter space around the expected best-fit region. In order to cover the rest of the full four-dimensional parameter space (∆m 2 41 ; |U e4 | 2 , |U µ4 | 2 , |U τ 4 | 2 ), we generated uniformly 21,000 more random points. We end up with a set P of 24,000 points for which we computed the IceCube χ 2 in an affordable time. In the complete fitting routine, we computed the IceCube contribution to the χ 2 in each point in the full four-dimensional parameter space with a linear interpolation of the χ 2 's of the nearest points in the set P obtained with a Delaunay triangulation.

Fit of the 2016 data set without MINOS and IceCube
In this subsection we present the results of the 3+1 global fit "Glo16A" of the appearance and disappearance SBL data available in 2016 15 , except MINOS [35] and IceCube [36], which will be added in subsection 3.3 in order to clarify their effects on the results of the analysis. In the Glo16A fit we also do not take into account the NEOS [37] data which have been available to us in the beginning of 2017 and will be considered in subsection 3.4.
The results of the Glo16A fit are shown by the first column of Tab. 5, by Fig. 8, and by the solid purple curves in Fig. 9, which gives the marginal ∆χ 2 as a function of the mixing parameters ∆m 2 41 , |U e4 | 2 , and |U µ4 | 2 , from which one can obtain the corresponding marginal allowed intervals at different confidence levels.
The global goodness of fit of 4.8% is acceptable, but there is a relevant appearancedisappearance tension quantified by a parameter goodness of fit of 0.13%. If one is willing to accept such appearance-disappearance tension, one can adopt the allowed regions of the oscillation parameters shown in Fig. 8.
The Glo16A fit is an update of the GLO fit presented in Ref. [12], with a similar set of data. It can also be compared with the global fit in Ref. [40], where a similar set of data was considered. With respect to Ref. [40], we find larger allowed regions for ∆m 2 2 eV 2 and we do not have the allowed region at ∆m 2 ≈ 6 eV 2 found in Ref. [40] at 99% CL. However, there is an approximate agreement of our results with those of Ref. [40], with a remarkable closeness of the best-fit point in the mixing parameter space.

Effects of MINOS and IceCube
In this subsection we present the 3+1 global fit "Glo16B" with the addition of the 2016 data of the MINOS [35] and IceCube [36] experiments. The results are shown by the 15 We consider all the νe andνe disappearance data discussed in section 2, with the exceptions of the T2K near detector constraints [94] on sin 2 2ϑee, which unfortunately cannot be included in the global fit because they have been obtained under the assumption |Uµ4| 2 = 0, and of the Mainz [83] and Troitsk [84,85] β-decay constraints, which are not needed because the value of ∆m 2 41 is constrained within a few eV 2 by the combination of appearance and disappearance data. second column of Tab. 5, by Fig. 10, and by the solid blue curves in Fig. 9.
Comparing Fig. 10(c) with Fig. 8(c), one can see that the addition of the MINOS and IceCube data leads to the exclusion of the low-∆m 2 41 -high-sin 2 2ϑ µµ part of the allowed region, as expected (see the discussion in section 1). On the other hand, the high-∆m 2 41low-sin 2 2ϑ µµ part of the allowed region is practically unaffected by the MINOS and IceCube constraints. As a consequence, also the low-∆m 2 41 -high-sin 2 2ϑ eµ part of the allowed region in Fig. 8(a) is excluded in Fig. 10(a), whereas the high-∆m 2 41 -low-sin 2 2ϑ eµ part of the allowed region is practically unaffected.
From Tab. 5 one can see that including the MINOS and IceCube data increases the appearance-disappearance tension by lowering the parameter goodness of fit from 0.13% to 0.075%. This is a consequence of the reduction of the upper limit of the allowed range of |U µ4 | 2 in the Glo16B fit with respect to the Glo16A fit shown in Fig. 9(c). Figure 11(a) shows the effect of adding to the data set of the Glo16A fit the MINOS and IceCube data separately and together. One can see that the IceCube data are slightly more effective than the MINOS data in reducing the low-∆m 2 41 -high-sin 2 2ϑ eµ part of the allowed region.
The MINOS and IceCube data give information not only on the 3+1 mixing parameters ∆m 2 41 , |U e4 | 2 , and |U µ4 | 2 that we have considered so far, but also on |U τ 4 | 2 . The sensitivity to |U τ 4 | 2 is due in MINOS to the neutral-current event sample [35] and in IceCube to the matter effects for high-energy neutrinos propagating in the Earth, which depend on all the elements of the mixing matrix [138][139][140][141][142][143][144][145]. Limits on the value of |U τ 4 | 2 have been obtained in the analyses of the atmospheric neutrino data of the Super-Kamiokande [146] and IceCube DeepCore [147] experiments, in the analysis of the data of the MINOS experiment [35,148,149], and in the phenomenological fits in Refs. [39,49]. There is also a bound on sin 2 2ϑ µτ = 4|U µ4 | 2 |U τ 4 | 2 given by the absence of a 3+1 excess of ν µ → ν τ oscillations in the OPERA experiment [150]. Figure 9(d) shows the marginal ∆χ 2 as a function of |U τ 4 | 2 in our Glo16B fit, from which one can see that we obtain the stringent upper bound |U τ 4 | 2 0.022 (0.071) at 2σ (3σ). (3.4) At 90% CL we have |U τ 4 | 2 0.014 and ϑ 34 7.4 • in the common parameterization of the 4 × 4 unitary mixing matrix used in Ref. [39]. This bound on ϑ 34 is about the same as that reported in Ref. [39] for ∆m 2 41 ≈ 6 eV 2 . However, we do not find a 90% CL allowed region of the mixing parameters at ∆m 2 41 ≈ 6 eV 2 and our bound on ϑ 34 applies for any value of ∆m 2 41 . Figure 11(b) shows the correlated bounds on |U τ 4 | 2 and ∆m 2 41 that we obtain considering the MINOS and IceCube data separately and together. One can see that the IceCube data give more stringent constraints on |U τ 4 | 2 than the MINOS data for ∆m 2 41 1.5 eV 2 .

Effects of NEOS
We finally consider also the NEOS [37] data and obtain the 3+1 global fit "Glo17" which includes all data available so far in 2017. The results are shown by the third column of Tab. 5, by Fig. 12, and by the solid orange curves in Fig. 9.
Comparing Fig. 12 with Fig. 10, it is evident that the inclusion of the NEOS constraints has a dramatic effect on the allowed regions, leading to the fragmentation of the allowed region in three islands with narrow ∆m 2 41 widths. The best-fit island is at ∆m 2 41 ≈ 1.7 eV 2 . There is an island allowed at 2σ at ∆m 2 41 ≈ 1.3 eV 2 , and an island allowed at 3σ at ∆m 2 41 ≈ 2.4 eV 2 . Moreover, the NEOS constraints shifts the 3σ allowed range of |U e4 | 2 from 0.014 − 0.051 in the Glo16B fit to 0.011 − 0.032 in the Glo17 fit, as shown in Fig. 9. Therefore, the appearance-disappearance tension is increased, as shown by the 0.019% parameter goodness of fit in Tab. 5. Since this low value of the appearance-disappearance parameter goodness of fit is hardly acceptable, we are led to consider, in the next subsection, the "pragmatic approach" proposed in Ref. [33].

Pragmatic fit
In this section we consider the "pragmatic approach" [33] in which the low-energy bins of the MiniBooNE experiment [50,51] which have an anomalous excess of (−) ν e -like events are omitted from the global fit. As shown in Fig. 1b of Ref. [38], the region allowed by the appearance data shifts towards larger values of ∆m 2 41 and smaller values of sin 2 2ϑ eµ when the MiniBooNE low-energy bins are omitted from the fit. As a result, the overlap of the appearance and disappearance allowed regions increases, relieving the appearancedisappearance tension.
One can question the scientific correctness of the data selection in the pragmatic approach, but we note that the MiniBooNE low-energy excess is widely considered to be suspicious 16 because of the large background. Some of this background can be due to photon events which are indistinguishable from ν e event in a liquid-scintillator detector. The suspicion that this photon background may be responsible for the MiniBooNE low-energy excess motivated the realization of the MicroBooNE experiment at Fermilab [52], which is able to distinguish between photon and (−) ν e events by using a Liquid Argon Time Projection Chamber (LArTPC). Waiting for the results of this experiment, we think that it is reasonable to adopt the pragmatic approach of omitting from the global fit the MiniBooNE low-energy data.
The results of the pragmatic 3+1 global fit "PrGlo17", which includes the MINOS, IceCube and NEOS data, are shown by the fourth column of Tab. 5, by Fig. 13, and by the dashed red curves in Fig. 9.
From Tab. 5 one can see that, as expected, the exclusion from the fit of the MiniBooNE low-energy data leads to an increase of the parameter goodness of fit from the unacceptable 0.019% of the Glo17 fit to the acceptable 2.7% of the PrGlo17 fit. There is still a mild appearance-disappearance tension, but the tolerable value of parameter goodness of fit leads us to consider the PrGlo17 fit as acceptable.
Comparing the allowed regions of the oscillation parameters in Fig. 13 for the PrGlo17 fit with those in Fig. 12 for the Glo17 fit and the corresponding marginal ∆χ 2 curves in Fig. 9, one can see that the differences are small. As a consequence of the larger overlap of the regions allowed by the fits of appearance and disappearance data, the PrGlo17 fit has a minimum χ 2 significantly smaller than the Glo17 fit, which leads to an increased preference of the best-fit island at ∆m 2 41 ≈ 1.7 eV 2 , to a small shrink of the island at ∆m 2 41 ≈ 1.3 eV 2 , and at a significant reduction of the island at ∆m 2 41 ≈ 2.4 eV 2 (the corresponding 3σ interval for one degree of freedom allowed by the marginal ∆χ 2 in Fig. 9(a) disappears). Table 6 gives the marginal allowed intervals of the mixing parameters |U e4 | 2 , |U µ4 | 2 , and |U τ 4 | 2 . The stringent upper bounds on |U τ 4 | 2 slightly improve those found in the Glo16B fit (see Eq. (3.4) and Fig. 9(d)). At 90% CL we have |U τ 4 | 2 0.011 and ϑ 34 6 • .

Conclusions
In this paper we updated the global fit of short-baseline neutrino oscillation data in the framework of 3+1 active-sterile neutrino mixing [7,12,32,33].
We considered first, in section 2, the data on ν e andν e disappearance which include the Gallium neutrino anomaly data [3][4][5][6][7] and the reactor antineutrino anomaly data [8]. The resulting allowed region in the sin 2 2ϑ ee -∆m 2 41 plane is rather wide, as shown in Fig. 5(a), but it is smaller than that found in our previous analysis [7], mainly as a result of the constraints given by the recent NEOS [37] experiment. The allowed region obtained with neutrino oscillation data alone has no upper bound for ∆m 2 41 , but it can be limited [32] using the constraints found in the Mainz [83] and Troitsk [84,85] β-decay experiments, as shown in Fig. 5(b). We found the upper limit ∆m 2 41 148 eV 2 at 3σ. Hence, as shown in Fig. 6, the ongoing reactor, source and β-decay experiments can clarify in a definitive way the existence of short-baseline (−) ν e disappearance due to active-sterile neutrino mixing. We presented also, in section 3, the results of global fits of all the available ν µ disappearance data, in addition to the (−) ν e disappearance data considered in section 2. We discussed the effects on the global fits of the recent data of the MINOS [35], IceCube [36], and NEOS [37] experiments. As expected, the MINOS, IceCube and NEOS data aggravate the appearance-disappearance tension, which becomes tolerable only in the pragmatic PrGlo17 fit discussed in subsection 3.5, which is our recommended result.
We found that, as expected [38,39], the MINOS and IceCube constraints on (−) ν µ disappearance disfavor the low-∆m 2 41 -high-sin 2 2ϑ µµ and the low-∆m 2 41 -high-sin 2 2ϑ eµ parts of the allowed region. The addition of the NEOS data has the more dramatic effect of reducing the allowed region to three islands with narrow ∆m 2 41 widths and 0.00048 sin 2 2ϑ eµ 0.0020 at 3σ. The best-fit island is at ∆m 2 41 ≈ 1.7 eV 2 . There is an island allowed at 2σ at ∆m 2 41 ≈ 1.3 eV 2 , and an island allowed at 3σ at ∆m 2 41 ≈ 2.4 eV 2 . However, as illustrated in Fig. 14, the ongoing and planned experiments have the possibility to cover all the allowed regions of the mixing parameters and we expect that they will reach in a few years a definitive conclusion on the existence of the short-baseline oscillations indicated by the LSND experiment and by the Gallium and reactor neutrino anomalies.
We did not consider the problem of the cosmological bounds on active-sterile neutrino mixing [14], which most likely must be solved with a non-standard effect as a large lepton asymmetry [167][168][169][170] or secret interactions of the sterile neutrino mediated by a massive vector or pseudoscalar boson [171][172][173][174][175][176][177], which suppress the thermalization of the sterile neutrino in the early Universe.
In conclusion, this paper gives information on what are the regions of the parameter space of 3+1 neutrino mixing which must be explored by new experiments in order to check the indications given by the LSND, Gallium and reactor anomalies. Let us emphasize the importance of an experimental confirmation of these oscillations, that would imply the existence of light sterile neutrinos. These are new particles with properties outside the realm of the Standard Model and their discovery would open a prodigious window on new low-energy physics. is the ratio of measured and predicted rates, σ exp a is the corresponding relative experimental uncertainty, σ cor a is the relative systematic uncertainty which is correlated in each group of experiments indicated by the braces, σ the a is the relative theoretical uncertainty which is correlated among all the experiments, and L a is the source-detector distance.   Table 4. Results of the fits of ν e andν e disappearance data: minimum χ 2 (χ 2 min ), number of degrees of freedom (NDF), goodness of fit (GoF), best fit values of ∆m 2 41 , sin 2 2ϑ ee , and |U e4 | 2 , χ 2 difference ∆χ 2 NO between the χ 2 of no oscillations and χ 2 min , and the resulting number of σ's (nσ NO ) for two degrees of freedom corresponding to two fitted parameters (∆m 2 41 and sin 2 2ϑ ee ). The columns correspond to the fits of the data of reactor rates (Rea:Rat), reactor spectra (Rea:Spe), reactor rates and spectra (Rea:Rat+Spe), reactor and Gallium data (Rea+Gal), ν e andν e disappearance data (ν e Dis), ν e andν e disappearance data and β decay constraints (ν e Dis+β).  Table 5.
Results of the 3+1 global Glo16A, Glo16B, Glo17, and PrGlo17 fits of SBL data discussed, respectively, in subsections 3.2, 3.3, 3.4, and 3.5. The first group of rows gives: the minimum χ 2 (χ 2 min ), the number of degrees of freedom (NDF), the goodness of fit (GoF), the best fit values of the mixing parameters ∆m 2 41 , |U e4 | 2 , |U µ4 | 2 , and of the oscillation amplitudes sin 2 2ϑ eµ , sin 2 2ϑ ee , sin 2 2ϑ µµ . The second group of rows gives the χ 2 difference ∆χ 2 NO between the χ 2 of no oscillations and χ 2 min and the resulting number of σ's (nσ NO ) for NDF NO degrees of freedom corresponding to the number of fitted parameters. The third and fourth group of rows give, respectively, the results of different 3+1 fits of appearance (App) and disappearance (Dis) data. The fifth group of rows gives the results for the appearance-disappearance parameter goodness of fit [178]: the χ 2 difference ∆χ 2 PG and the resulting goodness of fit GoF PG for NDF PG degrees of freedom.    Figure 2. Allowed regions in the sin 2 2ϑ ee -∆m 2 41 plane and marginal ∆χ 2 's for sin 2 2ϑ ee and ∆m 2 41 obtained from: (a) the combined fit of the rates of the reactor neutrino experiments in Tab. 1; (b) the combined fit of the spectra of Bugey-3 [53] and NEOS [37] reactor antineutrino experiments; The best-fit points corresponding to χ 2 min in Table 4 are indicated by crosses.  Table 4 are indicated by crosses.   Figure 5. Allowed regions in the sin 2 2ϑ ee -∆m 2 41 plane and marginal ∆χ 2 's for sin 2 2ϑ ee and ∆m 2 41 obtained from: (a) the combined fit of ν e andν e disappearance data; (b) the combined fit of ν e andν e disappearance data and the β-decay constraints of the Mainz [83] and Troitsk [84,85] experiments. The best-fit points corresponding to χ 2 min in Table 4 are indicated by crosses. (c) planes obtained in the 3+1 global fit "Glo16A" of the 2016 SBL data without the MINOS [35] and IceCube [36] data. There is a comparison with the 3σ allowed regions obtained from  (d) Figure 9. Marginal ∆χ 2 = χ 2 − χ 2 min as a function of the mixing parameters ∆m 2 41 (a), |U e4 | 2 (b), |U µ4 | 2 (c), and |U τ 4 | 2 (d). The black horizontal lines show the ∆χ 2 for one degree of freedom corresponding to the indicated confidence level (CL).   Figure 11. Comparison of (a) the 3σ allowed regions in the sin 2 2ϑ eµ -∆m 2 41 plane and (b) the 2σ allowed regions in the |U τ 4 | 2 -∆m 2 41 plane obtained by adding to the data set of the Glo16A fit the MINOS and IceCube data separately and together, and by adding also the NEOS data.  Figure 13. Allowed regions in the sin 2 2ϑ eµ -∆m 2 41 (a), sin 2 2ϑ ee -∆m 2 41 (b), and sin 2 2ϑ µµ -∆m 2 41 (c), planes obtained in the pragmatic 3+1 global fit "PrGlo17" of SBL data. There is a comparison with the 3σ allowed regions obtained from   Figure 14. Sensitivities of future experiments compared with the PrGlo17 allowed regions of Fig. 13.