Sterile Neutrinos or Flux Uncertainties? - Status of the Reactor Anti-Neutrino Anomaly

The $\sim 3\sigma$ discrepancy between the predicted and observed reactor anti-neutrino flux, known as the reactor anti-neutrino anomaly, continues to intrigue. The recent discovery of an unexpected bump in the reactor anti-neutrino spectrum, as well as indications that the flux deficit is different for different fission isotopes seems to disfavour the explanation of the anomaly in terms of sterile neutrino oscillations. We critically review this conclusion in view of all available data on electron (anti)neutrino disappearance. We find that the sterile neutrino hypothesis cannot be rejected based on global data and is only mildly disfavored compared to an individual rescaling of neutrino fluxes from different fission isotopes. The main reason for this is the presence of spectral features in recent data from the NEOS and DANSS experiments. If state-of-the-art predictions for reactor fluxes are taken at face value, sterile neutrino oscillations allow a consistent description of global data with a significance close to $3\sigma$ relative to the no-oscillation case. Even if reactor fluxes and spectra are left free in the fit, a $2\sigma$ hint in favour of sterile neutrinos remains, with allowed parameter regions consistent with an explanation of the anomaly in terms of oscillations.


Introduction
Calculations of the anti-neutrino fluxes emmitted from nuclear reactors performed in 2011 [1,2] have led to an increased flux prediction compared to previous estimates [3][4][5][6]. This implies a deficit in the observed reactor neutrino measurements compared to predictions, which is known under the name "reactor anti-neutrino anomaly" (RAA) [7]. Using the published systematic errors on the flux predictions, the significance of this anomaly is around 2.8σ. The anomaly can be explained by oscillations of electron anti-neutrinos into a light sterile neutrinos with a mass-squared difference of order 1 eV 2 [7]. For recent reviews on reactor neutrino flux calculations, and on possible caveats with these calculations that could account for the anomaly, see refs. [8][9][10][11][12]. The neutrino oscillation hypothesis is supported by an independent anomaly, namely a similar deficit of neutrinos in experiments using an intense radioative source in conjunction with gallium-based radiochemical detectors. This deficit is usually referred to as the "Gallium anomaly" [13,14].
In this work we re-consider the sterile neutrino interpretation of the reactor and gallium anomalies and update our analysis of this tension from refs. [15,16] (see also ref. [17,18]) in the light of the following recent experimental developments: 1. Precise measurements of the reactor anti-neutrino spectrum by modern experiments [19][20][21] have revealed a spectral feature ("bump") at neutrino energies around E ν ∼5 MeV, which is not predicted by the theoretical flux calculations. A compilation of results and a possible explanation in terms of detector energy scale non-linearity are presented in ref. [22]. The author of ref. [23] concludes from the data that the likely source of this feature is the anti-neutrino flux from 235 U fission. More discussions about possible origins of the bump can be found in ref. [24]. While the origin of the bump is under debate and sheds some doubt on the reliability of flux calculations (or their error estimates), its presence cannot explain the RAA.
2. Daya Bay [25], as well as the short-baseline reactor experiments NEOS [26] and DANSS [27,28] have presented new limits on sterile neutrino oscillations in the relevant parameter region. These new analyses rely on relative comparisons of measured spectra at different baselines and are therefore independent of flux predictions. While they find no clear evidence for oscillations, their observed spectra show some distortions which are consistent with the presence of a sterile neutrino in certain regions of parameter space. We will quantify the impact of these new results in relation to the previous RAA.
3. Using the time evolution of the observed anti-neutrino rate and the knowledge of the isotopic composition of the reactor cores, the Daya Bay collaboration was able to determine the individual anti-neutrino fluxes from the four most important fissible isotopes 235 U, 238 U, 239 Pu, 241 Pu [29]. Their results suggest that the flux from 235 U is the main source for the anomaly, while the one from 239 Pu is consistent with the prediction. Fluxes from 238 U and 241 Pu are numerically less important. Such a result would disfavour the sterile neutrino hypothesis, which predicts equal suppression of the fluxes from all isotopes. Below, we will quantify to what extent the Daya Bay measurement excludes a sterile neutrino explanation of the RAA in the context of global data.
The hypothesis that the anomaly is due to a mis-prediction of the 235 U flux has already received support in the global analysis from ref. [30], predating the Daya Bay results of ref. [29], and also in ref. [31], which includes the new Daya Bay data. On the other hand, the authors of ref. [12] demonstrate that, when comparing the data to a flux prediction based on nuclear data tables rather than measured beta decay spectra, the anomaly is of similar magnitude for all isotopes and therefore consistent with the sterile neutrino hypothesis. In ref. [32] a combined analysis of the new Daya Bay results [29] with previous measurements of the reactor neutrino rates has been performed, concluding that the sterile neutrino hypothesis gives a fit of comparable quality to the combined rate data as the 235 U-only hypothesis. Below we will present an analysis including previous rate measurements as well as recent energy-spectral data, reaching a similar conclusion as the authors of ref. [32].
The outline of the paper is as follows: in section 2, we repeat the statistical analysis of the Daya Bay data from ref. [29], comment on its interpretation, and carry out additional statistical tests. In section 3, we combine the Daya Bay measurement of the individual isotopic fluxes with the other globally available reactor data, paying special attention to the impact of the new NEOS and DANSS results. In section 4, we put our results in a wider context by also including ν e disappearance data from gallium radiochemical experiments, solar neutrinos, and accelerator experiments. We summarize and conclude in section 5. Technical details of our simulations are given in the appendix.
In this paper we restrict the analysis to the ( -) ν e disappearance sector, motivated by the reactor and gallium anomalies. The implications of these results for the LSNDν µ →ν e signal [33] in the context of global data on various oscillation channels will be presented in a forth-coming publication.
2 Daya Bay measurements of 235 U and 239 Pu fluxes In ref. [29], the Daya Bay collaboration has for the first time presented independent measurements of theν e fluxes from 235 U and 239 Pu fission. This analysis was enabled by the precise knowledge of the isotopic composition of the reactor fuel and its evolution with time, combined with the large statistics of the Daya Bay near detectors.
The fuel composition is parameterized by the fractional contribution F 239 of 239 Pu fissions to the total fission rate. There is an approximate 1-to-1 correspondence between F 239 and the fractional contributions of the other isotopes, F 235 for 235 U, F 238 for 238 U, and F 241 for 241 Pu, see figure 2 of [29]. 8 bins in F 239 are used. Data are reported as effective inverse beta decay (IBD) yields σ, given in units of cm 2 per fission. We write the predicted IBD rates in each F 239 bin as Here, the index i runs over the four fissible isotopes; σ HM i is the IBD rate according to the Huber & Mueller flux predictions [1,2]; F a i gives the effective contribution of isotope i to the total fission rate in the a-th F 239 bin (a = 1 . . . 8); ξ i are four pull parameters which allow each flux to deviate from the predictions; and P i osc is the averaged oscillation probability at the Daya Bay Experimental Halls 1 and 2 (EH1 and EH2). (Data from EH3 is not used in this analysis.) The predictions from ref. [2] are used for the isotopes 235 U, 239 Pu, 241 Pu, and those from ref. [1] are used for 238 U. P i osc depends on the oscillation parameters and has a small dependence on the isotope i due to the slightly different neutrino spectra for each isotope. In the region ∆m 2 41 0.05 eV 2 , oscillations are averaged out completely, P i osc becomes independent of the isotope i and acts just as a global normalization factor, P i osc ≈ 1 − 1 2 sin 2 2θ 14 . For smaller values of ∆m 2 41 we take into account the correct oscillatory behaviour.
As a test statistic for the analysis we use the χ 2 function Here, σ a obs and σ a pred are the observed and predicted IBD yields in the a-th F 239 bin. The covariance matrix V ab includes statistical and correlated systematic errors. The covariance matrix as well as σ a obs and F a i are taken from the supplementary material of ref. [29]. The term χ 2 flux (ξ i ) constrains the nuisance parameters ξ i , and depending on the analysis we adopt different assumptions for it. When we impose the Huber & Mueller flux predictions ("fixed fluxes"), this term takes into account the systematic uncertainties on the fluxes as published in refs. [1,2]. We will also perform a "free fluxes" analysis, where ξ 235 and ξ 239 are allowed to vary freely. In this analysis, χ 2 flux (ξ i ) still imposes a weak 1σ constraint of Analysis χ 2 min /dof gof sin 2 2θ bfp 14 ∆χ 2 (no osc) fixed fluxes + ν s 9.8/(8 − 1) 18% 0.11 3.9 free fluxes (no ν s ) 3.6/(8 − 2) 73% Table 1: Fits to the Daya Bay flux measurements. The "fixed fluxes" analysis assumes the Huber & Mueller flux predictions [1,2] (accounting for their quoted uncertainties) and includesν e disappearance into sterile neutrinos ν s . We assume ∆m 2 41 0.05 eV 2 , so that oscillations are in the averaging regime. For the "free fluxes" analysis, fluxes are allowed to vary freely, but θ 14 is set to zero. The goodness-of-fit (gof) p-values are calculated by Monte Carlo simulation and agree roughly with the χ 2 approximation.
10% relative to the Huber & Mueller predictions on the subleading isotopes (ξ 238 and ξ 241 ) to avoid unphysical results.
In table 1 we show the results of our fit to the Daya Bay flux data under the "fixed fluxes" and "free fluxes" assumptions. We assume ∆m 2 41 0.05 eV 2 so that the predictions become independent of ∆m 2 41 and the only relevant oscillation parameter is θ 14 . For an analysis including sterile neutrinos, the number of degrees of freedom is thus 8 − 1 = 7. For the "free fluxes" analysis assuming no sterile neutrino (θ 14 = 0), the number of degrees of freedom is 6, accounting for the two unconstrained pull parameters ξ 235 and ξ 239 . We have checked by explicit Monte Carlo simulation, that χ 2 min follows indeed a χ 2 -distribution with 7 and 6 dof, respectively, to very good accuracy. As is clear from table 1, the hypothesis of free fluxes gives a better fit to the data, with a goodness-of-fit (gof) p-value of 73%. However, the sterile neutrino hypothesis also has an acceptable gof with a p-value of 18% and therefore cannot be rejected at reasonable confidence. The best fit point has sin 2 2θ 14 = 0.11, and for fixed fluxes the no oscillation hypothesis is disfavoured with ∆χ 2 = 3.9, corresponding to about 2σ (1 dof) exclusion.
To quantify Daya Bay's preference for free fluxes compared to oscillations into sterile neutrinos, we construct a test statistic which compares the two hypotheses. Here, we call H 0 the hypothesis that the Huber & Mueller fluxes are correct, butν e can disappear due to oscillations at the eV 2 scale; H 1 is the hypothesis that the predicted normalization of the four fluxes is not trustworthy and should be left free. Hence, H 0 corresponds to the "fixed fluxes" analysis including oscillations, and H 1 to the "free fluxes" analysis without oscillations. Note that to a good approximation H 0 is a subset of H 1 , since oscillations basically act as a global normalization. Under this assumption we expect that T follows a χ 2 distribution with 1 dof. We have verified by Monte Carlo simulation that this is indeed approximately true and holds independently of the true value of θ 14 , as long as sin 2 2θ 14 0.6. We find where the p-value is evaluated by Monte Carlo simulation (we obtain 1.2% (2.5σ) in the approximation of a χ 2 distribution with 1 dof). Hence, the sterile neutrino hypothesis is disfavoured compared to the "free fluxes" hypothesis at the 99.3% confidence level. This is in qualitative agreement with the results of [29]. The reason for the slightly lower value for T obs in eq. (2.4) compared to the value of 7.9 obtained in [29] is that our "fixed fluxes" analysis includes the uncertainties in the Huber & Mueller flux prediction for the four isotopes (encoded in ξ i factors in eq. (2.1)), whereas ref. [29] does not. In summary, while the Daya Bay flux measurements favour the "free fluxes" hypothesis over the sterile neutrino hypothesis, the latter still provides a good fit to the data (gof of 18%). Assuming the Huber & Mueller flux predictions to be correct,ν e disappearance is favored over the no oscillation hypothesis at 2σ. Therefore, we proceed with the sterile neutrino analysis and combine the Daya Bay flux data with all other reactor data.

Data Sets Used and Analysis Procedure
Our analysis of reactor neutrino data is based on ref. [16], where technical details and further references can be found. Here we describe the main differences and updates with respect to that analysis. Table 2 summarizes the data sets included in our global fit. We distinguish experiments comparing the predicted and measured total neutrino fluxes and experiments that use spectral information in the analysis.
Compared to ref. [16], we have added a 4th Krasnoyarsk data point [34], see ref. [30] for details. For RENO and DoubleChooz, we include in our analysis the total rate measurements at the near detectors given in refs. [35] and [30], respectively. For RENO, we also include the ratio of total rates at the near and far sites from ref. [36]. We do not include the RENO and Double Chooz measurements of the neutrino spectrum, as a reliable interpretation of these measurements in the context of sterile neutrino models turned out to be difficult, and as their statistical power is far inferior to that of Daya Bay. The Daya Bay measurement of the isotope-dependent neutrino fluxes [29] discussed in section 2 are included as constraints on the flux normalizations in a consistent way, correlated between all reactor data. In other words, the nuisance parameters ξ i in eq. (2.1) are the same for each experiment, and the pull term χ 2 flux (ξ i ) in eq. (2.2) is added only once to the gloabl χ 2 . Oscillation effects in Daya Bay that affect the extraction of the fluxes are of course taken into account.
New data on the reactor anti-neutrino spectrum are included from the Daya Bay, NEOS, and DANSS experiments. In ref. [25] the Daya Bay collaboration has presented constraints on sterile neutrino mixing by fitting two ratios built out of the spectra recorded at the three experimental halls. We follow this strategy but use the larger data sample from ref. [46]. We fit the ratio of the binned spectra at EH3/EH1 and EH2/EH1. Details of the analysis are given in appendix A.1. The NEOS collaboration [26] has reported a high statistics measurement of the anti-neutrino spectrum at a distance of 24 m from the core of a 2.8 GW nuclear reactor. We include their results using the ratio of the measured spectrum to the shape predicted from the flux measured at the Daya Bay EH1 and EH2 detectors [21]. Details of the analysis are given in appendix A.2. The DANSS collaboration has reported preliminary results on the anti-neutrino event spectrum at distances of 10  rate measurements, data on the neutrino energy spectrum, and the very-long baseline experiment Kam-LAND. The column "# Data" gives the number of data points entering the corresponding χ 2 function. The total number of data points is 239. The acronym "EH" stands for "experimental hall" in Daya Bay, with EH1, EH2 being the two near detectors halls and EH3 the far detector hall. The last column highlights the most recent data sets (since summer 2016). In the text, we refer to these data sets as "new", to the previous ones as "old".
from a reactor core [28]. We include these measurements by fitting the bin-by-bin ratio of the two spectra, see appendix A.3 for details. In all cases we have verified that we can reproduce to good accuracy the results of the respective experimental collaborations, when data are analysed under the same assumptions.
As in section 2, we will in the following present two different global fits: one with fixed fluxes, in which we take the predicted anti-neutrino fluxes and their uncertainties at face value; and one with free fluxes, in which the flux from each fissible isotope is allowed to float independently. In the case of fixed fluxes, the predictions from ref. [2] are used for the isotopes 235 U, 239 Pu, 241 Pu, and those from ref. [1] are used for 238 U. In this analysis we always take into account the quoted systematic uncertainties on the fluxes [2], including correlations between isotopes and energy bins. These uncertainties are of order few %. In the fit with free fluxes, the normalizations of the 235 U and 239 Pu fluxes are left completely unconstrained, whereas for the subleading fluxes from 238 U and 241 Pu we impose a weak constraint of ±20% (1σ) in order to avoid unphysical values. (Note that this is more conservative than the ±10% uncertainty we used in section 2 to match Daya Bay's analysis.) Thanks to the Day Bay flux measurement [29] as well as the slightly different isotopic compositions of the different reactor cores at which experiments have been conducted, the data itself provides sufficient information on the flux normalizations, see e.g., refs. [30][31][32].
The predicted reactor neutrino spectra are used neither in the "fixed fluxes" nor in the "free fluxes" fit. Instead, when using spectral information, we always compare measured spectra at different baselines. Ignoring the predicted shape of the neutrino spectrum is motivated by the unexplained bump at E ν ∼ 5 MeV, which indicates that our theoretical understanding of the spectra is incomplete 1 . In Daya Bay, the spectral comparison is between the different experimental halls; in DANSS, it is between measurements at two different locations using the same (movable) detector; for NEOS, the observed spectrum is compared to the measured spectrum from Daya Bay; for Bugey-3, we have modifed our previous analysis [16,48] by introducing a free pull parameter for each energy bin and correlating it between the three detector distances of 15, 40, and 95 m (this leads to 60 − 25 = 35 degrees of freedom for Bugey-3). 2 Let us remark that our analysis of spectral data is somewhat overconservative because we allow the spectra to be deformed in an uncorrelated way between different experiments, i.e., we do not correlate the energy bins between different experiments (except for the NEOS-Daya Bay comparison).
Concerning the oscillation physics, we use the full 4-flavour disappearance probability. For our parameterization conventions see ref. [16]. The parameters ∆m 2 21 , θ 12 , ∆m 2 31 are fixed to the 3-flavour best fit points, while θ 13 is left free (since the used data determine it with good accuracy).

Results for the Combined Analysis of Reactor Data
In fig. 1 we illustrate the impact of the recent oscillation analyses from NEOS, DANSS, and the Daya Bay spectrum. In khaki, we show the 2σ allowed parameter region in the sin 2 2θ 14 vs. ∆m 2 41 plane, based on data predating the summer conferences 2016. The corresponding data sets are marked with "-" in the last column of table 2. The black and green dashed contours show the new exclusion limits from Daya Bay and DANSS, and the blue contours depict the limit from the combined NEOS and Daya Bay spectral analysis. Due to the relatively long baseline of the Daya Bay detectors, these data constrain the region of ∆m 2 41 0.3 eV 2 , while both NEOS and DANSS are most sensitive in the RAA region around few eV 2 .
As mentioned above, the NEOS analysis is based on the ratio of the spectra in the NEOS detector to the one extracted from Daya Bay EH1 and EH2 data. When taking into account the ∆m 2 41 dependence of the oscillations in the Daya Bay near detectors, NEOS data lead to a closed regions with a best fit point below ∆m 2 41 0.1 eV 2 , which is, however, excluded by the Daya Bay spectral data at the far detector (EH3). Therefore, we decided to show only the combined NEOS+Daya Bay (spectrum) constraint, in order to avoid the effect of the minimum in the excluded region. The complementarity of the two data sets is clearly visible, by comparing the blue and black curves. Figure 1: Allowed parameter regions at 2σ (2 dof) for the "flux-fixed" analysis, for the "old" data sample defined in table 2 (khaki regions), for the DANSS [28], Daya Bay spect. [46], the combined Day Bay spect. + NEOS [26] oscillation analyses, and the combined region of all data including also Daya Bay flux [29] (red regions). The cross marks the best fit of the combined region. From the energy spectral data shown in fig. 2 it is clear that both, NEOS and DANSS data prefer oscillations for flux-fixed compared to a constant re-scaling of fluxes. The effect is more pronounced for the reactor-only best fit (blue-solid curve) but still visible in the global disappearance analysis to be discussed in the next section (green-dashed curves).
Remarkably the wiggles in the exclusion curves in fig. 1 from NEOS and DANSS partially match onto each other, leaving large parts of the RAA region unconstrained. The red shaded regions in fig. 1 include all data sets from table 2. The fact that this region is pushed to larger mixing angles compared to the "old" region is due to the Daya Bay flux data, which prefer somewhat larger values of the mixing angle. Figure 3 shows the allowed regions at 1, 2, 3σ confidence level (2 dof) for the combined analysis of all reactor neutrino experiments listed in table 1 and compares results for the analyses with fixed and free fluxes. Note that "free fluxes" now includes also oscillations in Details of the observables and predictions for the two experiments can be found in the appendix. The large distortions of the NEOS data below 2 MeV are within the systematic error band (not shown, see [26]).
addition to leaving fluxes free in the fit, whereas before, it meant just rescaling of fluxes, but θ 14 = 0. In table 3, we summarize the best fit points, the corresponding χ 2 /dof values, and the ∆χ 2 between the best fit point and the no oscillation hypothesis. We observe that the significance of the RAA slightly increases from a p-value for no-oscillations of 0.91% (2.6σ) for "old" data to 0.36% (2.9σ) for combined reactor data. Clearly for the flux-free analysis the significance for oscillations decreases, but for the combined reactor data a hint for oscillations remains even for flux-free (p-value of 6.1%, 1.9σ), mostly driven by NEOS, cf. eq. (3.1). Note that in fig. 3 the preferred regions from the flux-fixed analysis are consistent with the flux-free exclusion limits. In table 4 we provide the values of the pull parameters ξ i , cf. eq. (2.1), obtained in the flux free analysis at the oscillation best fit point and for no oscillation. In the latter case, the relative rescaling of the two main flux contributors, ξ 235 and ξ 239 , qualitatively agree with the results in refs. [29,32] 3 , with 235 U being the main contributor to the flux deficit. At the oscillation best fit point, the suppression for the 235 U and 239 P fluxes due to the pull parameters decrease because of the suppression from oscillations. But as in the no oscillation case, the larger suppression corresponds to 235 U . Consistently, in the same way as the flux pulls decrease when including oscillations, the mixing angle at the best fit decrease when going from the flux fixed to the flux free analysis, cf. table 3.
In fig. 4 (left panel) we show the marginalized ∆χ 2 as a function of ∆m 2 41 for both the fixed fluxes and free fluxes analyses of the combined reactor data. We observe the two prominent minima around 1.7 and 3 eV 2 , both alowed at 1σ for fixed fluxes. For free fluxes,   we find the best fit at 1.7 eV 2 , with several other local minima below the 2σ threshold. Note that the maximal values of these curves correspond to the ∆χ 2 for no oscillations as given in table 3. The reason is that the no oscillation case can be obtained for any ∆m 2 41 when minimizing χ 2 with respect to the mixing angle.
We can also perform the same test as in section 2, comparing the two hypotheses fluxfixed + sterile neutrino versus flux-free without sterile neutrino. From the numbers in table 3 we obtain for the test statistic T defined in eq. (2.3): T obs = 2.9 (all reactors) .

(3.2)
The spectral distortions observed in NEOS and DANSS prefer sterile neutrino oscillations over flux rescaling and therefore the preference for the flux-free fit obtained by Daya Bay flux data decreases. We conclude that in light of the global data we cannot reject the sterile neutrino hypothesis compared to the flux-free hypothesis. Let us remark that due to spectral data, the sterile neutrino hypothesis is now no longer a subset of the flux-free hypothesis (without oscillations). Therefore the interpretation of above values for T obs in terms of pvalues is not straight forward and we limit our conclusion to qualitative statements relative to the result obtained for Daya Bay flux data alone in eq. (2.4). Note that when we use the same 10% prior uncertainty on the subleading fluxes 238 U and 241 Pu as in section 2 (instead of 20% adopted in eq. (3.2)), the value for T quoted in eq. (3.2) decreases further to 2.1, to be compared to 6.3 obtained for Daya Bay flux alone, see eq. (2.4). Interestingly, for the old reactor data, oscillations are even preferred over flux-free: For this data set the best fit for oscillations is obtained at a rather low value for the masssquared difference around 0.4 eV 2 . For this value, the observed rates at different baselines can be fitted better than in the case of constant flux reduction, which leads to a preferrence of oscillations. This result is in qualitative agreement with ref. [32], where further discussions of this effect can be found. We note that in the global analysis of all recent data, such low values of ∆m 2 41 are disfavoured at more than 3σ by spectral data, most importantly from Daya Bay, cf. fig. 4, and therefore the flux-free hypothesis gives still a slightly better fit than oscillations + the Huber & Mueller predictions.
Finally, let us note that DANSS result have not been published yet, and our analysis is based on preliminary results presented in a conference talk [28]. We comment briefly on the impact of those data on our result. Removing DANSS from the combined analysis decreases the observed value for T from 2.9 to 1.4, i.e., the sterile neutrino and the flux-free hypotheses become statistically equivalent. Furthermore, the ∆χ 2 of no-oscillations compared to the oscillation best fit point increases by about 1 unit when removing DANSS, both for the flux-fixed and flux-free analyses. While DANSS does show a weak preference for oscillations, there is a slight tension between the best fit points preferred with (∆m 2 41 ≈ 3 eV 2 ) and without (∆m 2 41 ≈ 1.7 eV 2 ) using DANSS, c.f. fig. 1. This leads to a slightly larger preference in favour of oscillations when DANSS is left out. However, quantitatively the impact is small and qualitatively the picture remains robust, irrespective of using preliminary DANSS data or not.
4 Global analysis of ( -) ν e disappearance data

Non-reactor data
In addition to the reactor neutrino data discussed before, there are other experiments which are sensitive to ( -) ν e disappearance and can therefore provide complementary information. In this work, we consider in particular the data listed in table 5: • Solar neutrino data. We update our previous analysis [16] by accounting for the 2055day energy and day/night asymmetry spectrum from Super-Kamiokande phase 4 [55]. We also include the new measurement of neutrinos from the proton-proton (pp) fusion chain in the Sun recently presented by Borexino [61]. In addition, we make use of the total rates from the radiochemical experiments Chlorine [49], GALLEX/GNO [50] and SAGE [51], the electron scattering data binned in energy and zenith angle from all the previous Super-Kamiokande phases [52][53][54], the individual data sets from the three phases of SNO [56][57][58], and the Borexino phase-I data samples consisting of 740.7 days of low-energy data [59] and 246 live days of high-energy data [60]. Thus the solar neutrino data used in our analysis consists of 325 data points. Details of the implementation of the oscillation probabilities and relevant parameters can be found in Appendix C  Table 5: Experimental data which we combine with the reactor data from table 2 in our global ν e /ν e disappearance analysis. In the last column we indicate updates with respect to ref. [16]. The total number of data points of non-reactor data is 361.
• Radioactive source experiments. The calibration of Gallium solar neutrino experiments has been tested by deploying radioactive sources in the GALLEX [50,62] and SAGE [63,64] detectors. Both experiments have updergone two calibration campaigns: one with 37 Ar and one with 51 Cr in the case of SAGE, and both with 51 Cr in the case of GALLEX. All four measurements have reported an event rate about 10% to 20% lower than expected, a fact commonly known as the "Gallium anomaly". A re-evaluation [69] of the poorly-known contribution of 71 Ge excited states to the relevant 71 Ga → 71 Ge nuclear cross-section presented in [70] has not settled the issue. The deficit may be explained by the hypothesis of ν e disappearance due to oscillations with a mass-squared difference at the eV 2 scale, and is therefore a relevant ingredent of our study. A detailed discussion of our implementation is provided in sec. 3.2 of [16].
• ν e scattering on 12 C. The LSND [68] and KARMEN [65,66] experiments have measured the reaction ν e + 12 C → e − + 12 N, where the 12 N decays back to 12 C + e + + ν e with a lifetime of 15.9 ms. The experimental signature for this process, characterized by the observation of a prompt electron followed by a delayed positron, allows for precise identification of signal events and efficient rejection of backgrounds. Both electron and positron energies are recorded, thus allowing to reconstruct the parent neutrino energy.
No deviation from the no-oscillation hypothesis is observed in either detector, which results in a limit on the sterile neutrino admixture to ν e [67]. Details of our implementation of LSND and KARMEN results on 12 C scattering are provided in Appendix E.1 of [16].
In our analysis, correlations among the various experimental results within each of the three classes of data listed above (solar, radioactive source, scattering on 12 C) are properly taken into account, whereas correlations between different classes are neglected. In principle, the GALLEX and SAGE experiments contribute both to the solar neutrino analysis and to the Gallium anomaly, thus introducing a correlation among these two sets. However, we have verified that the solar neutrino rate in GALLEX and SAGE is completely dominated by the ground-state cross-section, which has a small error. Conversely, the main source of uncertainty affecting the Gallium anomaly comes from the two excited levels GT175 and GT500 (see [16] for details), whose contribution to the solar neutrino interaction rate is only at the percent level. Therefore, a proper treatment of the correlations between the Gallium anomaly and solar neutrino data, despite introducing a non-trival complication, would add very little to the results of our study.

Results
The results of our global analysis of all ( -) ν e disappearance experiments are shown in fig. 5 for the "fixed fluxes" and "free fluxes" analyses. ∆χ 2 profiles as a function of ∆m 2 41 are shown in the right panel of fig. 4. Best-fit points and χ 2 values are reported in the last two rows of table 3. We observe that, both for free fluxes and for fixed fluxes, the combined fit is largely dominated by reactor neutrino data. The total number of data points in this analysis is 600, and the oscillation fit includes the six parameters ∆m 2 41 and the mixing angles θ 12 , θ 13 , θ 14 , θ 24 , θ 34 ; the other mass-squared differences and θ 23 are fixed to their 3-flavour best fit points. Although we do take into account the two complex phases on which solar oscillation probabilities formally depend, their impact on the χ 2 is negligible and we do not count them as degree of freedom in the fit, see appendices of ref. [16] for a discussion of complex phases.
For what concerns solar neutrino data, the mass-squared difference ∆m 2 41 implied by the reactor anomaly is virtually infinite in the calculation of the P ee survival probability, hence its specific value is not constrained by solar experiments. The bound on θ 14 is mainly driven by the good agreement between the theoretical expectation of the 8 B neutrino flux, which is predicted by the Standard Solar Model, and its precise determination in high-energy solar experiments. This includes direct measurements (through neutral current interactions in SNO) and indirect measurements (through the combination of charged current and elastic scattering data in SNO and SK). The inclusion of a sterile neutrino admixture with ν e implies an overall reduction of the flux of active neutrinos at the detector, thus spoiling such agreement. This results in an upper bound on θ 14 , which in the case of the "fixed flux" analysis is fully compatible with the entire region allowed by reactor data, thus adding little to the global analysis. In the "free fluxes" case, solar data help restricting θ 14 at ∆m 2 41 4 eV 2 , where reactor experiments lose sensitivity because the oscillation length becomes very short, implying a uniform suppression of the reactor neutrino flux in all reactor experiments, but no spectral features. Such a uniform suppression cannot be disentangled from a rescaling of the flux normalization. Similar arguments also apply to the LSND & KARMEN data on 12 C, which show no deviation from the standard oscillation scenario and therefore impose an upper bound on θ 14 in the large ∆m 2 41 region. The situation regarding the Gallium anomaly is somewhat different. As already explained, the GALLEX and SAGE experiments observe a deficit which can be interpreted in terms of sterile neutrino oscillations. However, its 2σ allowed region shows little overlap with the reactor region, except for a small area at large ∆m 2 41 . In general, Gallium data favor a larger value of θ 14 than reactor data. It should be noted, however, that while the lower bound on θ 14 from GALLEX and SAGE is rather weak, with the no-oscillation value θ 14 = 0 disfavoured only by ∆χ 2 = 8.72 with respect to the best-fit point (see sec. 3.2 of Ref. [16]), the upper bound on θ 14 from reactor data is pretty strong for ∆m 2 41 5 eV 2 . Therefore, a combination of reactor and gallium data naturally favors the reactor region, rather than the GALLEX and SAGE one, so that the net contribution of gallium data is vastly reduced.
Indeed, as can be seen also in table 3 and fig. 4, the results of the global analysis differ little from those of the reactor-only one, with a very similar ∆χ 2 for no oscillations in the "fixed fluxes" analysis of 11.3 versus 11.9 (p-values of 0.36% versus 0.26%). For free fluxes, the impact of Gallium data is somewhat larger, increasing the ∆χ 2 of no oscillations from 5.6 for reactors to 7.0 for the global data (p-values of 6.1% versus 3.1%). This corresponds to a hint in favour of oscillations at 2.2σ, resulting in closed regions at the 2σ level for the flux free analysis, visible in fig. 5 (right panel).
The test statistic T defined in eq. (2.3) for discriminating between the flux-fixed + oscillations versus flux-free + no oscillations decreases from 2.9 (reactor-only data) to 2.1 for the global disappearance data.

Discussion and Conclusions
In this paper we have investigated the status of the sterile neutrino hypothesis in the context of the global data on ( -) ν e disappearance, in the light of new results from reactor neutrino experiments. In particular, we have considered the impact of first results from the NEOS and DANSS short-baseline reactor experiments, as well as the recent determination of the inverse-beta decay rate induced by neutrinos from different fission isotopes by Daya Bay. In our reactor data analysis we have taken into account the disagreement of data and predictions in the spectral shape ("5 MeV bump") by using only relative spectra at different baselines.
We confirm that Daya Bay flux measurements alone favour the hypothesis that the source of the reactor anomaly is the flux of 235 U over the hypothesis of sterile neutrino oscillations. However, the sterile neutrino hypothesis also provides a good description of the data (p-value of 18%) and hence cannot be excluded. Therefore we combine this Daya Bay result with the remaining data from reactor experiments assuming the presence of a sterile neutrino. For the global reactor data, actually, the preference for re-scaling the 235 U flux over oscillations is reduced compared to the Daya Bay flux data alone, see eqs. (2.4) and (3.2), and the sterile neutrino hypothesis + Huber & Mueller flux predictions cannot be rejected. The main reason for this are features in the energy spectra reported by the DANSS and NEOS experiments, which prefer oscillations over a rescaling of fluxes.
Since the sterile neutrino hypothesis cannot be excluded, we present updated determinations of oscillation parameters. The shape of the exclusion curves from DANSS and NEOS spectral data leaves allowed parameter space consistent with each other and with the remaining reactor neutrino data. The combined analysis leads to islands in the allowed parameter space around (∆m 2 41 , sin 2 2θ 14 ) ∼ (1.7 eV 2 , 0.04) and (3 eV 2 , 0.1), and the significance of the sterile neutrino explanation of the reactor anomaly remains slightly below 3σ (the no-oscillation hypothesis has a p-value relative to the best fit point of 0.36%). We have also preformed an oscillation fit leaving the neutrino fluxes from the four main fission isotopes completely free. Although the significance of the anomaly decreases, a hint for oscillations remains at the 1.9σ level and the exclusion curves on the oscillation parameters are consistent with the best fit regions obtained in the analysis with fixed fluxes.
Finally, we have provided updated results from a global fit to ( -) ν e disappearance data, including the Gallium anomaly and constraints from solar neutrinos and from ν e -12 C scattering in LSND and KARMEN. The results of the global analysis are largely consistent with the reactor-only fit, and the indications for sterile neutrino oscillations remain at a significance close to 3σ (2σ) with respect to no oscillations in the case of flux-fixed (flux-free) analyses.
In conclusion, present data on ( -) ν e disappearance is still consistent with sterile neutrino oscillations at the eV scale with modest significance. To definitely clarify the question raised in the title, more data is needed, which can be expected in the near future from new shortbaseline reactor experiments as well as radioactive source experiments, see e.g., ref. [17] for references and a review of sensitivities of upcoming experiments. A Technical details on our analyses A.1 Daya Bay sterile neutrino fit Using the data from the different Daya Bay experimental halls and in the different energy bins, a χ 2 depending on θ 13 , θ 14 , ∆m 2 41 and pull parameters is computed as follows (∆m 2 31 , ∆m 2 21 , and θ 12 are included but kept fix at their 3-flavour best fit values): Here, O H i are the observed number of events in experimental hall H and energy bin i, and B H i are the corresponding predicted background. The measured event rates and background predictions are taken from the supplementary material of ref. [46]. N H i are the predicted event numbers in experimental hall H and bin i, see below. σ stat,HH i are the statistical errors of the ratios Finally, p is the vector of nuisance parameters accounting for the systematic uncertainties, and V p is the corresponding covariance matrix. It includes the uncertainties in the relative energy scale and the detection efficiency as well as their correlation, which can be found in table VIII of ref. [46]. Since we are using binby-bin ratios between detectors at different baselines, errors in the flux predictions and in the inverse beta decay cross sections will mostly cancel. The minimization over the pull parameters is done by linearizing the dependence of N H i on p and then solving a linear system of equations.
Since each experimental hall in Daya Bay houses several detectors, N H i is obtained by summing the contributions from all detectors in hall H. The predicted number of events in an individual detector d and energy bin i is where • the indices i, r, d and iso refer to energy bins, reactors, detectors, and fissible isotopes, respectively.
• d is the detector efficiency, taken from table VI in ref. [46]. We consider the efficiencies µ and m (corresponding to loss of events from the muon veto and multiplicity veto, respectively) as well as the variation in the number of target protons in each detector, ∆N p .
• L rd is the baseline between reactor r and detector d.
• E ν and E rec are the true neutrino energy and the energy reconstructed by the detector, respectively. The detector response function R(E rec , E ν ), taken from the supplementary material of ref. [46], describes the mapping between E ν and E rec .
• φ iso (E ν ) are the flux predictions from refs. [1,2], and f iso are the fission fractions. For each isotope, f iso is computed as the average of the fission fractions in table 9 of ref. [21].
• P rd νe→νe (E ν ) is the oscillation probability. • A H is a normalization factor, which is fixed by requiring the total predicted number of events in hall H without oscillations to match the corresponding number given in the supplementary material of ref. [46].

A.2 NEOS
Our fit to NEOS data is based on fig. 3(c) of ref. [26], where the data are presented as ratios between observed event rates in NEOS and a prediction based on the unfolded Daya Bay anti-neutrino spectrum from ref. [21]. We adopt the following χ 2 function: Here, O N i is the NEOS data point in energy bin i, and P N i is the theoretical prediction. To obtain the latter, we have to account for the fact that the unfolded Daya Bay spectrum is based on the assumption of three-flavour oscillations. We therefore have to unfold threeflavour oscillations (which are, however, small in Daya Bay and negligible in NEOS) and fold in four-flavour oscillations: where P Exp nν are the predicted event numbers in bin i for experiment Exp = NEOS, DB in the nν neutrino framework. The latter are obtained based on the Huber-Mueller fluxes [1,2]. The covariance matrix V ij in eq. (A.3) includes statistical errors extracted from fig. 3(c) in ref. [26] as well as the covariance matrix for the Daya Bay flux determination provided in ref. [21].
For the Daya Bay predictions we take into account an average of the near detectors (EH1 and EH2) as used for the Daya Bay unfolded spectrum in [21] and they are calculated as in appendix A.1. The number of events per bin in NEOS is computed in analogy to eq. (A.2). Since the NEOS detector is very close to the source, we also take into account the finite sizes of the reactor core and of the detector by averaging the oscillation probability weighted by 1/L 2 over L = (24 ± 1.5) m. Since no response function R(E rec , E ν ) is provided by the NEOS collaboration we adopt the model proposed in ref. [23] consisting of a a Gaussian for E rec > E p (E p = E ν − 0.8 MeV) and a rescaled Gaussian plus a constant value for E rec < E p to account for photons or positrons escaping the detector. In order to reproduce the NEOS spectrum from figure 3(b) in [26], we assume energy scale non-linearity effects based on the information on non-linearity provided by Daya Bay in the supplementary material of ref. [46].

A.3 DANSS
For the DANSS experiment, we use the preliminary data shown on slide 10 of ref. [28]. The data are given as ratios of observed event numbers between the two detector positions at L = 12.7 m (down) and L = 10.7 m (up) from the center of the reactor core. The data are divided into 30 energy bins of equal width, ranging from E p = 1.0 MeV to E p = 7.0 MeV.
Here, E p is the kinetic energy of the outgoing positron in inverse beta decayν e + p → n + e + . The χ 2 for DANSS is where the predicted down/up ratios P i are computed as ratios of oscillation probabilities, weighted by the geometric factor 1/L 2 . To account for the size and geometry of the detector and the reactor, we average the oscillation probabilities (divided by L 2 ) over L = L 0 ± 4.0 m. Here, L 0 is taken to be 12.85 m for the lower detector position and 10.9 m for the upper one. These numbers are slightly larger than the distances between the center of the reactor core and the center of the detector to account for the on average non-zero horizontal distance between the production and detection vertices. The energy resoluton of DANSS is modeled as a Gaussian with a width given by fig. 5 of ref. [72]. The covariance matrix V ij for DANSS includes only statistical uncertainties and a 2% systematic uncertainty on the down/up ratios.