Neutrino-4 anomaly: oscillations or fluctuations?

We present a deep study of the Neutrino-4 data aimed at finding the statistical significance of the large-mixing short-baseline neutrino oscillation signal claimed by the Neutrino-4 collaboration at more than $3\sigma$. We found that the results of the Neutrino-4 collaboration can be reproduced approximately only by neglecting the effects of the energy resolution of the detector. Including these effects, we found that the best fit is obtained for a mixing that is even larger, close to maximal, but the statistical significance of the short-baseline neutrino oscillation signal is only about $2.7\sigma$ if evaluated with the usual method based on Wilks' theorem. We show that the large Neutrino-4 mixing is in strong tension with the KATRIN, PROSPECT, STEREO, and solar $\nu_{e}$ bounds. Using a more reliable Monte Carlo simulation of a large set of Neutrino-4-like data, we found that the statistical significance of the Neutrino-4 short-baseline neutrino oscillation signal decreases to about $2.2\sigma$. We also show that it is not unlikely to find a best-fit point that has a large mixing, even maximal, in the absence of oscillations. Therefore, we conclude that the claimed Neutrino-4 indication in favor of short-baseline neutrino oscillations with very large mixing is rather doubtful.


I. INTRODUCTION
The possible existence of sterile neutrinos is one of the current hot research topics in the quest for new physics beyond the Standard Model. In particular, the possible existence of light sterile neutrinos at the eV mass scale has been motivated by the LSND [1], MiniBooNE [2,3], Gallium [4][5][6], and reactor [7] anomalies, that may be a signal of short baseline active-sterile neutrino oscillations (see the recent reviews in Refs. [8][9][10]). Recently, a strong indication in favor of such oscillations has been claimed by the Neutrino-4 collaboration [11][12][13]. However, this result is controversial and has been criticized in Ref. [14] (see also the answer in Ref. [15]) and in Ref. [16].
In this paper we present a detailed analysis of the Neutrino-4 data and we discuss in particular the effect of the energy resolution of the detector, that seems not to have been taken into account in the analysis of the Neutrino-4 collaboration (as independently noted in Ref. [16]). We also present the results of a Monte Carlo implementation of the statistical analysis of the Neutrino-4 data that is more reliable than the standard methods based on the χ 2 distribution predicted by Wilks' theorem [17], as discussed in Refs. [18][19][20][21][22][23].
In Section II we present our analysis method, in Section III we discuss the results that we obtained using the standard statistical method based on Wilks' theorem, in Section IV we present the results obtained with a Monte Carlo statistical analysis, and finally in Section V we summarize our results and draw our conclusions.

II. ANALYSIS METHOD
The Neutrino-4 experiment is operating since 2014 close to the SM-3 reactor in Dimitrovgrad (Russia). It is a 90 MW research reactor with a compact core (42 × 42 × 35 cm 3 ) using highly enriched 235 U fuel. The Neutrino-4 detector is made of liquid scintillator with 0.1% gadolinium concentration for a better detection of the neutron in theν e + p → n + e + detection process of the reactorν e 's. The detector is segmented in 50 sections with active size 22.5 × 22.5 × 70 cm 3 disposed in 10 rows. The detector can move horizontally, on the side of the reactor, allowing measurements of theν e at distances between 6.4 m and 11.9 m from the center of the reactor. The detector recorded about 300ν e events per day at the average distance of 8 m from the reactor core. The last published data [12,13] correspond to 720 days of reactor-on data taking and 417 days of reactor-off background measurements.
The data have been taken by the Neutrino-4 collaboration at n L = 24 distances between 6.4 m and 11.9 m from the center of the reactor, at intervals of 23.5 cm, and divided in n E = 9 bins of prompt energy from 1.5 to 6 MeV with uniform ∆E = 500 keV width. Let us remind that the neutrino energy E is related to the prompt energy E p by E = E p + m n − m p − m e E p + 0.78 MeV. (1) The 500 keV prompt energy bin width was motivated by the calibration of the detector with a 22 Na γ-source that showed an energy resolution of about 142 keV at the 511 keV 22 Na line and about 217 keV at the 1274 keV 22 Na line (see Figure 22 of Ref. [13]). It was also calibrated with the process n + p → d + γ that showed an energy resolution of about 276 keV at 2.2 MeV (see Figure 21 of Ref. [13]). The Neutrino-4 collaboration collected the data in the 216 ratios where i = 1, . . . , n E is the neutrino energy index, k = 1, . . . , n L is the reactor-detector distance index, L k is the average distance of the k th distance interval, and N exp ik is the rate of observedν e events in the i th energy bin and k th distance interval.
The results of short-baseline neutrino oscillation experiments as Neutrino-4 are typically analyzed in the framework of 3+1 active-sterile neutrino mixing. This is the simplest extension of standard three-neutrino mixing that can explain short-baseline neutrino oscillations with minimal perturbations to the standard three-neutrino mixing global fit of solar, atmospheric and long-baseline (accelerator and reactor) neutrino oscillation data [24][25][26]. In the 3+1 framework it is assumed that there is a non-standard massive neutrino ν 4 with mass m 4 1 eV such that m 1 , m 2 , m 3 m 4 , where m 1 , m 2 , m 3 are the masses of the three standard massive neutrinos ν 1 , ν 2 , ν 3 . In the flavor basis, besides the three standard active neutrinos ν e , ν µ , ν τ there is a sterile neutrino ν s , and the mixing is given by where U is the 4 × 4 unitary mixing matrix. In this 3+1 framework, the effective short-baseline (SBL) survival probability of electron neutrinos and antineutrinos is given by where ∆m 2 41 = m 2 4 − m 2 1 and sin 2 2ϑ ee = 4|U e4 | 2 (1 − |U e4 | 2 ), where ϑ ee is an effective mixing angle that coincides with ϑ 14 in the standard parameterization of the mixing matrix (see the recent reviews in Refs. [8][9][10]).
In the framework of 3+1 neutrino mixing, the expected rate of events in the i th energy bin and k th distance interval of the Neutrino-4 experiment is given by where N 0 i is the rate of expected events in the i th energy bin without short-baseline neutrino oscillations at a reference distance where L k = 1, and the oscillating terms are averaged over the appropriate neutrino energy range and the distance uncertainty of the i th energy bin and k th distance interval. Therefore, the theoretical prediction R the ik for the ratios in Eq. (2) is given by This expression shows the advantage of considering the ratios in Eq. (2), which must be confronted with the theoretical ratios R the ik that do not depend on the values of rates N 0 i , whose estimation has large uncertainties, mainly due to the uncertainties of the theoretical reactor neutrino fluxes (see the reviews in Refs. [27,28]).
The Neutrino-4 collaboration further processed the data by dividing the n tot = 216 R exp ik 's in groups of n g = 8 values that correspond to neighboring L/E intervals, arguing that the L/E dependence of the ratios R exp ik "allows the direct demonstration of the effect of oscillations". In this way, they obtained the n L/E = 27 averages Averaged oscillations with energy resolution smearing "500 keV" data Values of the "500 keV" experimental ratios R exp j (black points with error bars). The blue, red, and magenta histograms show, respectively, the values of the theoretical ratios R the j obtained for sin 2 2ϑee = 0.38 and ∆m 2 41 = 7.26 eV 2 without any averaging of the oscillating terms, with averaging over the energy and distance intervals without energy resolution, and with averaging over the energy and distance intervals with energy resolution. The green histogram corresponds to the best-fit values sin 2 2ϑee = 1 and ∆m 2 41 = 7.16 eV 2 obtained with averaging over the energy and distance intervals with energy resolution.
where g(j) is the set of indices i and k that corresponds to the j th group of neighboring L/E intervals. The Neutrino-4 collaboration published only the first n L/E = 19 values of the ratios R exp j in Figure 52 of Ref. [13] ("500 keV" data), corresponding to a range of L/E from about 1 to about 2.5 m/MeV, arguing that for L/E > 2.5m/MeV the oscillations are averaged-out by the energy resolution of the detector. We reproduced these data and their uncertainties in Figure 1.
The Neutrino-4 collaboration presented also the results of an average of the R exp j 's obtained with the 500 keV energy bin width discussed above and similar R exp j 's obtained with the 125 and 250 keV energy bin widths, arguing that the average of different sampling of the data suppresses fluctuations. We think that this procedure is rather ad-hoc and inappropriate, since the energy resolution of the detector is larger than about 200 keV. This procedure was criticized also in Ref. [16] on the basis of similar arguments. Therefore, in this paper we consider first the Neutrino-4 "500 keV" data and then we discuss how the results change by considering the "125-250-500 keV" averaged data.
The averages R exp j must be compared with the corresponding theoretical averages We calculated the averaged oscillation terms with where φν e (E) is the reactor neutrino flux, σν e p (E) is the detection cross section, E and   I. Results of our fits of the Neutrino-4 "500 keV" and "125-250-500 keV" data: minimum χ 2 (χ 2 min ), goodness of fit (GoF) for 17 degrees of freedom, best fit values of sin 2 2ϑee and ∆m 2 41 , and difference ∆χ 2 NO between the χ 2 of no oscillations and χ 2 min . The last four rows give the p-value and the number of σ's (σ-value) corresponding to ∆χ 2 NO obtained considering the χ 2 distribution (with two degrees of freedom corresponding to two fitted oscillation parameters sin 2 2ϑee and ∆m 2 41 ) and with a Monte Carlo estimation of the true distribution.
with energy resolution σ Ep given by The well-known dependence of the energy resolution on the square root of the energy is due to the Poisson statistics of photoelectrons. The coefficient is adapted to match the above-mentioned calibrations of the detector shown in Figures 21 and 22 of Ref. [13] and agrees with the independent estimation in Ref. [16].

III. WILKS' CONFIDENCE INTERVALS
In this Section we present the results that we obtained from the fits of the "500 keV" and "125-250-500 keV" data considering the standard χ 2 distribution for log-likelihood ratios, as predicted by Wilks' theorem [17]. We perform the fit of the Neutrino-4 data with the least-squares function where R exp j and ∆R exp j are the experimental data and their uncertainties. The "500 keV" data, that we consider first, are given in Figure 52 of Ref. [13] and reproduced in Figure 1.
The Neutrino-4 collaboration fitted the "500 keV" data with the least-squares function and obtained a 3.5σ evidence of neutrino oscillations with best-fit values sin 2 2ϑ ee = 0.38 and ∆m 2 41 = 7.26 eV 2 . The allowed regions in the (sin 2 2ϑ ee , ∆m 2 41 ) plane obtained by the Neutrino-4 collaboration from the "500 keV" data are shown in Figure 50 of Ref. [13]. One can see that there are five narrow-∆m 2 41 allowed regions at 3σ. We cannot reproduce the Neutrino-4 analysis of the "500 keV" data using the least-squares function (13), because the values of the R exp ik have not been published by the Neutrino-4 collaboration. However, we can reproduce approximately the Neutrino-4 analysis of the "500 keV" data using the least-squares function (12). In particular, it is clear that the best-fit values of the oscillation parameters found by the Neutrino-4 collaboration should fit well the experimental FIG. 2. Contours of the allowed regions in the (sin 2 2ϑee, ∆m 2 41 ) plane obtained with the analyses without and with energy resolution of the Neutrino-4 (a) "500 keV" and (b) "125-250-500 keV" data. The allowed regions are calculated with the standard ∆χ 2 method based on Wilks' theorem (∆χ 2 = χ 2 − χ 2 min = 2.3, 6.2, 11.8 for 1σ, 2σ, 3σ, respectively). The figures show also the 95% C.L. exclusion curves of KATRIN [29], PROSPECT [30], STEREO [31], and solar νe's that were calculated with the same statistical method.  (10) is almost twice of the 500 keV bin width.
As shown in the second column in Table I, the best fit of the "500 keV" data taking into account the energy resolution of the detector corresponds to maximal mixing (sin 2 2ϑ ee = 1). This value is much larger than the best fit value sin 2 2ϑ ee = 0.38 found by the Neutrino-4 collaboration. The corresponding theoretical R the j values are shown by the green histogram in Figure 1.
We can fit the "500 keV" data with the best-fit values of the oscillation parameters found by the Neutrino-4 collaboration only by neglecting the energy resolution of the detector, i.e. by considering Under this assumption, we obtained the theoretical R the j values represented by the red histogram in Figure 1, that gives an acceptable fit of the data. For comparison, in Figure 1 we have also drawn the blue histogram corresponding to unaveraged oscillations, that is obtained by replacing sin 2 ∆m 2 . In this case, in the calculation of R the j there is only the averaging due to the mean of the eight R the ij with different L k and E i that contribute to each bin of L/E.
From the red histogram in Figure 1 one can see that the oscillations averaged over the sizes of the energy bins and distance intervals without the energy resolution smearing are strongly suppressed at large values of L/E with respect to unaveraged oscillations. On the other hand, for small values of L/E the suppression is rather small and one can fit the data with the best-fit values of the oscillation parameters claimed by the Neutrino-4 collaboration.
However, the magenta histogram in Figure 1 shows that the energy smearing due to the energy resolution of the detector suppresses strongly the oscillations at all values of L/E and the data cannot be fitted with the relatively small best-fit value of the effective mixing angle claimed by the Neutrino-4 collaboration (sin 2 2ϑ ee = 0.38). In order to fit the data taking into account the strong suppression of the oscillations at all values of L/E due to the energy resolution of the detector, we need a value of the effective mixing angle that is close to our maximal mixing best-fit value (sin 2 2ϑ ee = 1), as shown by the green histogram in Figure 1.
The first two columns of Table I give the results of our fits of the "500 keV" data without and with the detector energy resolution. Figure 2(a) shows a comparison of the corresponding allowed regions in the (sin 2 2ϑ ee , ∆m 2 41 ) plane. We calculated these allowed regions with the standard ∆χ 2 method based on Wilks' theorem (∆χ 2 = χ 2 − χ 2 min = 2.3, 6.2, 11.8 for 1σ, 2σ, 3σ, respectively), assuming, as usual, that the least-squares function (12) is proportional to the log-likelihood.
From Figure 2(a), one can see that the allowed regions and best-fit values of the oscillation parameters that we obtained without taking into account the detector energy resolution reproduce approximately the results of the Neutrino-4 collaboration [13]. Unfortunately, these results are not correct, because the smearing due to the energy resolution of the detector has a strong effect, that moves the best fit to maximal mixing and enlarges dramatically the allowed regions. From Figure 2(a) one can see that the 3σ allowed regions with energy resolution are not bounded for small values of sin 2 2ϑ ee , in contrast to those obtained without energy resolution. Therefore, as one can see from Table I, taking into account the detector energy resolution lowers the significance of the indication in favor of neutrino oscillations obtained from the "500 keV" data from 3.2σ to 2.7σ.
These σ values have been obtained from the p-values given in Table I, that have been calculated considering for the difference of the value of χ 2 without neutrino oscillations and the value of χ 2 min in the case of neutrino oscillations assuming a χ 2 distribution with two degrees of freedom (corresponding to the two fitted oscillation parameters sin 2 2ϑ ee and ∆m 2 41 ), according to Wilks' theorem [17]. A more reliable Monte Carlo estimation of the p-values and corresponding σ-values is discussed in Section IV. Figure 2(a) shows also the 95% C.L. exclusion curves of KATRIN [29] (see also Ref. [32]), PROSPECT [30], STEREO [31], and solar ν e 's that are calculated with the same method of the Neutrino-4 allowed regions. One can see that all these exclusion curves are in tension with the Neutrino-4 allowed regions. In particular, both the PROSPECT and STEREO exclusion curves disfavor the Neutrino-4 2σ regions obtained taking into account the energy resolution of the detector. The KATRIN exclusion region disfavors the Neutrino-4 1σ region surrounding the best-fit point at maximal mixing.
We calculated a new solar neutrino bound on sin 2 2ϑ ee [33][34][35][36][37][38] considering only solar neutrino data, without the KamLAND constraint, that depends on the absolute reactor neutrino fluxes. The degeneracy of the effects of ϑ ee = ϑ 14 and ϑ 13 is resolved by using the model-independent constraint on ϑ 13 obtained independently form the absolute reactor neutrino fluxes by the Daya Bay [39], RENO [40], and Double Chooz [41] reactor neutrino experiments through the comparison of the data measured at near and far detectors. Therefore, the limit that we obtained (sin 2 2ϑ ee < 0.22 at 95% C.L.) is model-independent and robust. One can see from Figure 2(a) that the solar neutrino bound disfavors the large-mixing Neutrino-4 allowed regions. Indeed, such large values of 3+1 active-sterile neutrino mixing would not be a minimal perturbation of the standard three-neutrino mixing and would spoil the standard three-neutrino mixing global fit of solar, atmospheric and long-baseline (accelerator and reactor) neutrino oscillation data [24][25][26].
As already mentioned in Section II, the Neutrino-4 collaboration considered also the values of R exp j averaged over 125, 250 and 500 keV energy bin widths, arguing that the average of different samplings of the data suppresses the fluctuations. They presented the results obtained in this way as the most reliable results of the experiment. As we already remarked in Section II, we disagree with this method because the procedure is rather ad-hoc and it is not appropriate to consider energy bin widths that are much smaller than the energy resolution of the detector. However, we are obliged to consider also this procedure in order to check the corresponding results of the Neutrino-4 collaboration.
As one can see from the fourth column in Table I and from Figure 2(b), also using the "125-250-500 keV" data we obtained a large best-fit value of the effective mixing angle, that is smaller than the maximal best-fit value obtained from the "500 keV" data, but still very large, with maximal mixing allowed within 1σ. This is what should be expected if one considers energy bin widths that are smaller than the energy resolution of the detector: there is slightly less averaging of the oscillations due to the smaller energy bin widths, but the dominant averaging effect due to the energy resolution of the detector still requires the effective mixing angle to be very large in order to fit the data.
On the other hand, from the analysis of the "125-250-500 keV" data the Neutrino-4 collaboration obtained the best-fit values sin 2 2ϑ ee = 0.26 and ∆m 2 41 = 7.25 eV 2 and the allowed regions in the (sin 2 2ϑ ee , ∆m 2 41 ) plane were shown in Figure 56 of Ref. [13]. However, according to our analysis this relatively small value of the effective mixing angle cannot fit the data if the energy resolution of the detector is taken into account, as illustrated by the magenta histogram in Figure 3. Here we have a situation similar to what we have found for the "500 keV" data: the "125-250-500 keV" data can be fitted with sin 2 2ϑ ee 0.26 only by neglecting the energy resolution of the detector. Indeed, as shown in Table I and Figure 2(b), neglecting the energy resolution of the detector we obtained the best fit at sin 2 2ϑ ee = 0.27 and ∆m 2 41 = 8.8 eV 2 , and there is an allowed region at 1σ at sin 2 2ϑ ee = 0.2 and ∆m 2 41 7.2 eV 2 , corresponding to the best fit of the Neutrino-4 collaboration.
Table I and Figure 2(b) show also that the indication in favor of neutrino oscillations obtained with the "125-250-500 keV" data is less significant than that obtained with the "500 keV" data if the energy resolution of the detector is neglected (2.7σ instead of 3.2σ), but it is almost the same when the energy resolution of the detector is taken into account (2.8σ and 2.7σ).  Table I).
The similarity of the statistical significances in favor of neutrino oscillations that we obtained from the analyses of the "125-250-500 keV" data without and with the energy resolution of the detector (2.7σ and 2.8σ, respectively) is due to the similarity of the quality of the best fit in the two cases, that correspond to almost equal values of χ 2 min , as shown in Table I. From Figure 3 one can see that the fit of the data in the two cases is different, but in both cases it is not very good for some data points. Figure 2(b) shows that also considering the "125-250-500 keV" data the PROSPECT, STEREO, and solar exclusion curves disfavor the Neutrino-4 2σ regions obtained taking into account the energy resolution of the detector, and the KATRIN exclusion region disfavors the Neutrino-4 1σ region surrounding the best-fit point at maximal mixing.
In conclusion of this Section, it is useful to summarize our findings: 1. Both the "500 keV" and the "125-250-500 keV" data of the Neutrino-4 experiment can be fitted well with the values of the effective mixing angle claimed by the Neutrino-4 collaboration (sin 2 2ϑ ee = 0.38 for the "500 keV" and sin 2 2ϑ ee = 0.26 for the "125-250-500 keV") by neglecting the energy resolution of the detector.
2. If the energy resolution of the detector is taken into account, the best fit value of the effective mixing angle is maximal (sin 2 2ϑ ee = 1) for the "500 keV" data and close to maximal (sin 2 2ϑ ee = 0.93) for the "125-250-500 keV" data.
3. If the energy resolution of the detector is taken into account, the indication in favor of short-baseline oscillations of the Neutrino-4 data is about 2.7 − 2.8σ if one assumes the approximate validity of Wilks' theorem.
4. If the energy resolution of the detector is taken into account, there is a strong tension between the Neutrino-4 allowed regions in the (sin 2 2ϑ ee , ∆m 2 41 ) plane and the 95% C.L. exclusion curves of KATRIN [29], PROSPECT [30], STEREO [31], and solar ν e 's.
In Section IV we present a more reliable Monte Carlo estimation of the statistical significance of the indication in favor of short-baseline oscillations of the Neutrino-4 data and the corresponding estimation of the allowed regions in the (sin 2 2ϑ ee , ∆m 2 41 ) plane.  Figure 2 of the allowed regions calculated assuming a χ 2 distribution for ∆χ 2 = χ 2 − χ 2 min .

IV. MONTE CARLO CONFIDENCE INTERVALS
It is well known that the conditions for the validity of Wilks' theorem are not always satisfied in neutrino oscillation experiments [18][19][20][21][22][23]. In this case, a Monte Carlo statistical analysis gives a more reliable estimation of the confidence Probabilities to obtain the best-fit in four ranges of sin 2 2ϑee in the absence of neutrino oscillations with the four methods of data analysis discussed in the paper.
intervals for the oscillation parameters. In this Section we present the results of Monte Carlo statistical analyses of the "500 keV" and the "125-250-500 keV" data of the Neutrino-4 experiment without and with the effect of the energy resolution of the detector. For each point on a grid in the (sin 2 2ϑ ee , ∆m 2 41 ) plane we generated a large number of random data sets (of the order of 10 5 ) with the uncertainties of the Neutrino-4 data set. For each random data set we calculated the value of χ 2 corresponding to the generating values of sin 2 2ϑ ee and ∆m 2 41 and we found the minimum value of χ 2 in the (sin 2 2ϑ ee , ∆m 2 41 ) plane. Denoting these two quantities by χ 2 MC (sin 2 2ϑ ee , ∆m 2 41 ) and χ 2 MC,min (sin 2 2ϑ ee , ∆m 2 41 ), we obtained the distribution of ∆χ 2 MC (sin 2 2ϑ ee , ∆m 2 41 ) = χ 2 MC (sin 2 2ϑ ee , ∆m 2 41 ) − χ 2 MC,min (sin 2 2ϑ ee , ∆m 2 41 ), that allows us to determine if the value of ∆χ 2 (sin 2 2ϑ ee , ∆m 2 41 ) = χ 2 (sin 2 2ϑ ee , ∆m 2 41 ) − χ 2 min (sin 2 2ϑ ee , ∆m 2 41 ) obtained with the analysis of the actual Neutrino-4 data is included or not in a region with a fixed confidence level. Figure 4 shows the comparison of the contours of the Monte Carlo allowed regions with the Wilks ∆χ 2 contours that were obtained in Figure 2 assuming a χ 2 distribution for ∆χ 2 = χ 2 − χ 2 min . One can see that in all the four analyses the Monte Carlo allowed regions are much larger than the Wilks ∆χ 2 allowed regions. In particular, for a fixed value of the confidence level smaller values of the effective mixing are allowed by the Monte Carlo calculation and the absence of neutrino oscillation is compatible with the data at 3σ in all the four analyses. This result is reflected in the decrease of the statistical significance of the Neutrino-4 indication in favor of neutrino oscillations shown in Table I. Considering only the more accurate analyses that take into account the energy resolution of the detector, one can see that the Monte Carlo calculation lowers the statistical significance of the indication in favor of neutrino oscillations of the "500 keV" data from 2.7σ to 2.2σ and that of the "125-250-500 keV" data from 2.8σ to 2.2σ.
Therefore, although our analyses of the Neutrino-4 data show that taking into account the energy resolution of the detector leads to a best-fit value of the mixing close to maximal, that is in contradiction with the bounds of other experiments, the allowed range of the mixing is large and includes the absence of mixing, and hence of oscillations, at 2.2σ, independently on the choice of the "500 keV" or "125-250-500 keV" data set, with a reliable Monte Carlo calculation. Hence, it is unclear if the Neutrino-4 anomaly is due to neutrino oscillations or a statistical fluctuation of the data.
A further indication can be obtained by finding the distribution of the best-fit points assuming the absence of oscillations, as done in Figure 1 of Ref. [14] for the STEREO experiment and in Figure 2 of Ref. [23] for a toy reactor neutrino oscillation experiment. Figure 5 shows the distributions of best-fit points that we obtained for the four analysis methods of the Neutrino-4 data that we considered in this paper. One can see that the actual best fit points obtained from the Neutrino-4 data, shown by the blue crosses, lie in regions with high probability. Therefore, they are not unlikely in the absence of oscillations, in spite of their large mixing.
From Figure 5 one can also see that, although the data are generated assuming the absence of mixing, the probability of finding zero-mixing best fit is negligible. Instead, there is a large probability to find large values of the mixing, in the range 0.1 sin 2 2ϑ ee ≤ 1. Table II shows the values of the probabilities in four intervals of sin 2 2ϑ ee . One can see that there is a maximum of about 60-70% of probability to find the best fit for 0.1 < sin 2 2ϑ ee < 0.5, but the probability to find a best fit with sin 2 2ϑ ee > 0.9 is far from negligible, being about 16-18% in the analyses taking into account the energy resolution of the detector. Therefore, it is clear that the large mixing of the best-fit points obtained from the analysis of the Neutrino-4 data is compatible with the absence of oscillations.

V. CONCLUSIONS
In this paper we have presented the results of detailed analyses of the Neutrino-4 data aimed at finding the significance of the large-mixing short-baseline neutrino oscillation signal claimed by the Neutrino-4 collaboration at more than 3σ [11][12][13]. We found that the results of the Neutrino-4 collaboration can be reproduced approximately only by neglecting the effects of the energy resolution of the detector. Including these effects, we found that the best-fit point and the surrounding 1σ allowed region in the (sin 2 2ϑ ee , ∆m 2 41 ) plane lie at even larger values of the mixing than that claimed by the Neutrino-4 collaboration. However, the 3σ allowed region calculated with the standard ∆χ 2 method based on Wilks' theorem is much larger than that claimed by the Neutrino-4 collaboration and include the case of zero mixing, i.e. the absence of oscillations. The corresponding statistical significance of short-baseline neutrino oscillations is about 2.7σ.
We have also shown that there is a strong tension between the Neutrino-4 allowed regions in the (sin 2 2ϑ ee , ∆m 2 41 ) plane and the 95% C.L. exclusion curves of KATRIN [29], PROSPECT [30], STEREO [31], and solar ν e 's. We have further studied the statistical significance of the Neutrino-4 indication in favor of short-baseline neutrino oscillations with a more reliable Monte Carlo evaluation of the distribution of ∆χ 2 = χ 2 −χ 2 min in the (sin 2 2ϑ ee , ∆m 2 41 ) plane. We found that the allowed regions extend to lower values of the mixing, as expected [20,22,23], and the statistical significance of short-baseline neutrino oscillations decreases to about 2.2σ.
We have also shown with a Monte Carlo simulation of a large set of Neutrino-4-like data that it is not unlikely to obtain a best-fit point that has a large mixing, even maximal, in the absence of oscillations. Therefore, we conclude that the claimed Neutrino-4 indication in favor of short-baseline neutrino oscillations with very large mixing is rather doubtful.