Reactor Antineutrino Anomaly Reanalysis in Context of Inverse-Square Law Violation

Abstract

We discuss a possibility that the so-called reactor antineutrino anomaly (RAA), which is a deficit of the ${\overline{\nu }}_e$ rates in the reactor experiments in comparison to the theoretical expectations, can at least in part be explained by applying a quantum field-theoretical approach to neutrino oscillations, which in particular predicts a small deviation from the classical inverse-square law at short (but still macroscopic) distances between the neutrino source and detector. An extensive statistical analysis of the current reactor data on the integrated ${\overline{\nu }}_e$ event rates vs. baseline is performed to examine this speculation. The obtained results are applied to study another long-standing puzzle—gallium neutrino anomaly (GNA), which is a missing ${\nu }_e$ flux from ${}^{37}$Ar and ${}^{51}$Cr electron-capture decays as measured by the gallium–germanium solar neutrino detectors GALLEX and SAGE.

1. Introduction

Nuclear reactors produce a clean and intense flux of electron antineutrinos and thus are very good sources for experiments in neutrino physics. The ${\overline{\nu }}_e$ spectrum is composed of thousands of spectral components formed by the $\beta$ decay of the fission products of four main isotopes, ${}^{235}$U, ${}^{238}$U, ${}^{239}$Pu, and ${}^{241}$Pu (see [1,2,3] for comprehensive reviews). Calculating this spectrum is not an easy and well-defined task. Rather sophisticated (decade ago) calculations [4,5] yielded a net 3–3.5% upward shift in the predicted spectrum-averaged event rate with respect to the previously expected flux used in the earlier short baseline (SBL) reactor experiments (ILL [6,7,8], SRP [9,10,11,12,13,14], Gösgen [15,16], Krasnoyarsk [17,18,19,20], Rovno [21,22,23,24,25,26], Bugey [27,28,29,30,31,32,33]), as well as in the medium and long baseline (MBL, LBL) experiments Palo Verde [34,35,36,37], CHOOZ [38,39,40], and KamLAND [41,42]. The flux normalization uncertainty in the calculation by Mueller et al. [5] was declared to be $\pm 2.7$%. This implies that the measured event rates in the reactor experiments is about 6% less than previously thought. This deficit known as the “reactor antineutrino anomaly” (RAA) still remains an unresolved problem of particle and nuclear physics.

Figure 1 illustrates the current state-of-the-art of the RAA; the early and more recent reactor data in this figure are compared with the prediction based on the ${\overline{\nu }}_e$ flux by Mueller et al. [5]. Here, and below, we use the best-fit values for the neutrino mass-squared splittings and mixing angles from the recent global analysis of the neutrino oscillation data by Esteban et al. [43]; here and below, we assume the normal neutrino mass ordering and no $CP$ violation (irrelevant to the issue under consideration). The inverse $\beta$ decay (IBD) cross section is calculated by using the results of [44] (see Section 6.3 for more details). In Figure 1 and similar plots shown below, all the curves correspond to a reactor with pure ${}^{235}$U fuel; in all subsequent calculations, we explicitly take into account the particular fuel composition in each experiment, although the corresponding effect is very small (see Section 6.5 for explanation). The original data (references are listed in the caption to the figure) are rescaled according to Mention et al. [45] and then renormalized to the new world average value of the neutron mean life [46], as explained in Section 6.1. Note that the RENO and Daya Bay datasets from, respectively, [47,48], are relative measurements; they are simply normalized to the curve to demonstrate that these measurements are in excellent agreement in shape with the standard $3\nu$ oscillation scenario. The newer high-precision measurements [49,50] are absolute (in Figure 1; they are placed at the effective flux-weighted baselines). In order not to complicate the figure, we do not show very recent relative SBL measurements of Neutrino-4 (SM-3 HEU research reactor, Dimitrovgrad), STEREO (ILL HEU research reactor, Grenoble) [51], PROSPECT (HFIR HEU research reactor, Oak Ridge) [52,53], and DANSS (a $3.1$ GW${}_{\mathrm{th}}$ LEU industrial reactor, Kalinin NPP) [54,55], all of which, however, are critical for the present study and will be discussed in detail in Section 6.

It can be seen from the figure that the theory is in very poor agreement with most of the absolute measurements: on average, only about 94–95% of the emitted antineutrino flux was detected at short and medium baselines and maybe even less at very short baselines (the measurements at $L\lesssim 15$ m, where L is the distance between the reactor core and detector). Several exceptions (Nucifer, SRP-II, Palo Verde, CHOOZ) do not formally contradict the general trend within the data uncertainties. Consequently, both early and new measurements at different baselines indicate either “new physics” or incorrect inputs, primarily related to the reactor antineutrino flux calculations. Both these possibilities are currently under intensive experimental and theoretical studies, with the search for light sterile neutrinos as the most popular explanation of the anomaly.

In this article, which is a continuation of our previous studies [61,62,63], we discuss an explanation of the RAA alternative (fully or partially) to the sterile neutrino hypothesis and complementary to the more prosaic one, associated with the difficulties in the ${\overline{\nu }}_e$ flux modeling. Our hypothesis is based on the possibility of violation of the classical (geometric) inverse-square law (ISL) at short but macroscopic distances between the (anti)neutrino source and the detector. The ISL violation (ISLV) is one of the predictions of a field-theoretical approach to neutrino oscillations, which, however, cannot predict the spatial scale at which the ISLV effect may become measurable. We therefore use the reactor data to determine or confine this scale (assuming that the ISLV is related to the RAA). We emphasize that the present study is not a test of the ISLV as a mathematical statement, but only an elucidation of whether the ISLV scale (not determined by the theory) can be related to the SBL reactor data. Since, as we hope to show, the answer to this question is positive, we are also trying to apply the same hypothesis to another still unsolved puzzle of neutrino physics—the so-called “gallium neutrino anomaly” (GNA). However, that issue is somewhat more speculative and the conclusions are less reliable.

2. Extra Neutrinos or Miscalculated Flux?

Most if not all efforts to explain the reactor anomaly rely on the hypothesis of the existence of one or more light (eV mass scale) sterile neutrinos that are fundamental neutral fermions with no standard model interactions, except those induced by mixing with the standard (active) neutrinos (see [64,65,66,67,68,69] for reviews and further references). The active-to-sterile neutrino mixing would lead to a distance-dependent spectral distortion and overall reduction of the reactor ${\overline{\nu }}_e$ flux.

In Figure 2, we show, as an example, the results of calculations carried out in the framework of the simplest “3 + 1” phenomenological model with one sterile (anti)neutrino, ${\nu }_4$, by using the three pairs of the ${\nu }_4-{\nu }_1$ mixing parameters, $(\Delta m_{41}^2,{sin}^{2}2{\theta }_{14})$, listed in the legend of the figure.

These values were derived in [70] from detailed statistical analyses of all the neutrino oscillation data available to date. The “SBL rates only” fit includes the SBL reactor data except the points “Krasnoyarsk-IV” and “Rovno 92” and very new data from Nucifer (at the OSIRIS research reactor, CEA-Saclay), and MiniCHANDLER (Mobile Neutrino Lab deployed at the North Anna NPP), see Figure 1. The “SBL + Bugey 3 spectrum” fit includes the same SBL data set and spectral data from Bugey 3 [32]. The “Global ${\nu }_e$ disappearance” fit involves the data from the reactor experiments (including the spectral data from Bugey 3), as well as solar neutrinos (261 data points from Homestake, SAGE, GALLEX/GNO, Super-Kamiokande, and SNO experiments), radioactive source experiments at SAGE and GALLEX, and the LSND/KARMEN ${\nu }_e$ disappearance data from ${\nu }_e-{}^{12}$C scattering (see [70] for the full list of references and further details). It should be mentioned that these fits operate with somewhat lower (to within roughly 1%) values for the IBD event rates and with a bit different covariance matrix, as compared to those used in the present analysis (see Section 6.1). Moreover, some of the inputs are a little out of date (see, e.g., [71] for a more recent global analysis). This is, however, not very important for our purposes, since the considered “3 + 1” model is not used in further analysis and is only needed to demonstrate that it leads to an effect very similar to a banal renormalization of the ${\overline{\nu }}_e$ flux.

The solid curve in Figure 2 represents the same $3\nu$ oscillation prediction as in Figure 1, but shifted down by the normalization factor $N_0$ derived from a fit to all the reactor data. In this fit, we take into account all known correlations between the data, including the overall normalization uncertainty, which is taken to be 2.7% [45]. The details of such fits are discussed below, in Section 6.5.1. The obtained normalization factor $N_0=0.{953}_{-0.026}^{+0.027}$ (${\chi }^2/\mathrm{NDF}\approx 0.97$) does not contradict the adopted flux uncertainty within ≈2$\sigma$ but is somewhat different from the results of previous calculations [44,45,72], which used different data samples and input parameters. All curves in Figure 2 but “SBL rates only” are in agreement, within the errors, with the new absolute measurements of RENO and Daya Bay, but are in tension with several data points. The “SBL rates only” curve is expectedly in agreement with most of the SBL data but is in slight tension with the Palo Verde, CHOOZ, RENO, and Daya Bay absolute data points.

As first stated in [73], the true uncertainty in the ${\overline{\nu }}_e$ flux predictions may be as large as $5\%$, and the spectral shape uncertainties may be much larger due to a poorly known structure of the forbidden decays. This finding has been in essence confirmed by several new precision measurements and theoretical studies. In the recent MBL experiments RENO [74], Daya Bay [75,76], and Double Chooz [77], using well-calibrated detectors and fairly different industrial reactors, an unexpected excess (so-called “5 MeV bump” or “shoulder”) was observed in the IBD events at ${\overline{\nu }}_e$ energies within the 4.8–7.3 MeV. Similar excess has been also seen in the reactor SBL ($L=23.7$ m) experiment NEOS (Hanbit-5 LEU reactor, Yeonggwang) [78]. Post factum, the same bump was recognized in the data from the earlier experiment, carried out at three distances from the $2.8$ GW reactor at Gösgen NPP [79], and also in a series of experiments performed at five (shorter) distances from the $1.375$ GW WWER-440 reactor at Rovno NPP [22] (see also [80]). The statistical significance in all mentioned experiments is beyond doubt and thus the observed excess appears to be baseline independent. It is currently unclear which physics are responsible for this bump, although several possibilities have been proposed to explain it [81,82,83,84] (see also [85,86,87] for further discussion and references). This problem complicates the study of the IBD spectrum distortions mediated by sterile neutrinos but might have little or no relation to the reactor rate anomaly. The study of the spectral reactor data obtained in the earlier and new experiments goes far beyond the scope of this work. However, it can be inferred from the above that the efforts to explain the reactor anomaly with the sterile neutrino hypothesis may be somewhat premature due to still unsolved (or yet unidentified) problems related to nuclear physics rather than “new physics”.

In this connection, model-independent methods that do not require knowledge of the non-oscillating ${\overline{\nu }}_e$ spectrum at the reactor should be mentioned. Such methods were, in particular, used in the experiment Neutrino-4 [88,89] (see also [90,91,92] for earlier results and further references) and in the recent combined analysis of RENO and NEOS data [93] (since the RENO and NEOS detectors share the same reactor complex, this analysis removes the ${\overline{\nu }}_e$ source dependency in the previous NEOS sterile neutrino search [78]). The results of these experiments provide intriguing evidence in favor of the “3 + 1” scenario, although they appear to be in deep conflict with each other. In the two consecutive analyses of the Neutrino-4 data [88,89] (eprint version 7), the following allowed values for the mass-squared splitting and mixing angle were obtained by using the so-called coherent data summation method:

$$\begin{array}{c}\Delta m_{41}^2=\left(7.25\pm 0.{13}_{}\mathrm{stat}\pm 1.{08}_{}\mathrm{syst}\right){\mathrm{eV}}^2,{sin}^{2}2{\theta }_{14}=0.26\pm 0.{12}_{}\mathrm{stat}\pm 0.{05}_{}\mathrm{syst}\left(3.0\sigma \right),\\ \Delta m_{41}^2=\left(7.30\pm 0.{13}_{}\mathrm{stat}\pm 1.{16}_{}\mathrm{syst}\right){\mathrm{eV}}^2,{sin}^{2}2{\theta }_{14}=0.36\pm 0.{12}_{}\mathrm{stat}\left(2.9\sigma \right),\end{array}$$

These results do not agree with the constraints based on the combined analysis of the MINOS, Daya Bay, and Bugey-3 data [94], and with the limits obtained in the recent SBL experiments DANSS [95,96,97], STEREO [51], NEOS [78] (see also [98], and PROSPECT [53]; this issue is now broadly discussed in the literature [99,100,101,102,103]. On the other hand, the $68\%$ C.L. allowed the region with the best fit of

$$|\Delta m_{41}^2|=\left(2.37\pm 0.03\right){\mathrm{eV}}^2,{sin}^{2}2{\theta }_{14}=0.09\pm 0.03,$$

obtained in [93] formally (with minor reservations), does not contradict the mentioned constraints. Figure 3 shows a comparison of the ratios similar to those shown in Figure 2, but with the sterile neutrino parameters obtained in [88,89,93].

It is seen that the Neutrino-4 results contradict most of the data on the integrated event rate, and especially dramatically differ from the precision absolute measurements of Double Chooz [59], Daya Bay [50], and RENO [49]. Comparison with Figure 2 suggests that the Neutrino-4 results disagree with the outcome of the global-fit analyses of all oscillation experiments. Quite to the contrary, the “RENO + NEOS” best-fit parameters are in agreement with the “global ${\nu }_e$ disappearance” $4\nu$ fit. In addition, lastly, both the “RENO + NEOS” and (to a much greater extent) Neutrino-4 allowed regions that are in strong tension with the cosmological constraints based on “3 + 1” model analyses of the modern CMB and BAO data [104,105].

One can in particular see from Figure 2 and Figure 3 that the proper renormalization of the ${\overline{\nu }}_e$ flux is hardly distinguishable from the “global ${\nu }_e$ disappearance” $4\nu$ fit and almost fully coincides with the “SBL + Bugey 3 spectrum” and “RENO + NEOS” best fits. This coincidence is in fact accidental because the normalization factor is highly dependent on the ${\overline{\nu }}_e$ spectrum model. The shapes of the oscillation curves also depend on the shape of the spectrum and, therefore, on the model of the spectrum, but to a much lesser extent. What is more important is that mixing with the sterile neutrinos of mass $m_4\lesssim 2-3$ eV is very similar to, if not indistinguishable from, an overall downward shift of the curve $R=R\left(L\right)$ (relative to the standard $3\nu$ curve), except very short baselines, where the fine structure of the “3 + 1” curves becomes potentially measurable. To put it differently, the expected effect from sterile neutrinos for the integrated event rate is very similar to an overall ${\overline{\nu }}_e$ flux renormalization, and the latter is defined by the mixing angle ${\theta }_{14}$. Below, the overall normalization factor will be used as an adjustable parameter and will have a double meaning: as a factor compensating the uncertainty in the reactor antineutrino spectrum calculations, or (if desired) as an indicator of the sterile neutrino effect.

It may be even more interesting that the steady decrease of the event rate at very short baselines, mentioned in our previous works [62,63], if real, cannot be described by a flux renormalization alone, as well as by the eV-scale sterile neutrinos with the parameters following from the “global ${\nu }_e$ disappearance” fits. Thus, it seems pertinent to consider an alternative or complementary explanation, namely the ISLV hypothesis. It will be shown below that this explanation is almost insensitive to the features of the ${\overline{\nu }}_e$ energy spectra (of course, when we deal only with the integrated event rates) and thus we can (perhaps temporarily) decouple the nuclear physics problems from possible allusions to the effects of new physics. Finally, we note that even the existence of the light sterile neutrinos does not at all exclude the possibility of the presence of the ISLV effect as well, since both effects may work together.

3. Why Quantum Field Theory?

Before digging into the inverse-square law violation that is based on the QFT approach to neutrino oscillations, let us briefly discuss the standard quantum-mechanical (QM) approach. According to the latter, the probability of the transformation of the neutrino with flavor $\alpha$ to the neutrino with flavor $\beta$ ($\alpha ,\beta =e,\mu ,\tau$) on the way from source to detector, separated by distance L, is described by the following well-known expression:

$${\mathcal{P}}^{\mathrm{QM}}\left({\nu }_{\alpha }\to {\nu }_{\beta }\right)\equiv {\mathcal{P}}_{\alpha \beta }^{\mathrm{QM}}=\sum _{jk}{\left(V_{\alpha k}{V}_{\beta j}\right)}^{*}{V}_{\alpha j}{V}_{\beta k}exp\left(i\frac{\Delta m_{jk}^2}{2E_{\nu }}L\right).$$

(1)

Here, $\Delta m_{jk}^2=m_j^2-m_k^2$, $V_{\alpha j}$ are the elements of the neutrino flavor mixing matrix in vacuum (“PMNS” matrix) defined to link the neutrino mass eigenstates $|{\nu }_j\rangle$ ($j=1,2,3$) with definite masses $m_j$ and momenta ${\mathit{p}}_j$ to the neutrino flavor eigenstates $|{\nu }_{\alpha }\rangle$ ($\alpha =e,\mu ,\tau$) with definite lepton numbers,

$$|{\nu }_{\alpha }\rangle =\sum _j{V}_{\alpha j}^{*}|{\nu }_j\rangle ⟺|{\nu }_j\rangle =\sum _{\alpha }{V}_{\alpha j}|{\nu }_{\alpha }\rangle$$

(both sets of the states are orthonormal: $\langle {\nu }_j|{\nu }_k\rangle ={\delta }_{kj}$, $\langle {\nu }_{\alpha }|{\nu }_{\beta }\rangle ={\delta }_{\alpha \beta }$), and $E_{\nu }$ is the neutrino energy. From here on, we use natural units with $\hslash =c=1$. To obtain Equation (1), one has to use several assumptions that seem intuitively obvious, bordering on commonplace:

(i)

Massive neutrino states originating from reaction or decay of any kind have the same (definite) momentum: ${\mathit{p}}_j={\mathit{p}}_{\nu }$.

(ii)

Neutrino masses $m_j$ are so small that, in essentially all experimental circumstances (or, more precisely, in a wide class of reference frames), the neutrinos are ultrarelativistic ($|{\mathit{p}}_j{|}^2\gg m_j^2$) and hence $E_j=\sqrt{|{\mathit{p}}_j{|}^2+m_j^2}\simeq |{\mathit{p}}_{\nu }|+m_j^2/(2|{\mathit{p}}_{\nu }|)$; this is a cornerstone of the QM theory of neutrino oscillations. One can nevertheless neglect that the flavor states (assumed to be “physical”) have no definite masses and momenta and suppose they have (can be characterized by = can be interpreted as) the common energy $E_{\nu }=|{\mathit{p}}_{\nu }|$.

(iii)

Due to the same reason (ultrarelativism), the travel time T of the neutrino from the source to the detector can be replaced by the distance between them, $T\approx L$.

Formula (1) is widely used in the analyses of all earlier and current experiments studying vacuum neutrino oscillations. However, from the theoretical point of view, the assumptions used for its derivation are rather doubtful.

First of all, if a quantum state has definite momentum, its spatial position is completely undefined owing to Heisenberg’s uncertainty principle. Thus, assumption (i) implies that the spatial coordinates of the neutrino production (${\mathit{X}}_s$) and detection (${\mathit{X}}_d$) are fully uncertain and, as a result, the distance $L=|{\mathit{X}}_d-{\mathit{X}}_s|$ is uncertain too. The same is true for the travel time T, although more cumbersome reasoning is required here. In a more sophisticated theory, the neutrino momentum uncertainty $\delta |{\mathit{p}}_{\nu }|$ must be at least the minimum of the inverse characteristic dimensions of the source and detector devices along the neutrino beam. Consequently, more realistic neutrino states should rather be wave packets (albeit with a small momentum dispersion), depending in the general case on the quantum states of the particles involved in the neutrino production and detection processes, but not the states with definite momentum. Well, one may view this inference as a formalistic cavil, and the definite momentum assumption as a reasonable approximation.

Is the “equal momenta assumption” another reasonable approximation? No, this key assumption is unphysical as it is reference-frame dependent. Indeed, according to the Lorentz transformations,

$$E_j^{\prime }={\Gamma }_{\mathit{u}}\left[E_j-\left(\mathit{u}{\mathit{p}}_{\nu }\right)\right],{\mathit{p}}_j^{\prime }={\mathit{p}}_{\nu }+{\Gamma }_{\mathit{u}}\left[\frac{{\Gamma }_{\mathit{u}}}{{\Gamma }_{\mathit{u}}+1}\left(\mathit{u}{\mathit{p}}_{\nu }\right)-E_j\right]\mathit{u},{\Gamma }_{\mathit{u}}=\frac{1}{\sqrt{1-{\mathit{u}}^2}},$$

if the equality ${\mathit{p}}_j={\mathit{p}}_{\nu }$ holds true (for all j) in some reference frame, then it is not true in another frame moving with the velocity $\mathit{u}$ relative to the first one, namely,

$${\mathit{p}}_j^{\prime }-{\mathit{p}}_k^{\prime }=\left(E_k^{\prime }-E_j^{\prime }\right)\mathit{u}={\Gamma }_{\mathit{u}}\left(E_k-E_j\right)\mathit{u}\approx \frac{\left(m_k^2-m_j^2\right)\mathit{u}}{2|{\mathit{p}}_{\nu }|\sqrt{1-{\mathit{u}}^2}}.$$

Hence, for relativistic frames (typical in, e.g., astrophysical environments), violation of condition (i) may lead to a large phase shift, considering that the oscillation phase is approximately invariant ($L^{\prime }/E_{\nu }^{\prime }\approx L/E_{\nu }$) as ${\Gamma }_{\mathit{u}}\ll min(E_j/m_j)$ and approximations (ii) and (iii) are valid. Neutrino masses are so small that even very relativistic boosts (${\Gamma }_{\mathit{u}}\gg 1)$) leave neutrinos ultrarelativistic, that is, condition (ii) is satisfied, but condition (i) is clearly violated. Of course, this violation is negligible if the velocity $\mathit{u}$ is nonrelativistic, but, in general, this is not the case. Interpreting the Lorentz transformation as active, we see that condition (i) cannot be satisfied simultaneously for all neutrinos emitted due to inelastic collisions or decay of particles moving at different velocities (from subrelativistic to ultrarelativistic), as, e.g., in a broad-spectrum pion beam from an accelerator, relativistic astrophysical jet, or cosmic-ray collisions with the Earth atmosphere. As a result, we conclude that the equal momenta assumption is, strictly speaking, meaningless.

Moreover, even if condition (i) is satisfied approximately in some circumstances (e.g., for a narrow-energy beam of parent particles), the velocities of neutrinos of different masses are different:

$$v_j=|{\mathit{v}}_j|=\frac{|{\mathit{p}}_{\nu }|}{E_j}\approx 1-\frac{m_j^2}{2|{\mathit{p}}_{\nu }{|}^2}⟹v_j-v_k\approx \frac{m_k^2-m_j^2}{2E_{\nu }^2}.$$

During the time T, the neutrino ${\nu }_j$ travels the distance $L_j=|{\mathit{v}}_j|T$. As the neutrino velocities are different, there must be a spread in distances of each neutrino pair:

$$\delta L_{jk}=\left|L_j-L_k\right|\approx \frac{\left|\Delta m_{jk}^2\right|}{2E_{\nu }^2}L,$$

where we put (denote) $L=T$. If, for example, $L=1$ AU and $E_{\nu }=10$ MeV, then $\delta L_{12}\approx 5.6\times {10}^{-6}$ cm. Is such a spread small or large to preserve the quantum coherence necessary for neutrino flavor transitions? If one considers the physical neutrino states as wave packets, it is obvious that the oscillations’ pattern must disappear when the neutrino packets diverge far enough from each other. However, how far? In other words: how the effective size of the neutrino wave packet depends on energy and on features of the neutrino production and detection processes? The standard QM approach does not provide information on this.

Another objection to the QM approach relates to hypothetical heavy (sterile or superweakly interacting) neutrinos. Consider a toy experiment in which a neutrino beam is generated by decay at flight of very high-energy particles P whose mass is less than the mass of heavy neutrino and, therefore, the latter cannot be produced in decays of particles P (regardless of the value of their coupling with P). However, according to the naive QM approach, which ignores the physics of neutrino production and detection, the heavy neutrinos can appear in the beam through mixing with the ordinary light neutrinos and, moreover, be ultrarelativistic. The Lorentz boost into the rest frame of particle P immediately shows the impossibility of such a phenomenon, but the standard QM theory does not prohibit this. It is worth noting that this objection is not in the least directed against the existence of heavy neutrinos.

There are many other, sometimes interrelated questions that cannot be answered within the QM approach, for example:

•

Do charged leptons oscillate, and if not (as experiment suggests), then why?

•

Do relic (“C$\nu$B”) neutrinos (some of which are definitely nonrelativistic) oscillate? Or, in a more general form: how to extend the theory to cover nonrelativistic neutrinos?

•

How does the relative motion of the neutrino source and detector affect the survival and transition probabilities?

Further discussions of some of these and related issues and numerous references can be found in [106,107,108,109]. In the following consideration, we adhere to an approach based on perturbative quantum field theory (QFT), which explicitly accounts for the neutrino production and detection processes, does not use the ad hoc assumptions of the QM approach and is free of paradoxes. It allows one to interpret the QM formula (1) properly, determine its area of applicability, and derive corrections to this formula. It also makes it possible to answer the above questions, predicts potentially observable effects (e.g., loss of coherence and dispersion distorting the standard QM neutrino oscillation pattern), and leaves room for further generalizations and extensions. Various aspects of the QFT approach have been discussed by many authors [109,110,111,112,113,114,115,116,117,118,119,120,121,122,123,124,125,126,127,128,129,130,131,132,133,134,135,136,137,138,139,140,141,142,143] (see also references in these article), who used rather different methods and approximations, leading to similar (although not always identical) conclusions and sometimes to disputes (see, e.g., [144,145,146]). The experimental activity in the field is still at the very beginning [147], but the reactor experiments (main subject of the present article) have the potential to study the decoherence and dispersion effects predicted by the QFT approach [148,149,150,151].

Below, we will very briefly review the essentials of the covariant “diagrammatic” formalism proposed in [134,135] (see also [109] for a more detailed presentation) and dwell on another prediction of the theory, not related to decoherence or dispersion. The key ingredient of the formalism are covariant wave packets used to describe the quantum states of the initial and final particles involved in the neutrino production and detection processes. Many studies have been devoted to the development of this aspect of the theory (see, e.g., [152,153,154,155,156,157] and references therein), which extends well beyond the narrow topic of neutrino oscillations.

4. A Sketch of the QFT Approach

The “neutrino oscillation” phenomenon in the S-matrix QFT approach is nothing else than a result of interference of the macroscopic Feynman diagrams perturbatively describing the lepton number violating processes with the massive neutrino eigenfields that are internal lines (propagators) connecting vertices of interactions in which the neutrino is produced (together with a charged lepton) and absorbed (producing another charged lepton). In other words, the massive neutrinos, ${\nu }_j$ ($j=1,2,3$), are treated as virtual intermediate states, while the neutrinos with definite flavors, ${\nu }_{\alpha }$ ($\alpha =e,\mu ,\tau$), do not participate in the formalism at all, although they can be used to compare the predictions with those from other approaches, in particular, with the QM approach based on the concept of flavor mixing. The spatial interval between the Feynman diagram vertices (hereinafter referred to as “source” and “detector”) can be arbitrarily large.

A generic example of such a macrodiagram is shown in Figure 4, which also introduces notation used below. The external lines of the macrodiagrams are assumed to be asymptotically free quasi-stable wave packets (WP) rather than the conventional one-particle Fock’s states $|\mathit{k},s\rangle$ with definite 3-momenta $\mathit{k}$ and spin projections s.

If the states $I_{s,d}$ and $F_{s,d^{\prime }}$ consist exclusively of hadronic WPs and $\alpha \ne \beta$, the lepton charges $L_{\alpha }$ and $L_{\beta }$ are not conserved in the process $I_s\oplus I_d\to F_s^{\prime }{\ell }_{\alpha }^{+}\oplus F_d^{\prime }{\ell }_{\beta }^{-}$ (where sign “⊕” indicates that the regions of interaction of the corresponding WPs are macroscopically separated in space and/or time) that is only possible via exchange of massive neutrinos (no matter whether they are Dirac or Majorana particles). Formally, the sets of initial (I) and final (F) states interacting in the vertices of the macrodiagram may include any number of WPs (particles or nuclei), but more often the I states include one or two WPs as, e.g., in the process ${\pi }^{+}\oplus n\to {\mu }^{+}\oplus {\tau }^{-}p$. The particular “decryption” of the neutrino production/absorption mechanism in Figure 4 assumes the standard model charged current interaction of quarks and leptons, although this is not necessary for the formalism under discussion; the more general case can include nonstandard interactions (e.g., flavour-changing neutral currents), new generation neutrinos and charged leptons, and so on.

According to [135], the free external WP states are constructed as covariant space-time point dependent linear superpositions of the one-particle states,

$$|\mathit{p},s,x\rangle =\int \frac{d\mathit{k}\varphi (\mathit{k},\mathit{p})e^{i(k-p)x}}{{\left(2\pi \right)}^{3}2E_{\mathit{k}}}|\mathit{k},s\rangle$$

(2)

(where p and k are the 4-momenta; $k=(E_{\mathit{k}},\mathit{k})$, $E_{\mathit{k}}=\sqrt{{\mathit{k}}^2+m^2}$, and $x=(x_0,\mathit{x})$), satisfying the correspondence principle which demands that $|\mathit{p},s,x\rangle$ turns into $|\mathit{k},s\rangle$ in the plane-wave limit (PWL), $|\mathit{p},s,x\rangle \stackrel{\mathrm{PWL}}{⟼}|\mathit{k},s\rangle$, which is equivalent to the following condition for the relativistic invariant form factor function $\varphi$:

$$\varphi (\mathit{k},\mathit{p})\stackrel{\mathrm{PWL}}{⟼}{\left(2\pi \right)}^{3}2E_{\mathit{p}}\delta (\mathit{k}-\mathit{p}).$$

(3)

The detail properties of the WP states Equation (2) are discussed in [109,153]. Very shortly, the function $\varphi (\mathit{k},\mathit{p})$ can parametrically depend on a set of constants or Lorenz-invariant combinations of “hidden” variables (4-momenta of the WP states participated in the creation of the state Equation (2)), but in the simplest case, the momentum spread of the packet Equation (2) can be characterized by a single parameter $\sigma$ (momentum spread) defined in such a way that the limit Equation (3) occurs as $\sigma =0$. Thus, it is natural to assume $\sigma$ to be small, namely $0\le \sigma \ll m$ (WP states for massless particles require a more special consideration). In this simplest case, the WP is spherically symmetric in its center of inertia frame of reference (where $\mathit{p}=0$) and, in arbitrary RF, the vector $\mathit{p}$ has the physical meaning of the most probable 3-momentum.

In our approach, the amplitude of the process $I_s\oplus I_d{\ell }_{\alpha }^{+}\to F_s\oplus F_d{\ell }_{\beta }^{-}$ is constructed as ordinary,

$${\mathcal{A}}_{\beta \alpha }\stackrel{\mathrm{def}}{=}\frac{\langle \mathrm{out}|\mathbb{S}|\mathrm{in}\rangle }{{\left(\langle \mathrm{in}|\mathrm{in}\rangle \langle \mathrm{out}|\mathrm{out}\rangle \right)}^{1/2}},$$

(4)

but the standard Feynman rules have to be modified to account that the standard initial and final one-particle states with definite momenta have to be replaced by the WP states Equation (2). Then, after applying several natural model assumptions and technical simplifications, it is proved [134,135] that the neutrino-induced event rate in an ideal detector can be expressed (somewhat symbolically) in the following form:

$$\frac{dN}{d\tau }=\frac{1}{V_{\mathcal{D}}{V}_{\mathcal{S}}}{\int }_{V_{\mathcal{S}}}d\mathit{x}{\int }_{V_{\mathcal{D}}}d\mathit{y}\int d{\mathfrak{F}}_{\nu }\int d{\sigma }_{\nu \mathcal{D}}{\mathcal{P}}_{\alpha \beta }\left(E_{\nu },\left|\mathit{y}-\mathit{x}\right|\right).$$

(5)

Here, $\tau$ is the detector exposure time, $E_{\nu }$ is the neutrino energy, ${\mathcal{P}}_{\alpha \beta }$ is the QFT generalization of the standard quantum-mechanical neutrino flavor transition probability, and the differential form $d{\sigma }_{\nu \mathcal{D}}$ represents the differential cross section of the neutrino scattering from the whole detector device; $d{\mathfrak{F}}_{\nu }$ is the differential neutrino flux incident on the detector from a stationary source device (e.g., a fission reactor core). The integrations in Equation (5) are over the source ($\mathcal{S}$) and detector ($\mathcal{D}$) fiducial volumes $V_{\mathcal{S}}$ and $V_{\mathcal{D}}$. In the conventional ultrarelativistic approximation, the generalized “flavor transition probability” includes the standard QM oscillation phase factor and corrections of different kinds to it, but the derivation of these results does not require the unphysical assumptions discussed in Section 3. The main requirements explicitly used in the derivation of Equation (5) are, instead, rather natural and apparent:

  (i)

spatial dimensions of $\mathcal{S}$ and $\mathcal{D}$ along the (anti)neutrino beam are small compared to the distance between them but large compared to the effective dimensions ($\sim {\sigma }_{\varkappa }^{-1}$) of all WPs $\varkappa$ colliding, decaying, or appearing in $\mathcal{S}$ and $\mathcal{D}$;

 (ii)

the statistical distributions of the incoming WPs $(a\in I_{s,d})$ over their mean 3-momenta, discrete quantum numbers, and mean spatial coordinates in $\mathcal{S}$ and $\mathcal{D}$ can be described by stationary one-particle distribution functions or, more generally, by density matrices;

(iii)

the spatial distributions of the incoming WPs are smooth, that is, the scale of their variations within $\mathcal{S}$ and $\mathcal{D}$ is much larger than the effective dimensions of the WPs themselves;

 (iv)

the edge effects can be neglected, and the time intervals needed to turn $\mathcal{S}$ and $\mathcal{D}$ on and off are small compared to the periods of their steady operation;

  (v)

the experiment catches only the secondaries b arising in $\mathcal{D}$, whereas the background events caused by the (long-range) secondaries $b\in I_s$ falling into $\mathcal{D}$ from $\mathcal{S}$ can be fully ignored;

 (vi)

the neutrino oscillation lengths $L_{jk}=4\pi E_{\nu }/|\Delta m_{jk}^2|$ are large compared to ${\sigma }_{\varkappa }^{-1}$.

Not all of these constraints are mutually independent and are not equally important. It should also be emphasized that essentially all of them are not necessary in the general formalism but are usually obeyed (or are believed to be obeyed) in most (anti)neutrino oscillation experiments (see, however, Section 7).

The theory explicitly predicts that, at sufficiently long spaces between the neutrino production and absorption points $\mathit{x}$ and $\mathit{y}$, the neutrino flux decreases with increasing $L=|\mathit{y}-\mathit{x}|$ in compliance with the classical inverse-square law (ISL):

$$d{\mathfrak{F}}_{\nu }\propto {|\mathit{y}-\mathit{x}|}^{-2}.$$

(6)

This reasonably expected result is unrelated to the lepton numbers violation; it was derived by using a very general theorem, called the Grimus–Stockinger (GS) theorem [112], which defines the asymptotic behaviour of the amplitude Equation (4) at $L\to \infty$, and this is the crucial point in the context of the problem under investigation. As it follows from the formalism, the L-dependence of the amplitude described by the macrodiagram, like that shown in Figure 4, is defined by the neutrino propagator modified (“dressed”) by the external WPs,

$$G_j\left(x\right)=\int \frac{d^{4}q}{{\left(2\pi \right)}^4}\frac{{\tilde{\delta }}_s(q-q_s){\tilde{\delta }}_d(q+q_d)(\widehat{q}+m_j)e^{-iqx}}{q^2-m_j^2+iϵ},$$

(7)

where $x=(y_0-x_0,\mathit{y}-\mathit{x})\equiv (T,\mathit{L})$, $q_s$ and $q_d$ are the 4-momentum transfers,

$$q_s=\sum _{a\in I_s}{p}_a-\sum _{b\in F_s}{p}_b,q_d=\sum _{a\in I_d}{p}_a-\sum _{b\in F_d}{p}_b,$$

(8)

$p_{\varkappa }$ are the most probable (on-shell) 4-momenta of the external packets $\varkappa \in I_s\oplus I_d\oplus F_s\oplus F_d$, and $m_j$ is the mass of the neutrino mass eigenfield${\nu }_j$. The functions ${\tilde{\delta }}_s(q-q_s)$ and ${\tilde{\delta }}_d(q+q_d)$ are the “smeared” $\delta$ functions parametrically dependent on the (most probable) 4-momenta $p_{\varkappa }$, masses $m_{\varkappa }$ ($m_{\varkappa }^2=p_{\varkappa }^2$), and momentum dispersions ${\sigma }_{\varkappa }$ of the external in and out WPs; it is assumed that $0\le {\sigma }_{\varkappa }^2⋘m_{\varkappa }^2$. The explicit form of the functions ${\tilde{\delta }}_{s,d}$ is determined by the Lorenz-invariant form factor functions ${\varphi }_{\varkappa }={\varphi }_{\varkappa }({\mathit{k}}_{\varkappa },{\mathit{p}}_{\varkappa };{\sigma }_{\varkappa })$ which describe the external WPs, but, regardless of the explicit form of the functions ${\varphi }_{\varkappa }$, in the plain-wave limit (${\sigma }_{\varkappa }\to 0$, $\forall \varkappa$), the functions ${\tilde{\delta }}_{s,d}$ should turn into the ordinary $4D$ Dirac $\delta$ functions,

$${\tilde{\delta }}_s\left(q-q_s\right)\stackrel{\mathrm{PWL}}{⟼}\delta \left(q-q_s\right),{\tilde{\delta }}_d\left(q+q_d\right)\stackrel{\mathrm{PWL}}{⟼}\delta \left(q+q_d\right),$$

thus leading to the exact energy–momentum conservation in the vertices of the macrodiagram shown in Figure 4; as a result, the function Equation (7) becomes, up to a multiplier, the standard fermion propagator.

If, however, the momentum spreads ${\sigma }_{\varkappa }$ are finite, and the space-time behavior of the function Equation (7) is nontrivial. In particular, its spatial dependence at sufficiently large distances L is given by the above-mentioned GS theorem [112], according to which

$$J(\mathit{L},\kappa )=\int \frac{d\mathit{q}}{{\left(2\pi \right)}^3}\frac{\Phi \left(\mathit{q}\right)e^{i\mathit{qL}}}{{\mathit{q}}^2-{\kappa }^2-iϵ}=\frac{e^{i\kappa L}\Phi (-\kappa \mathit{l})}{4\pi L}\left[1+\mathcal{O}\left(\frac{1}{\sqrt{L}}\right)\right]\left(\mathit{l}=\frac{\mathit{L}}{L}\right)$$

as $L=|\mathit{L}|\to \infty$. This offers the QFT explanation of the ISL behavior Equation (6) but does not, however, provide the spatial scale above which the distance L may be considered as “sufficiently large”. It is assumed that the complex-valued function $\Phi \left(\mathit{q}\right)$ itself and its first and second derivatives decrease at least like $1/{|\mathit{q}|}^2$ as $|\mathit{q}|\to \infty$ and ${\kappa }^2>0$. In [61], an extended version of the GS theorem has been proved, which parametrically defines such a scale by using the asymptotic expansion of the integral $J(\mathit{L},\kappa )$ in terms of inverse powers of L at large L. To be more precise, the theorem in its simplest form states that, for any function $\Phi \left(\mathit{q}\right)$ in the Schwartz space $S\left({\mathbb{R}}^3\right),$

$$J(\mathit{L},\kappa )=\frac{e^{i\kappa L}}{4\pi L}{\left[\Phi \left(\mathit{q}\right)+\sum _{n\ge 1}{\displaystyle \frac{{(-i)}^n{D}_n\Phi \left(\mathit{q}\right)}{L^n}}\right]}_{\mathit{q}=-\kappa \mathit{l}},L\to \infty ,$$

(9)

where $D_n$ are recursively defined differential operators in the momentum space; the lowest order operators are

$$\begin{array}{cc}\hfill D_1=& \frac{\kappa }{2}\left[{\mathbf{\nabla }}_{\mathit{q}}^2-{\left(\mathit{l}{\mathbf{\nabla }}_{\mathit{q}}\right)}^2\right]-\left(\mathit{l}{\mathbf{\nabla }}_{\mathit{q}}\right),\hfill \\ \hfill D_2=& \frac{{\kappa }^2}{8}{\left[{\mathbf{\nabla }}_{\mathit{q}}^2-{\left(\mathit{l}{\mathbf{\nabla }}_{\mathit{q}}\right)}^2\right]}^2-\kappa \left(\mathit{l}{\mathbf{\nabla }}_{\mathit{q}}\right)\left[{\mathbf{\nabla }}_{\mathit{q}}^2-{\left(\mathit{l}{\mathbf{\nabla }}_{\mathit{q}}\right)}^2\right]-\frac{1}{2}\left[{\mathbf{\nabla }}_{\mathit{q}}^2-3{\left(\mathit{l}{\mathbf{\nabla }}_{\mathit{q}}\right)}^2\right],\hfill \\ \hfill D_3=& \frac{{\kappa }^3}{48}{\left[{\mathbf{\nabla }}_{\mathit{q}}^2-{\left(\mathit{l}{\mathbf{\nabla }}_{\mathit{q}}\right)}^2\right]}^3-\frac{3{\kappa }^2}{8}\left(\mathit{l}{\mathbf{\nabla }}_{\mathit{q}}\right){\left[{\mathbf{\nabla }}_{\mathit{q}}^2-{\left(\mathit{l}{\mathbf{\nabla }}_{\mathit{q}}\right)}^2\right]}^2-\frac{3\kappa }{8}\left[{\mathbf{\nabla }}_{\mathit{q}}^2-{\left(\mathit{l}{\mathbf{\nabla }}_{\mathit{q}}\right)}^2\right]\hfill \\ \hfill & \times \left[{\mathbf{\nabla }}_{\mathit{q}}^2-5{\left(\mathit{l}{\mathbf{\nabla }}_{\mathit{q}}\right)}^2\right]+\frac{1}{2}\left(\mathit{l}{\mathbf{\nabla }}_{\mathit{q}}\right)\left[3{\mathbf{\nabla }}_{\mathit{q}}^2-5{\left(\mathit{l}{\mathbf{\nabla }}_{\mathit{q}}\right)}^2\right],\hfill \\ \hfill D_4=& \frac{{\kappa }^4}{384}{\left[{\mathbf{\nabla }}_{\mathit{q}}^2-{\left(\mathit{l}{\mathbf{\nabla }}_{\mathit{q}}\right)}^2\right]}^4-\frac{{\kappa }^3}{12}\left(\mathit{l}{\mathbf{\nabla }}_{\mathit{q}}\right){\left[{\mathbf{\nabla }}_{\mathit{q}}^2-{\left(\mathit{l}{\mathbf{\nabla }}_{\mathit{q}}\right)}^2\right]}^3-\frac{{\kappa }^2}{8}\left[{\mathbf{\nabla }}_{\mathit{q}}^2-{\left(\mathit{l}{\mathbf{\nabla }}_{\mathit{q}}\right)}^2\right]\hfill \\ \hfill & \times \left[{\mathbf{\nabla }}_{\mathit{q}}^2-7{\left(\mathit{l}{\mathbf{\nabla }}_{\mathit{q}}\right)}^2\right]-\frac{\kappa }{2}\left(\mathit{l}{\mathbf{\nabla }}_{\mathit{q}}\right)\left[{\mathbf{\nabla }}_{\mathit{q}}^2-{\left(\mathit{l}{\mathbf{\nabla }}_{\mathit{q}}\right)}^2\right]\left[3{\mathbf{\nabla }}_{\mathit{q}}^2-7{\left(\mathit{l}{\mathbf{\nabla }}_{\mathit{q}}\right)}^2\right]\hfill \\ \hfill & +\frac{1}{8}\left[3{\mathbf{\nabla }}_{\mathit{q}}^4-30{\mathbf{\nabla }}_{\mathit{q}}^2{\left(\mathit{l}{\mathbf{\nabla }}_{\mathit{q}}\right)}^2+35{\left(\mathit{l}{\mathbf{\nabla }}_{\mathit{q}}\right)}^4\right].\hfill \end{array}$$

In [158], a closed formula for the $1/L$ expansion of $J(\mathit{L},\kappa )$ was obtained, which is equivalent to Equation (9) but is a bit more mathematically transparent. We do not discuss this result, since below we will only use the leading order (LO) terms of expansion Equation (9), which is sufficient for our present purposes. An analysis of Equation (9) shows that the $1/L$ behavior of the amplitude Equation (4) (and thus the ISL behavior of the event rate) is violated at the distances $L\lesssim {\mathfrak{L}}_0$, where

$${\mathfrak{L}}_0\sim \kappa {\Sigma }_{\mathrm{eff}}^{-2}\approx 20\left(\frac{\kappa }{1\mathrm{MeV}}\right){\left(\frac{{\Sigma }_{\mathrm{eff}}}{1\mathrm{eV}}\right)}^{-2}\mathrm{cm}$$

(10)

and the function ${\Sigma }_{\mathrm{eff}}={\Sigma }_{\mathrm{eff}}(\kappa ;\{{\mathit{p}}_{\varkappa },m_{\varkappa },{\sigma }_{\varkappa }\})$ represents a cumulative effect of the overlapping of the external (in and out) states $\varkappa$ dependent on the mean velocities, masses, and momentum spreads of these states, and on the mean neutrino momentum $\kappa$, defined by these quantities. The explicit form of the function ${\Sigma }_{\mathrm{eff}}$ can be found after specification of a particular model for the external WP states. A simple example is discussed in [61] within the so-called contracted relativistic Gaussian packet (CRGP) model [152,153]. It is in particular shown that the “overlap function” ${\Sigma }_{\mathrm{eff}}$ is defined through the space-time and transverse (with respect to the neutrino propagation direction $\mathit{l}$) components of the so-called inverse overlap tensors (${\tilde{\Re }}_s$ and ${\tilde{\Re }}_d$), which determine the effective space-time ($4D$) overlap volumes of the in and out WP states around the impact points $X_s$ and $X_d$ in the vertices of the macrodiagram Equation (Figure 4). The most general formulas for the tensors ${\tilde{\Re }}_{s,d}$ are derived in [61] (see also [109] for more details). It is significant that they are nearly independent of the masses of ordinary neutrinos (assuming these masses to be small with respect to the neutrino energy and thus $\kappa \simeq E_{\nu }$). Within the CRGP model, it can also be shown that the magnitude of ${\Sigma }_{\mathrm{eff}}$ is strongly affected by the hierarchy of the momentum spreads (${\sigma }_{\varkappa }$) of the external WSs $\varkappa$, but in a simple (though not very realistic) case, when all these spreads are similar in order of magnitude, the overlap function ${\Sigma }_{\mathrm{eff}}$ can also be of the same order, up to the Lorentz factors of the most relativistic particles. Generally, the spatial scale Equation (10) can be macroscopically large for sufficiently small values of the function ${\Sigma }_{\mathrm{eff}}$ and/or for sufficiently high neutrino energies, thus leading to a potentially measurable ISL violation (ISLV). Verification of this possibility is the main goal of the subsequent analysis of the reactor antineutrino data.

It is shown in [152] that the series Equation (9) modifies the formula for the event rate Equation (5) in such a way that the differential neutrino flux Equation (6) is multiplied by the asymptotic series in even powers of L,

$$\begin{array}{c}\hfill 1+\sum _{n\ge 1}\frac{{\mathfrak{C}}_n}{{|\mathit{y}-\mathit{x}|}^{2n}},\end{array}$$

(11)

with the coefficient functions ${\mathfrak{C}}_n$ explicitly defined from the expansion Equation (9). A critical property of the series Equation (11) is that its first term is negative. The proof of this property is not very simple, and one of the most complicated steps in this proof is the proper integration with respect to $q_0$, which is required to obtain the modified neutrino propagator $G_j\left(x\right)$ defined by Equation (7). Using the CRGP model and the third order saddle-point asymptotic expansion (see, e.g., [159]), it can be proved that

$${\mathfrak{C}}_1=-\frac{E_{\nu }^2}{{\Sigma }_{\mathrm{eff}}^4},\frac{1}{{\Sigma }_{\mathrm{eff}}^4}=\frac{1}{8}\left[{\rho }^2-\frac{2\kappa }{\mathfrak{r}}+\frac{{\left({\upsilon }_1^2+{\upsilon }_2^2\right)}^2}{{\mathfrak{r}}^2}\right]>0,$$

(12)

where

$$\begin{array}{c}\hfill \begin{array}{cc}\hfill \mathfrak{r}=& {\tilde{\Re }}_{00}-2{\tilde{\Re }}_{03}+{\tilde{\Re }}_{33},\hfill \\ \hfill {\rho }^2=& {\tilde{\Re }}_{11}^2+2{\tilde{\Re }}_{12}^2+{\tilde{\Re }}_{22}^2,\hfill \end{array}\begin{array}{cc}\hfill \kappa =& {\upsilon }_1^2{\tilde{\Re }}_{11}+2{\upsilon }_1{\upsilon }_2{\tilde{\Re }}_{12}+{\upsilon }_2^2{\tilde{\Re }}_{22},\hfill \\ \hfill {\upsilon }_i=& {\tilde{\Re }}_{0i}-{\tilde{\Re }}_{3i}(i=1,2),\hfill \end{array}\end{array}$$

$${\tilde{\Re }}^{\mu \nu }={\tilde{\Re }}_s^{\mu \nu }+{\tilde{\Re }}_d^{\mu \nu },$$

and ${\tilde{\Re }}_{s,d}^{\mu \nu }$ are the components of the aforementioned inverse overlap tensors in the source and detector vertices; the purely transversal term ${\rho }^2/8$ in Equation (12) is the first-order contribution; the second and third order corrections are not in general small and describe nontrivial effects of the in-in, out-out, and in-out WP overlaps in time and space in both source and detector vertices. It should be pointed out that the scale of the functions ${\tilde{\Re }}_{s,d}^{\mu \nu }$ is defined not only by the momentum spreads of the external WPs, ${\sigma }_{\varkappa }$, but also by their masses and momenta, some of which can be ultrarelativistic. Therefore, the value of the function ${\Sigma }_{\mathrm{eff}}$ may generally differ from any of ${\sigma }_{\varkappa }$ by orders of magnitude. We also note that expression Equation (12) is derived for the coordinate system whose third axis is directed along the neutrino velocity and thus the boost covariance is not explicit. Note that the Lorenz-invariant function

$${\Sigma }_j=1/\left({\Gamma }_j\sqrt{\mathfrak{r}}\right)=m_j{\left({\tilde{\Re }}^{\mu \nu }{p}_{\mu }{p}_{\nu }\right)}^{-1/2}$$

represents (at long but not very long distances) the momentum spread of the effective WP of neutrino ${\nu }_j$, which, in turn, determines the effects of decoherence (see [61,134] and also [154] for rigorous consideration). This, in particular, means that the function ${\Sigma }_{\mathrm{eff}}$ is only mediately related to the neutrino momentum dispersion and thus it cannot be measured at long distances. The ISLV effect is ultimately associated with incomplete overlap in space and time of the noncollinearly interacting incoming and outgoing WPs (e.g., parent and daughter nuclei) involved in the processes of neutrino production and absorption, which must lead to a decrease of the cumulative probability of the WP interaction. The negativeness of the coefficient function ${\mathfrak{C}}_1$ is therefore quite expected.

Let us add that a similar reason (partial overlap of wave packets) is basically responsible also for the very long-range decoherence. In this case, due to a kind of duality, the behavior of the neutrino propagator is naturally translated into the language of the effective neutrino WPs; the effect consists of a distortion of the standard oscillation pattern, and—in the very long distance limit—to the disappearance of the oscillations. The explanation in this language is rather apparent: WPs of neutrinos with different masses move at different mean velocities and, therefore, the probability of their overlap decreases with time and distance. Simultaneously, the neutrino WPs spread with time, thereby increasing the probability of their mutual overlapping. These two competing processes govern the evolution of the oscillation pattern at very long spatial ranges and definitely must work at extraterrestrial (astrophysical) distances.

At short distances (as in the short baseline reactor experiments), the ordinary neutrino oscillations by themselves do not play any role and therefore the ISLV effect has nothing to do with the neutrino masses and mixing, as well as, and even more so, with the overlap or spreading of the effective neutrino WPs. The single macrodiagram Equation (Figure 4) describes the cumulative, “already accomplished” result of the interaction of the external WPs in its vertices, and the counting rate Equation (5) is the result of averaging of the squared amplitude Equation (4) over the periods of operation of the source and exposure of the detector, and these periods bear no relation to the time of neutrino propagation from the source to the detector (this is, by the way, one of the main differences of our formalism from others). Consequently, the spatial scale of the ISLV does not have to be comparable with the effective size of the neutrino WP.

Using Equation (11) in the first approximation (thereby assuming that ${\left|\mathit{y}-\mathit{x}\right|}^2\gg |{\mathfrak{C}}_2/{\mathfrak{C}}_1|$ and hence the ISLV correction is small), substituting it into Equation (5), and taking into account the inequality ${\mathfrak{C}}_1<0$, we arrive at the modified formula for the neutrino event rate,

$$\frac{dN}{d\tau }=\frac{1}{V_{\mathcal{D}}{V}_{\mathcal{S}}}{\int }_{V_{\mathcal{S}}}d\mathit{x}{\int }_{V_{\mathcal{D}}}d\mathit{y}\int d{\mathfrak{F}}_{\nu }\int d{\sigma }_{\nu \mathcal{D}}\left(1-\frac{{\mathfrak{L}}_0^2}{{\left|\mathit{y}-\mathit{x}\right|}^2}\right){\mathcal{P}}_{\alpha \beta }\left(E_{\nu },\left|\mathit{y}-\mathit{x}\right|\right),$$

(13)

which represents the phenomenological signature of the ISLV effect. Needless to say, at present, the function ${\mathfrak{L}}_0$ cannot be derived from the first-principle calculations, but it can be measured (if it is not too small) in the spectral measurements in the experiments with appropriately short but macroscopic baselines. In further analysis, we assume that the decoherence and dispersion effects that are present in the general QFT expression for the survival probability ${\mathcal{P}}_{ee}$ are insignificant at the distances under consideration, and thus the standard QM formula (1) is applicable. Additional simplifications and concretization of formula (13) are discussed below in Section 6.3.

5. Antineutrino Energy Spectra

Reactor antineutrino fluxes are one of the main ingredients and the main source of uncertainties of subsequent analysis. In a typical nuclear reactor, almost all (>99%) antineutrinos are produced through thousands of ${\beta }^{-}$ decay branches of fission fragments from ${}^{235}$U, ${}^{238}$U, ${}^{239}$Pu, and ${}^{241}$Pu. There are two main approaches to calculate the antineutrino fluxes produced by these isotopes. One can either employ the ab initio method [81,160,161,162,163,164,165,166,167,168,169,170,171,172,173] by a direct summation of the ${\overline{\nu }}_e$ energy spectra form the $\beta$ decay of each fission fragment using information from the relevant nuclear databases, or apply a conversion procedure [4,174,175,176,177,178,179,180,181,182,183] based on measurements of the integral electron energy spectra for the fissionable isotopes; the conversion-based models have relatively little dependence on the nuclear databases. Both approaches can be used in a “summation + conversion” combination, as in [5,184]. In [80], a direct method for measuring the IBD reaction products was applied (using the high-statistics data collected with the high-aperture RONS spectrometer at the Rovno NPP [185]).

It might be good to point out here that, if the ISLV hypothesis turns out to be true, then the conversion method based on the $\beta$-spectrum measurements at very short distances from the reactor core may be either inapplicable at all or of limited utility. Unfortunately, our formalism is not yet capable of making even qualitative predictions for the very small baselines ($L\lesssim L_0\sim \langle {\mathfrak{L}}_0\rangle$), where the LO ISLV correction is obviously insufficient. Thus, let us postpone this problem for later. In any case, in this article, we did not set ourselves the daunting task of studying all the competing models for the reactor antineutrino spectra. Instead, we take a pragmatic approach, using a representative set of the ${\overline{\nu }}_e$ energy spectra calculated by different methods, namely we use the spectra calculated by Huber [4], Mueller et al. [5], Fallot et al. [166], and Silaeva and Sinev [173], and compare these with the original conversions of the $\beta$-spectrum measurements performed with the spectrometers BILL (at the ILL in Grenoble) [176,177] and RONS (at the Rovno NPP) [80]. In addition, we combine theoretical or conversion-based spectra with the cumulative ${\overline{\nu }}_e$ spectrum from ${}^{238}$U fission obtained with the scientific neutron source FRM II in Garching [180], instead of the corresponding author’s contributions. Thus, we use ten different models in total.

Comparisons of the models under consideration are shown in Figure 5, Figure 6 and Figure 7, where the spectra are plotted as the ratios to the Mueller parametrization

$$S_k\left(E\right)=exp\left[{10}^{-3}\sum _{n=0}^5{a}_{kn}{\left(\frac{E}{1\mathrm{MeV}}\right)}^n\right]\frac{1}{\mathrm{fission}\mathrm{MeV}},$$

(14)

where $k={}^{235}\mathrm{U},{}^{238}\mathrm{U},{}^{239}\mathrm{Pu},{}^{241}\mathrm{Pu}$ and the coefficients $a_{kn}$ are given in the following matrix:

$$\begin{array}{ccccccc}n& 0& 1& 2& 3& 4& 5\\ {}^{235}\mathrm{U}& 3217& -3111& 13954& -369.0& 44.45& -2.053\\ {}^{238}\mathrm{U}& 483.3& 192.7& -128.3& -6.762& 2.233& -0.1536\\ {}^{239}\mathrm{Pu}& 6413& -7432& 3535& -882.0& 102.5& -4.550\\ {}^{241}\mathrm{Pu}& 3251& -3204& 1428& -367.5& 42.54& -1.896\end{array}$$

The corresponding spectra are plotted in Figure 8. To facilitate visual comparison of the models displayed in Figure 5, Figure 6 and Figure 7 with each other, we also show, as a reference model, the parametrization of the spectra calculated by Vogel and Engel [165] (of course also normalized to the Mueller parametrization). To simplify and speed up the subsequent calculations, we use the second degree B-spline interpolations of the spectra from all the models being discussed. Some examples of this interpolation are shown in Figure 5, Figure 6 and Figure 7. In case of the FRM II data, the spline is smoothly sewn with Mueller’s parameterization outside the measurement region. We have verified that the use, when available, of the author’s parameterization (similar to Equation (14)) has practically no effect on the results of the fits.

We emphasize that our goal is not at all to find out which of the spectrum models or method is the best according to one or another criterion, but it is to study the sensitivity of the parameters $N_0$ and $L_0$ to variations of the spectra. As can be seen from the figures, the spectra predicted in different models can differ by an amount that sometimes exceeds the declared uncertainties of the models. The recent Silaeva–Sinev calculation [173] demonstrates a spectral feature that resembles the notorious 5 MeV bump for the ${\overline{\nu }}_e$ fluxes from ${}^{238}$U and ${}^{241}$Pu and thus it can be used in particular for studying the interplay between the overall renormalization (dependent on the bump) and the ISLV effect.

In the present analysis, we do not take into account the correlated and uncorrelated uncertainties of the spectra predicted in each individual model, but instead we accumulate all the uncertainties in the single normalization factor $N_0$. We also neglect differences in the detection thresholds for the ${\overline{\nu }}_e$ energy, adopted (explicitly or implicitly) in different experiments and use integration from the IBD reaction threshold. Our estimates suggest that this simplification would only lead to marginal effects in our numerical analysis, when the experimental cuts are not too far from the kinematic boundaries. Similarly, it would be also an excess of accuracy to take into account the ${\beta }^{-}$ decay branches from ${}^{106}$Rh and ${}^{144}$Pr and antineutrinos from nonfuel materials in reactors [188], which contribute slightly above the IBD threshold. On the other hand, large discrepancies between the high-energy spectral shapes usually have a little effect on the outputs because of very small contributions of these tails. It should be, however, pointed out that the coefficients in the asymptotic expansion Equation (9) are functions of all kinematic variables of all particles and nuclei involved in the processes of (anti)neutrino production and detection, and, as a consequence, the effective parameter $L_0$ is in general sensitive to the detection threshold. Thus, the differences in the experimental detection thresholds and, more generally, in the experimental cuts add uncontrollable uncertainty to the global analysis. Considering that the expected ISLV effect itself is very small, all of these issues must be the subject of future investigations.

In addition, there are other sources of uncertainty that are potentially relevant to our analysis. One of them is the change in the IBD detection rate arising from evolution of the fuel content in the reactor, which may be especially important for the experiments based on commercial LEU reactors experiencing significant changes in fission rates during their fuel cycles (see, e.g., [189,190,191,192] for extended discussion). Therefore, it is important to study the robustness of the parameters $N_0$ and $L_0$, extracted from the fits, also to this effect. To simplify this job, we just compare the dependencies of the IBD yields calculated for several fuel compositions. Looking ahead, we note that the corresponding effect turns out to be either small or negligible at the current level of accuracy.

The stability of $L_0$ to these kinds of uncertainties is quite expected: the main influence of the ISLV effect is associated with short baselines, where the standard three-flavor oscillations are completely insignificant and nonstandard distortions (in particular, due to the hypothetical eV-scale sterile neutrinos) are, according to our assumption, either absent or small. As a result, the shapes of the ${\overline{\nu }}_e$ energy spectra and fuel composition are almost or fully inessential in the ratio Equation (16). In contrast, the normalization factor $N_0$ is equally sensitive to the data measured at all baselines. Moreover, the high-precision measurements from Daya Bay and RENO provide essential contributions to ${\chi }^2$. Since the two fitted parameters are strongly correlated, they both must be in general sensitive to the above-mentioned effects. The sensitivity to the subleading effects mentioned is relatively low, and our simplifications are adequate to current accuracy of the experimental data and theoretical inputs. However, as measurements and calculations improve (and if our hypothesis will not be disconfirmed), it will be important to use the most accurate up-to-date models of the antineutrino energy spectra and take into account all potentially significant sources of uncertainty, including bin-to-bin correlations.

6. Statistical Analysis

In our analysis, we use the data of three types which are described below.

6.1. The Main Dataset

The full set of data on the ratios, R, of the observed and predicted counting rates of the ${\overline{\nu }}_e$-induced IBD events and relevant quantities, required for the subsequent analyses, is presented in

Table A1. The main experimental dataset used in the subsequent analysis. The first column enumerates the experiments or (for segmented or movable detectors) baseline bins. The numbers in second column enumerate the groups of the cross-correlated measurements. The fractions of the main fissile isotopes are listed in columns 4–7. The ratios, R, of the observed to predicted integral production rates (used the analysis by Mention et al. [45] and corrected on updated neutron lifetime from [46]) are given in column 8. Columns 9–14 sequentially show the total experimental error, correlated error, baseline (or effective baseline), year of final publication or author’s communication, reference(s), and type of measurements (absolute or relative) for each experiment (bin). The overall uncertainty of the ${\overline{\nu }}_e$ flux is subtracted from the systematic errors of all absolute data. Some correlation errors have been revised. The two special cases—groups 4 and 13—are described in Section 6.4 and the sub-matrices ${\mathbf{V}}_4$ and ${\mathbf{V}}_{13}$ are given by Equations (22) and (23). The results of the experiment DANSS (group 17) are preliminary (the data are kindly provided by the DANSS collaboration [200]). The data for groups 16, 18, and 19 were obtained by digitizing the measured dependencies of the counting rates (presented in arbitrary units) on distance L, divided by the $\propto 1/L^2$ approximations of these data. The block-diagonal covariance matrix for the experiment groups 1–10, 12, 14, 15, 20, and 21 (experiments/bins # 1–26, 29, 38, 39, 153, and 154) is shown in Figure 9, and the matrices for groups 11 and 13 (relative data) are shown as inserts.

#	Gr.	Experiment	${}^{235}$U	${}^{239}$Pu	${}^{238}$U	${}^{241}$Pu	Ratio	${\mathit{\sigma }}^{\mathbf{exp}}$ (%)	${\mathit{\sigma }}_{\mathbf{ij}}^{\mathbf{corr}}$ (%)	L (m)	Year	Ref.	Type
1	1	ILL	1.0	0	0	0	0.796	9.01	4.80	8.76	1981	[7,8]	abs
2	1	Gösgen-I	0.620	0.274	0.074	0.042	0.959	5.36	4.80	37.9	1986	[16]	abs
3	1	Gösgen-II	0.584	0.298	0.068	0.050	0.985	5.33	4.80	45.9	1986	[16]	abs
4	1	Gösgen-III	0.543	0.329	0.070	0.058	0.918	6.80	4.80	64.7	1986	[16]	abs
5	2	Krasnoyarsk-I	1.0	0	0	0	0.929	6.00	4.84	32.8	1987	[17]	abs
6	2	Krasnoyarsk-II	1.0	0	0	0	0.946	20.40	4.76	92.3	1987	[17]	abs
7	2	Krasnoyarsk-III	1.0	0	0	0	0.940	4.22	3.40	57.0	1994	[19]	abs
8	2	Krasnoyarsk-IV	1.0	0	0	0	1.092	17.75	3.40	231.4	1994	[19]	abs
9	2	Krasnoyarsk-V	1.0	0	0	0	0.958	3.00	0	34.0	1999	[20]	abs
10	3	Rovno 88-1I	0.607	0.277	0.074	0.042	0.910	6.39	2.20	18.00	1988	[22]	abs
11	3	Rovno 88-2I	0.603	0.276	0.076	0.045	0.941	6.39	2.20	17.96	1988	[22]	abs
12	3	Rovno 88-1S	0.606	0.277	0.074	0.043	0.965	7.34	2.20	18.15	1988	[22]	abs
13	3	Rovno 88-2S	0.557	0.313	0.076	0.054	0.952	7.34	2.20	25.17	1988	[22]	abs
14	3	Rovno 88-3S	0.606	0.274	0.074	0.046	0.931	6.77	2.20	18.18	1988	[22]	abs
15	4	Rovno91	0.614	0.274	0.074	0.038	0.933	2.80		18.19	1991	[24]	abs
16	4	Rovno92	0.614	0.274	0.074	0.038	0.911	3.80	${\mathbf{V}}_4$	12.15	1992	[26]	abs
17	4	Bugey-4	0.538	0.328	0.078	0.056	0.935	1.40		14.88	1994	[29]	abs
18	5	Bugey-3-I	0.538	0.328	0.078	0.056	0.939	4.45	4.00	15	1995	[32]	abs
19	5	Bugey-3-II	0.538	0.328	0.078	0.056	0.945	4.54	4.00	40	1995	[32]	abs
20	5	Bugey-3-III	0.538	0.328	0.078	0.056	0.869	15.10	4.00	95	1995	[32]	abs
21	6	SRP-I	1.0	0	0	0	0.945	2.40	1.50	18.18	1996	[14]	abs
22	6	SRP-II	1.0	0	0	0	1.010	2.53	1.50	23.82	1996	[14]	abs
23	7	CHOOZ	0.496	0.351	0.087	0.066	0.954	3.42	–	1052	2003	[40]	abs
24	8	Palo Verde	0.600	0.270	0.070	0.060	0.968	5.42	–	833.2	2001	[36]	abs
25	9	Double Chooz	0.520	0.333	0.087	0.060	0.939	1.40	–	400.3	2020	[59]	abs
26	10	KamLAND	0.571	0.295	0.078	0.056	0.586	15.87	–	180$\times {10}^3$	2005	[41]	abs
27	11	RENO	0.556	0.326	0.071	0.047	0.999	1.78	1.50	409	2012	[47]	rel
28	11	RENO	0.556	0.326	0.071	0.047	0.930	1.92	1.50	1444	2012	[47]	rel
29	12	RENO	0.556	0.326	0.071	0.047	0.957	2.00	–	410.6	2020	[49]	abs
30	13	Daya Bay	0.564	0.303	0.076	0.056	0.995	0.70		517.78	2017	[48]	rel
31	13	Daya Bay	0.564	0.303	0.076	0.056	0.995	0.70		513.05	2017	[48]	rel
32	13	Daya Bay	0.564	0.303	0.076	0.056	0.994	0.63		554.39	2017	[48]	rel
33	13	Daya Bay	0.564	0.303	0.076	0.056	0.992	0.64	${\mathbf{V}}_{13}$	557.80	2017	[48]	rel
34	13	Daya Bay	0.564	0.303	0.076	0.056	0.938	0.68		1566.78	2017	[48]	rel
35	13	Daya Bay	0.564	0.303	0.076	0.056	0.941	0.68		1569.13	2017	[48]	rel
36	13	Daya Bay	0.564	0.303	0.076	0.056	0.942	0.68		1571.08	2017	[48]	rel
37	13	Daya Bay	0.564	0.303	0.076	0.056	0.948	0.69		1572.93	2017	[48]	rel
38	14	Daya Bay	0.564	0.304	0.076	0.056	0.944	2.83	–	578	2019	[50]	abs
39	15	Nucifer	0.926	0.061	0.008	0.005	1.024	10.65	–	7.21	2016	[56]	abs
40	16	Neutrino-4	1.0	0	0	0	0.859	4.83	0	7.000	2020	[88]	rel
41	16	Neutrino-4	1.0	0	0	0	0.936	6.05	0	8.400	2020	[88]	rel
42	16	Neutrino-4	1.0	0	0	0	0.801	9.17	0	9.840	2020	[88]	rel
43	16	Neutrino-4	1.0	0	0	0	0.861	8.56	0	11.230	2020	[88]	rel
44	16	Neutrino-4	1.0	0	0	0	0.755	13.24	0	6.41	2020	[88]	rel
45	16	Neutrino-4	1.0	0	0	0	0.717	14.23	0	6.64	2020	[88]	rel
46	16	Neutrino-4	1.0	0	0	0	0.912	10.42	0	6.88	2020	[88]	rel
47	16	Neutrino-4	1.0	0	0	0	0.700	14.29	0	7.11	2020	[88]	rel
48	16	Neutrino-4	1.0	0	0	0	0.822	12.77	0	7.35	2020	[88]	rel
49	16	Neutrino-4	1.0	0	0	0	0.913	13.14	0	7.58	2020	[88]	rel
50	16	Neutrino-4	1.0	0	0	0	0.913	12.05	0	7.82	2020	[88]	rel
51	16	Neutrino-4	1.0	0	0	0	0.900	12.78	0	8.06	2020	[88]	rel
52	16	Neutrino-4	1.0	0	0	0	0.938	16.52	0	8.29	2020	[88]	rel
53	16	Neutrino-4	1.0	0	0	0	0.812	20.32	0	8.53	2020	[88]	rel
54	16	Neutrino-4	1.0	0	0	0	0.849	20.02	0	8.76	2020	[88]	rel
55	16	Neutrino-4	1.0	0	0	0	0.681	26.43	0	9.00	2020	[88]	rel
56	16	Neutrino-4	1.0	0	0	0	0.645	27.91	0	9.23	2020	[88]	rel
57	16	Neutrino-4	1.0	0	0	0	0.672	29.02	0	9.47	2020	[88]	rel
58	16	Neutrino-4	1.0	0	0	0	0.551	36.30	0	9.70	2020	[88]	rel
59	16	Neutrino-4	1.0	0	0	0	0.780	26.28	0	9.94	2020	[88]	rel
60	16	Neutrino-4	1.0	0	0	0	0.756	22.49	0	10.17	2020	[88]	rel
61	16	Neutrino-4	1.0	0	0	0	0.973	16.44	0	10.41	2020	[88]	rel
62	16	Neutrino-4	1.0	0	0	0	0.886	20.88	0	10.64	2020	[88]	rel
63	16	Neutrino-4	1.0	0	0	0	0.352	62.50	0	10.88	2020	[88]	rel
64	16	Neutrino-4	1.0	0	0	0	1.072	24.25	0	11.12	2020	[88]	rel
65	16	Neutrino-4	1.0	0	0	0	0.727	39.89	0	11.35	2020	[88]	rel
66	16	Neutrino-4	1.0	0	0	0	0.950	38.95	0	11.58	2020	[88]	rel
67	16	Neutrino-4	1.0	0	0	0	1.246	32.91	0	11.82	2020	[88]	rel
68	17	DANSS	0.580	0.300	0.070	0.050	1.006	0.36	0	10.405	2020	[200]	rel
69	17	DANSS	0.580	0.300	0.070	0.050	1.012	0.43	0	11.395	2020	[200]	rel
70	17	DANSS	0.580	0.300	0.070	0.050	1.017	0.50	0	12.385	2020	[200]	rel
71	17	DANSS	0.580	0.300	0.070	0.050	1.010	0.38	0	10.710	2020	[200]	rel
72	17	DANSS	0.580	0.300	0.070	0.050	1.013	0.45	0	11.700	2020	[200]	rel
73	17	DANSS	0.580	0.300	0.070	0.050	1.009	0.44	0	12.690	2020	[200]	rel
74	17	DANSS	0.580	0.300	0.070	0.050	1.008	0.40	0	11.015	2020	[200]	rel
75	17	DANSS	0.580	0.300	0.070	0.050	1.021	0.47	0	12.005	2020	[200]	rel
76	17	DANSS	0.580	0.300	0.070	0.050	1.006	0.56	0	12.995	2020	[200]	rel
77	18	PROSPECT	0.990	0.003	0.003	0.003	0.994	5.92	0	6.955	2021	[53]	rel
78	18	PROSPECT	0.990	0.003	0.003	0.003	0.865	6.02	0	7.064	2021	[53]	rel
79	18	PROSPECT	0.990	0.003	0.003	0.003	0.978	6.11	0	7.158	2021	[53]	rel
80	18	PROSPECT	0.990	0.003	0.003	0.003	1.044	5.58	0	7.246	2021	[53]	rel
81	18	PROSPECT	0.990	0.003	0.003	0.003	1.003	5.10	0	7.346	2021	[53]	rel
82	18	PROSPECT	0.990	0.003	0.003	0.003	0.923	4.69	0	7.351	2021	[53]	rel
83	18	PROSPECT	0.990	0.003	0.003	0.003	0.933	5.77	0	7.361	2021	[53]	rel
84	18	PROSPECT	0.990	0.003	0.003	0.003	1.037	6.39	0	7.399	2021	[53]	rel
85	18	PROSPECT	0.990	0.003	0.003	0.003	0.947	5.10	0	7.448	2021	[53]	rel
86	18	PROSPECT	0.990	0.003	0.003	0.003	0.957	5.51	0	7.450	2021	[53]	rel
87	18	PROSPECT	0.990	0.003	0.003	0.003	0.964	5.72	0	7.457	2021	[53]	rel
88	18	PROSPECT	0.990	0.003	0.003	0.003	1.063	6.71	0	7.489	2021	[53]	rel
89	18	PROSPECT	0.990	0.003	0.003	0.003	1.100	7.58	0	7.513	2021	[53]	rel
90	18	PROSPECT	0.990	0.003	0.003	0.003	0.899	4.70	0	7.551	2021	[53]	rel
91	18	PROSPECT	0.990	0.003	0.003	0.003	0.954	5.26	0	7.652	2021	[53]	rel
92	18	PROSPECT	0.990	0.003	0.003	0.003	1.010	6.26	0	7.656	2021	[53]	rel
93	18	PROSPECT	0.990	0.003	0.003	0.003	0.919	5.01	0	7.661	2021	[53]	rel
94	18	PROSPECT	0.990	0.003	0.003	0.003	1.027	5.27	0	7.673	2021	[53]	rel
95	18	PROSPECT	0.990	0.003	0.003	0.003	0.958	6.54	0	7.715	2021	[53]	rel
96	18	PROSPECT	0.990	0.003	0.003	0.003	1.059	9.61	0	7.744	2021	[53]	rel
97	18	PROSPECT	0.990	0.003	0.003	0.003	0.925	6.50	0	7.754	2021	[53]	rel
98	18	PROSPECT	0.990	0.003	0.003	0.003	0.939	6.72	0	7.755	2021	[53]	rel
99	18	PROSPECT	0.990	0.003	0.003	0.003	0.885	5.48	0	7.758	2021	[53]	rel
100	18	PROSPECT	0.990	0.003	0.003	0.003	0.942	4.47	0	7.768	2021	[53]	rel
101	18	PROSPECT	0.990	0.003	0.003	0.003	0.902	5.43	0	7.783	2021	[53]	rel
102	18	PROSPECT	0.990	0.003	0.003	0.003	1.021	6.66	0	7.804	2021	[53]	rel
103	18	PROSPECT	0.990	0.003	0.003	0.003	0.934	4.55	0	7.865	2021	[53]	rel
104	18	PROSPECT	0.990	0.003	0.003	0.003	1.004	4.60	0	7.877	2021	[53]	rel
105	18	PROSPECT	0.990	0.003	0.003	0.003	1.089	7.29	0	7.946	2021	[53]	rel
106	18	PROSPECT	0.990	0.003	0.003	0.003	0.984	6.29	0	7.958	2021	[53]	rel
107	18	PROSPECT	0.990	0.003	0.003	0.003	1.036	8.01	0	7.959	2021	[53]	rel
108	18	PROSPECT	0.990	0.003	0.003	0.003	0.950	4.84	0	7.962	2021	[53]	rel
109	18	PROSPECT	0.990	0.003	0.003	0.003	0.909	4.35	0	7.972	2021	[53]	rel
110	18	PROSPECT	0.990	0.003	0.003	0.003	1.052	5.14	0	7.987	2021	[53]	rel
111	18	PROSPECT	0.990	0.003	0.003	0.003	0.887	5.75	0	8.007	2021	[53]	rel
112	18	PROSPECT	0.990	0.003	0.003	0.003	0.915	6.92	0	8.032	2021	[53]	rel
113	18	PROSPECT	0.990	0.003	0.003	0.003	0.951	5.02	0	8.062	2021	[53]	rel
114	18	PROSPECT	0.990	0.003	0.003	0.003	1.055	8.56	0	8.063	2021	[53]	rel
115	18	PROSPECT	0.990	0.003	0.003	0.003	0.923	4.49	0	8.068	2021	[53]	rel
116	18	PROSPECT	0.990	0.003	0.003	0.003	0.996	4.69	0	8.081	2021	[53]	rel
117	18	PROSPECT	0.990	0.003	0.003	0.003	0.937	5.67	0	8.098	2021	[53]	rel
118	18	PROSPECT	0.990	0.003	0.003	0.003	0.897	6.16	0	8.120	2021	[53]	rel
119	18	PROSPECT	0.990	0.003	0.003	0.003	1.061	6.00	0	8.162	2021	[53]	rel
120	18	PROSPECT	0.990	0.003	0.003	0.003	0.876	4.92	0	8.167	2021	[53]	rel
121	18	PROSPECT	0.990	0.003	0.003	0.003	0.959	4.77	0	8.176	2021	[53]	rel
122	18	PROSPECT	0.990	0.003	0.003	0.003	0.909	9.16	0	8.181	2021	[53]	rel
123	18	PROSPECT	0.990	0.003	0.003	0.003	0.921	5.35	0	8.190	2021	[53]	rel
124	18	PROSPECT	0.990	0.003	0.003	0.003	1.005	6.43	0	8.235	2021	[53]	rel
125	18	PROSPECT	0.990	0.003	0.003	0.003	0.921	8.46	0	8.265	2021	[53]	rel
126	18	PROSPECT	0.990	0.003	0.003	0.003	0.921	6.49	0	8.266	2021	[53]	rel
127	18	PROSPECT	0.990	0.003	0.003	0.003	0.976	5.85	0	8.273	2021	[53]	rel
128	18	PROSPECT	0.990	0.003	0.003	0.003	0.913	6.00	0	8.284	2021	[53]	rel
129	18	PROSPECT	0.990	0.003	0.003	0.003	0.981	7.38	0	8.301	2021	[53]	rel
130	18	PROSPECT	0.990	0.003	0.003	0.003	1.086	6.85	0	8.323	2021	[53]	rel
131	18	PROSPECT	0.990	0.003	0.003	0.003	1.024	7.81	0	8.350	2021	[53]	rel
132	18	PROSPECT	0.990	0.003	0.003	0.003	0.817	6.39	0	8.380	2021	[53]	rel
133	18	PROSPECT	0.990	0.003	0.003	0.003	1.007	12.88	0	8.382	2021	[53]	rel
134	18	PROSPECT	0.990	0.003	0.003	0.003	1.050	6.93	0	8.437	2021	[53]	rel
135	18	PROSPECT	0.990	0.003	0.003	0.003	0.923	6.80	0	8.477	2021	[53]	rel
136	18	PROSPECT	0.990	0.003	0.003	0.003	0.980	7.95	0	8.505	2021	[53]	rel
137	18	PROSPECT	0.990	0.003	0.003	0.003	1.025	6.70	0	8.526	2021	[53]	rel
138	18	PROSPECT	0.990	0.003	0.003	0.003	0.890	10.80	0	8.552	2021	[53]	rel
139	18	PROSPECT	0.990	0.003	0.003	0.003	1.087	8.92	0	8.584	2021	[53]	rel
140	18	PROSPECT	0.990	0.003	0.003	0.003	0.796	6.69	0	8.616	2021	[53]	rel
141	18	PROSPECT	0.990	0.003	0.003	0.003	0.965	9.79	0	8.640	2021	[53]	rel
142	18	PROSPECT	0.990	0.003	0.003	0.003	0.983	7.11	0	8.708	2021	[53]	rel
143	18	PROSPECT	0.990	0.003	0.003	0.003	0.971	9.46	0	8.729	2021	[53]	rel
144	18	PROSPECT	0.990	0.003	0.003	0.003	0.855	8.10	0	8.802	2021	[53]	rel
145	18	PROSPECT	0.990	0.003	0.003	0.003	1.075	8.38	0	8.820	2021	[53]	rel
146	18	PROSPECT	0.990	0.003	0.003	0.003	1.050	9.87	0	8.912	2021	[53]	rel
147	19	STEREO	1.0	0	0	0	1.015	2.13	0	9.42	2020	[51]	rel
148	19	STEREO	1.0	0	0	0	1.011	2.17	0	9.77	2020	[51]	rel
149	19	STEREO	1.0	0	0	0	0.993	2.37	0	10.13	2020	[51]	rel
150	19	STEREO	1.0	0	0	0	1.018	2.48	0	10.48	2020	[51]	rel
151	19	STEREO	1.0	0	0	0	0.978	2.64	0	10.84	2020	[51]	rel
152	19	STEREO	1.0	0	0	0	1.031	3.02	0	11.18	2020	[51]	rel
153	20	STEREO	1.0	0	0	0	0.954	3.53	–	10.3	2020	[57]	abs
154	21	MiniCHANDLER	0.580	0.300	0.070	0.050	0.824	18.20	–	25	2020	[58]	abs

(see Appendix A). Content of the table is detailed in the caption, but some additional comments need to be made.

The predicted rates are calculated using the results of Mention et al. [45] based on the reactor ${\overline{\nu }}_e$ flux prediction by Mueller et al. [5]. Essentially, all reactor experiments used Vogel’s IBD cross section [193], which is inversely proportional to the neutron mean life, ${\tau }_n$, and the latter is subject to continuous refinement. Thus, we systematically rescaled the ratios to account for the latest value, ${\tau }_n=879.4\pm 0.6$ s, recommended by the Particle Data Group [46]. Typically, this reduces the ratios within ≲1%; the specific value of the shift depends on the particular value of ${\tau }_n$ used in the related experiment and/or in the calculations of Mention et al. [45]. Conversion of the data to another spectrum model also reduces to a proper renormalization of R defined by the ratio of the theoretical rates calculated with the two models.

The diagonal and correlation errors shown in columns 9 and 10 of the table do not fully determine the correlations in the respective entries of the covariance matrix.

The latter is shown in Figure 9 for experiment groups 1–10, 12, 14, 15, 20, and 21 (experiments/bins # 1–26, 29, 38, 39, 153, and 154). Although it is constructed for the Mueller’s flux, dependency on the flux model is insignificant. All the data are formally correlated through the uncertain ${\overline{\nu }}_e$ flux normalization. This uncertainty, common to all absolute measurements, is removed from the systematic errors in order to use this poorly known quantity as a variable parameter in test calculations. Consequently, it is also subtracted from the covariance matrix shown in Figure 9.

The data from certain earlier experiments listed in Table A1 are sometimes quite different from those presented in similar compilations (see, e.g., [70,72,76,194,195]). In particular, in the recent study [195], the ratios for the most of the earlier data were reanalysed using the IBD rates calculated with the Huber–Mueller [4,5], Estienne et al. (ab initio) [172], and Hayen et al. [182] ${\overline{\nu }}_e$ flux predictions. Comparison with these results indicates sometimes significant differences from the results of Mention et al. [45] based on Mueller’s spectrum. We intend to include these results in future analysis, when new and improved current reactor data will become available.

For the experiments, Double Chooz, Daya Bay, RENO, and KamLAND, which operate with the ${\overline{\nu }}_e$ fluxes from several or many reactors, we use the effective (flux-weighted) baselines, as provided by the authors. For example, in the case of KamLAND, ∼80% of the total ${\overline{\nu }}_e$ events come from commercial reactors at $175\pm 35$ km from the detector site [60], more distant reactors increase the effective baseline to about 180 km. As another example, we remark that the flux-weighted baselines for the RENO experiment, given in Table A1 according to [47], differ slightly from newer estimates: $410.6$ m for the near detector and $1445.7$ m for the far detector [197]. Fortunately, this difference has no sizable effect on the global fits. Let us mention another uncertainty related to RENO. We did not find definite information on correlation between points # 27 and 28. We, however, verified that the fits are almost insensitive to this correlation. For definiteness, we use ${\sigma }_{\mathrm{corr}}^{27,28}=1.5\%$—the value typical for the RENO measurements prior to 2013. The unpublished covariance matrix for the Daya Bay data (insert in Figure 9) is courtesy of the Daya Bay collaboration [196].

The zeros in the tenth column (${\sigma }_{\mathrm{corr}}^{ij}$) for groups 16–19 (Neutrino-4, DANSS, PROSPECT, STEREO) do not in any way mean that correlations are absent or small, rather they indicate that the information about correlations is unavailable, and thus we are forced to ignore them. However, these datasets definitely do not correlate with each other and with other datasets and therefore the normalization factors can be calculated individually for each of them (and similarly for groups 11 and 13).

The Neutrino-4 data # 40–43 and 44–67 [88,89] are not from independent measurements, but are two representations of the same data sample for, respectively, wide and narrow baseline binnings. In our analysis, we adopt the datasets # 40–43 as more statistically valid, while the datasets # 44–67 are used for test calculations (see below). Let us mention in passing that, in the earlier eprint versions of [88,89], another binning was used (three data points). We verified that the two representations of the Neutrino-4 data sample (with three and four points) yield the same (within the specified precision) result in our analysis and so we will no longer dwell on this issue. Recall that we do not take up the intriguing results of the spectral Neutrino-4 measurements, (see [99,100,101,102,103] for the comments and authors’ responses) arbitrarily assuming that the integral event rate is not related to the sterile neutrinos.

The high-precision absolute measurements of RENO (# 29) [47], Daya Bay (# 38) [50,198], and STEREO (# 153) [57] are the results of more advanced (compared to the pertinent relative measurements) data processing and use additional averaging over the detectors (RENO, Daya Bay) or detector cells (STEREO). Thus, the corresponding values of L given in the table are effective baselines (e.g., for the Daya Bay measurement # 38, the latter is the flux-weighted baseline for the two near halls). Although the effective baseline is slightly dependent on the spectrum model and the dataset used (compare, e.g., the data on the absolute rate and baseline reported in [50,199]), this has no visible effect on the results of our analysis. In order to avoid misunderstandings, it should also be noted that the reactor antineutrino yields reported in [50,198,199] were normalized to the Huber–Mueller theoretical prediction, and the resulting ratio was corrected for the expected oscillation effect. In Table A1, the oscillation effect has been removed back from the Daya Bay absolute data (the oscillations are of course included in our analysis). In our calculations, we regard the data from groups 11 and 12 (RENO) and 13 and 14 (Daya Bay) to be independent and decrease the full NDF by 1 for each group of the relative data. Even though this approach is rather naive, we have no reliable way of assessing the correlation between the relative and absolute data obtained in the same experiment. On the other hand, both types of the data are equally significant: the first is responsible for the form of the function $R=R\left(L\right)$, and the second is for the normalization.

The data from the DANSS experiment at Kalinin NPP (group 17 in Table A1) are of particular interest because of their unprecedented accuracy due to the huge count rate, clear IBD signature, reliable shielding, and good control of systematics. These data are not final and have not yet been published but were reported at several meetings, e.g., at the TAUP’2019 conference [54] and at ICHEP’2020 virtual conference [55]. To obtain this sample, the detector fiducial volume was divided into three vertical sections, and, to account for the poorly known efficiencies of these, the measured rates were renormalized to the top (closest to reactor) section. The simplest model of the reactor with a uniform core was used for fitting. Currently, the statistics in DANSS are still being collected, the selection criteria and analysis procedure are being clarified, and backgrounds are being refined; the final results are expected soon [200].

The data point # 154 from MiniCHANDLER [58] (a 80 kg prototype of the full detector CHANDLER [201]) is currently of methodological interest, demonstrating the potential of the novel technology.

6.2. “Position l to Position 2 Ratios”

Another type of relative data are given by measurements of the IBD count rates ($N_{1,2}$) at two distances ($L_{1,2}$) from the reactor core, made with one movable detector or with two identical detectors. It is clear that the solid-angle corrected “position l to position 2” ratio

$$R_{12}=\frac{N_1}{N_2}{\left(\frac{L_2}{L_1}\right)}^2$$

(15)

is less dependent (compared to the absolute ratios R) on knowledge of the initial spectrum, cross sections, and some other systematic effects because most of the uncertainties cancel out. Assuming classical ISL, the ratio Equation (15) should be equal to 1 in a perfect experiment.

Table A2. Solid-angle corrected ratios ($R_{12}$) of the IBD events measured at two distances ($L_1$, $L_2$) from reactor core by using movable or identical detectors. Some data are recalculated from the originally published inverse ratios ($R_{21}$), event number or cross section ratios. Subscript “ det+react” means that the error includes statistics, uncertainty associated with the detection efficiency, and insufficient knowledge of the reactor power. Data marked with tick (√) are used in the global statistical analysis. The rest data are not involved because they are either not independent, or have already been included in a different form, or are out of date. All the data, except those marked with an asterisk (*), are shown in Figure 20. The earlier measurements with high ${\overline{\nu }}_e$ energy thresholds (e.g. [10]) are not included. We show, for contrast, obviously outdated data from experiments with the Irvine mobile neutrino oscillation detector in the SRP Neutrino Laboratory [12,13] and with the movable Bugey detector [27].

#	Experiment	Year	Ref.	${\mathit{L}}_{\mathbf{1}}$ (m)	${\mathit{L}}_{\mathbf{2}}$ (m)	${\mathit{R}}_{\mathbf{12}}$	Comment
1	Gösgen	1984	[15]	$37.9$	$45.9$	$0.99\pm {0.03}_{}\mathrm{stat}\pm {0.02}_{\mathrm{syst}}$	from inverse ratio	√
2	Krasnoyarsk	1987	[17]	$32.8$	$92.3$	$0.984\pm {0.168}_{}\mathrm{stat}+\mathrm{syst}$	from event number ratio
						$0.984\pm {0.209}_{}\mathrm{stat}+\mathrm{syst}$	from cross section ratio	*
3	Krasnoyarsk	1990	[18]	$32.8$	$57.0$	$1.01\pm {0.07}_{}\det +\mathrm{react}$	from cross section ratio
4				$32.8$	$57.6$	$1.02\pm {0.07}_{}\det +\mathrm{react}$	from cross section ratio
5				$57.0$	$57.6$	$1.01\pm {0.03}_{}\mathrm{stat}+\mathrm{syst}$	from event number ratio	√
						$1.01\pm {0.06}_{}\det +\mathrm{react}$	from cross section ratio	*
6				$57.6$	$231.4$	$0.86\pm {0.22}_{}\mathrm{stat}+\mathrm{syst}$	from event number ratio
7				$57.0$	$231.4$	$0.87\pm {0.22}_{}\mathrm{stat}+\mathrm{syst}$	from event number ratio
8	Krasnoyarsk	1994	[19]	$57.0$	$231.4$	$0.86\pm {0.15}_{}\mathrm{stat}+\mathrm{syst}$	from event number ratio
9	Rovno	1986	[21]	18	25	$0.986\pm {0.040}_{}\mathrm{stat}\pm {0.029}_{}\mathrm{syst}$
10	Rovno	1988	[22]	18	25	$0.987\pm {0.039}_{}\mathrm{stat}+\mathrm{syst}$	from inverse ratio	√
11	Rovno	1991	[24]	$325.66$	$336.38$	$0.993\pm {0.042}_{}\mathrm{stat}+\mathrm{syst}$	from cross section ratio	√
12	Rovno	1992	[26]	$12.15$	$18.34$	$0.976\pm {0.020}_{}\mathrm{stat}\pm {0.015}_{}\mathrm{syst}$		√
13	SRP	1986	[12]	$18.2$	$23.7$	$1.038\pm {0.0137}_{}\mathrm{stat}$
14	SRP	1987	[13]	$18.2$	$23.8$	$1.034\pm {0.010}_{}\mathrm{stat}\pm {0.020}_{}\mathrm{syst}$
15	SRP	1996	[14]	$18.18$	$23.82$	$0.938\pm {0.012}_{}\mathrm{stat}\pm {0.017}_{}\mathrm{syst}$	from event number ratio
16	Bugey	1984	[27]	$13.63$	$18.30$	$1.102\pm {0.014}_{}\mathrm{stat}\pm {0.026}_{}\mathrm{syst}$
17	Bugey	1988	[30]	$13.63$	$18.30$	$0.997\pm {0.020}_{}\mathrm{stat}\pm {0.021}_{}\mathrm{syst}$	integrated spectra
18				$13.63$	$18.30$	$1.008\pm {0.039}_{}\mathrm{stat}\pm {0.021}_{}\mathrm{syst}$	standard shielding data
19				$13.63$	$18.30$	$1.007\pm {0.021}_{}\mathrm{stat}\pm {0.021}_{}\mathrm{syst}$	data with lowest threshold
20	Bugey	1995	[32]	15	40	$1.004\pm {0.015}_{}\mathrm{stat}$	from inverse spectrum ratio	√
						$0.994\pm {0.011}_{}\mathrm{stat}$	from double ratio	*
21				15	95	$0.997\pm {0.146}_{}\mathrm{stat}$	from inverse spectrum ratio	√
						$1.080\pm {0.156}_{}\mathrm{stat}$	from double ratio	*

collects available to our knowledge data on this quantity. Some of the data are recalculated from the originally published inverse ratios ($R_{21}$), event number, or ratios of measured cross section. Selected experimental details are given in the penultimate column. We do not show outdated results that have been revised in subsequent analyses (see, e.g., [12,13,27,28] for the earlier SRP and Bugey measurements). Only a small part of the data from Table A2 is used in the subsequent fits.

6.3. Theoretical Model for Analysis

To check the assumption that the ISLV effect could actually be, at least in part, responsible for the RAA, we perform a statistical analysis of the earlier and new reactor data. Since in this paper, we use only the spectrum-averaged event rates, the distance-dependent factor in Equation (13) can be replaced by the simple common factor $(1-L_0^2/L^2)$, in which $L_0\sim \langle {\mathfrak{L}}_0\rangle \ll L$ is an energy independent effective parameter, which characterizes the ISLV spatial scale and which is the main subject of the present statistical analysis. Recall that the function ${\mathfrak{L}}_0$ can be, strictly speaking, different for different reactor/detector pairs. We, however, neglect possible differences and apply the mean value theorem to reduce the problem to the single length-dimension parameter, common to all reactor experiments; it is clear that this simplification can only be justified a posteriori.

It is implicitly assumed that the dimensions of the reactor core and detector (or detector section) are small (although not necessarily very small) compared to the distance between them. To check the applicability of this simplification to averaging over the source/detector volumes, let us consider a simplified example: neglecting the dimensions of the detector, we take into account the dimensions of the source (reactor core). Let the latter be a uniformly fissile fuel filled cylinder with height h and radius r, the end of which is directed towards the coaxially placed detector. Then, the ratio, $\zeta$, of the exactly volume-averaged counting rate and our approximation can be estimated as

$$\zeta =\underset{0}{\overset{r}{\int }}\underset{-h/2}{\overset{h/2}{\int }}\frac{\rho d\rho dz}{{(L-z)}^2+{\rho }^2}\left[1-\frac{L_0^2}{{(L-z)}^2+{\rho }^2}\right]{\left\{\left(1-\frac{L_0^2}{L^2}\right)\left[\underset{0}{\overset{r}{\int }}\underset{-h/2}{\overset{h/2}{\int }}\frac{\rho d\rho dz}{{(L-z)}^2+{\rho }^2}\right]\right\}}^{-1},$$

where L is the distance between the source and detector centers, which is assumed to be short enough to neglect the oscillation effect, but large compared to r and h. Then, expanding $\zeta$ in powers of $\omega =r/L$ and $\eta =h/L$, we obtain

$$\zeta =1+\frac{1}{2}{\left(\frac{L^2}{L_0^2}-1\right)}^{-1}\left({\omega }^2-\frac{7}{6}{\eta }^2-\frac{5}{6}{\omega }^4+\frac{11}{6}{\omega }^2{\eta }^2-\frac{11}{24}{\eta }^4+\dots \right).$$

Plainly, $|1-\zeta |\ll 1$ owing to a partial compensation of contributions from longer and shorter distances. For a numerical example, let us consider the DANSS experiment operating a big industrial reactor with $r=1.56$ m and $h=3.55$ m (cf. $0.2\mathrm{m}r\times 0.5\mathrm{m}h$ for the HFIR reactor core in the PROSPECT experiment). The distance from the center of the reactor core to the nearest position of the top section of the (relatively) compact neutrino spectrometer is about $11.4$ m, therefore, taking $L_0=3$ m as a conservative upper limit (see below), we see that $0<1-\zeta \lesssim 0.005$; and $1-\zeta \lesssim 0.002$ as $L_0<2$ m. Such a correction (although easily accounted for) can be safely ignored at the current level of the data accuracy. Thus, we use the following theoretical model to fit the reactor data:

$$R^{\mathrm{th}}(L;L_0)=\left(1-\frac{L_0^2}{L^2}\right){\displaystyle \frac{\underset{0}{\overset{\infty }{\int }}d{E}_{\nu }\sum _k{f}_k{\mathcal{P}}_{\mathrm{surv}}^{3\nu }(L,E_{\nu }){\sigma }_{\mathrm{IBD}}\left(E_{\nu }\right)S_k\left(E_{\nu }\right)}{\underset{0}{\overset{\infty }{\int }}d{E}_{\nu }\sum _k{f}_k{\sigma }_{\mathrm{IBD}}\left(E_{\nu }\right)S_k\left(E_{\nu }\right)}}.$$

(16)

Here, $f_k$ is the reactor fissile isotope fraction, $S_k\left(E_{\nu }\right)$ is the ${\overline{\nu }}_e$ energy spectrum, ${\sigma }_{\mathrm{IBD}}\left(E_{\nu }\right)$ is the IBD cross section, and ${\mathcal{P}}_{\mathrm{surv}}^{3\nu }(L,E_{\nu })\equiv {\mathcal{P}}_{ee}^{\mathrm{QM}}(L,E_{\nu })$ is the ${\overline{\nu }}_e$ survival probability in the standard $3\nu$ mixing scheme:

$$\begin{array}{cc}\hfill {\mathcal{P}}_{\mathrm{surv}}^{3\nu }(L,E_{\nu })=& 1-{sin}^2\left(2{\theta }_{13}\right)\left({cos}^2{\theta }_{12}{sin}^2{\Delta }_{31}+{sin}^2{\theta }_{12}{sin}^2{\Delta }_{32}\right)\hfill \\ \hfill & -{cos}^4{\theta }_{13}{sin}^2\left(2{\theta }_{12}\right){sin}^2{\Delta }_{21},\hfill \end{array}$$

(17)

where ${\Delta }_{ij}=\Delta m_{ij}^{2}L/\left(4E_{\nu }\right)$ and, as the current reactor data are not sensitive to the $CP$-violating phase and neutrino mass hierarchy, the normal hierarchy is assumed (thus $\Delta m_{32}^2=\Delta m_{31}^2-\Delta m_{21}^2$). In our calculations, we use the central (best-fit) values for the neutrino mixing angles and mass-squared splittings obtained in a comprehensive up-to-date global analysis (“NuFIT 5.0”) of the three-flavor neutrino oscillations performed by Esteban et al. (NuFIT team) [43]:

$$\begin{array}{c}{\theta }_{12}={\left(33.{44}_{-0.74}^{+0.77}\right)}^{\circ },{\theta }_{13}={\left(8.{57}_{-0.12}^{+0.12}\right)}^{\circ },\\ \Delta m_{21}^2=\left(7.{42}_{-0.20}^{+0.21}\right)\times {10}^{-5}{\mathrm{eV}}^2,\Delta m_{31}^2=\left(2.{517}_{-0.028}^{+0.026}\right)\times {10}^{-3}{\mathrm{eV}}^2.\end{array}$$

(18)

The fit involved essentially all available data of solar, atmospheric, reactor, and accelerator neutrino oscillation experiments. We have verified that the results of our fits of the reactor data are very weakly sensitive to variations of the parameters Equation (18) within $1\sigma$ ranges. The NuFIT 5.0 results are in good agreement with the results of recent independent global-fit analyses by Capozzi et al. [202] (an update of [203]) and by de Salas et al. [204], which, along with the oscillation data, include other relevant neutrino physics inputs. Consequently, there is no point in including the oscillation parameters as fitting variables into the statistical analysis, at least at today’s level of the reactor data accuracy.

The “3 + 1” curves shown in Figure 2 and Figure 3 are calculated as

$$R^{3+1}\left(L\right)={\displaystyle \frac{\underset{0}{\overset{\infty }{\int }}d{E}_{\nu }\sum _k{f}_k{\mathcal{P}}_{\mathrm{surv}}^{4\nu }(L,E_{\nu }){\sigma }_{\mathrm{IBD}}\left(E_{\nu }\right)S_k\left(E_{\nu }\right)}{\underset{0}{\overset{\infty }{\int }}d{E}_{\nu }\sum _k{f}_k{\sigma }_{\mathrm{IBD}}\left(E_{\nu }\right)S_k\left(E_{\nu }\right)}},$$

(19)

where the ${\overline{\nu }}_e$ survival probability is given by

$$P_{\mathrm{surv}}^{4\nu }=1-{cos}^4{\theta }_{14}\left(1-P_{\mathrm{surv}}^{3\nu }\right)-{sin}^2\left(2{\theta }_{14}\right){sin}^2{\tilde{\Delta }}_{41}$$

(20)

with

$${sin}^2{\tilde{\Delta }}_{41}={cos}^2{\theta }_{13}\left({cos}^2{\theta }_{12}{sin}^2{\Delta }_{41}+{sin}^2{\theta }_{12}{sin}^2{\Delta }_{42}\right)+{sin}^2{\theta }_{13}{sin}^2{\Delta }_{43}\stackrel{\mathrm{VSBL}}{\simeq }{sin}^2{\Delta }_{41},$$

${\Delta }_{42}={\Delta }_{41}-{\Delta }_{21}$, ${\Delta }_{43}={\Delta }_{41}-{\Delta }_{31}$; the last equality assumes that $m_4^2\gg {max}_k\left(m_k^2\right)$ ($k=1,2,3$). At very short baselines, Equation (20) is simplified to $P_{\mathrm{surv}}^{4\nu }\simeq 1-{sin}^2\left(2{\theta }_{14}\right){sin}^2{\Delta }_{41}$.

The IBD cross section, ${\sigma }_{\mathrm{IBD}}\left(E_{\nu }\right)$, is calculated by using the detailed analytical results of Ivanov et al. [44], which take into account

  (i)

the contributions of weak magnetism and the neutron recoil to next-to-leading order in the large nucleon mass ($M=(m_p+m_n)/2$) expansion,

 (ii)

the radiative corrections to order $\alpha /\pi$, caused by single virtual photon exchanges and the radiative IBD, calculated to leading order in the $1/M$ expansion,

(iii)

the radiative corrections defined by W- and Z-boson exchanges and QCD corrections, taken as a common factor $1+{\Delta }_R$ [205] (see [44] for details and references).

The uncertainties in the multiplicative corrections are naturally absorbed in the flux normalization factors.

6.4. ${\chi }^2$ Scheme, Covariances

Taking into account the large uncertainty in the ${\overline{\nu }}_e$ flux normalization, we will consider the overall normalization factor $N_0$ to be a fitting parameter. To find the best-fit parameters $N_0$ and $L_0$, we minimize the standard ${\chi }^2$ for the correlated and uncorrelated data:

$${\chi }^2={\chi }_{\mathrm{abs}}^2+{\chi }_{\mathrm{rel}}^2+{\chi }_{12}^2.$$

(21)

Here,

$${\chi }_{\mathrm{abs}}^2={\left({\mathit{R}}^{\exp }-N_0{\mathit{R}}^{\mathrm{th}}\right)}^T{\mathit{W}}^{-1}\left({\mathit{R}}^{\exp }-N_0{\mathit{R}}^{\mathrm{th}}\right),$$

${\mathit{R}}^{\exp }$ is the vector constructed from the experimental ratios ($R_i^{\exp }$) marked in Table A1 as absolute measurements, ${\mathit{R}}^{\mathrm{th}}$ is the vector of the relevant theoretical predictions Equation (16) ($R_i^{\mathrm{th}}=R_i^{\mathrm{th}}(L_i;L_0)$),

$$\mathit{W}=∥Cov(R_i^{\exp },R_j^{\exp })+{\delta }_{\mathrm{flux}}^2{R}_i^{\mathrm{th}}{R}_j^{\mathrm{th}}∥$$

is the full covariance matrix shown (with ${\delta }_{\mathrm{flux}}=0$) in Figure 9, where ${\delta }_{\mathrm{flux}}$ is the relative uncertainty in the theoretical flux prediction, presumably common to all absolute measurements. The actual value of ${\delta }_{\mathrm{flux}}$ is not well known, but it can be shown that $N_0$ extracted in a one-parameter fit with fixed $L_0$ and both $N_0$ and $L_0$ extracted from a two-parameter fit are insensitive to ${\delta }_{\mathrm{flux}}$. The latter, however, affects the error in the determination of $N_0$ and hence the confidence contours in the $(N_0,L_0)$ plane (see below). It also slightly affects the value of $L_0$ in a one-parameter fit with fixed $N_0$ (for ${\delta }_{\mathrm{flux}}\lesssim 5\%$), but such a fit is not very interesting and we do not discuss it in what follows. For definiteness, we use ${\delta }_{\mathrm{flux}}=2.7\%$ everywhere, as suggested in [45]; this value has been subtracted from the systematic errors of the data on R. It is needless to say that the error in determining $N_0$ automatically absorbs constant or weakly energy-dependent uncertainties of various kinds (like the uncertainty in the neutron lifetime or factorizable corrections to the IBD cross section) common to all experiments.

A few technical details about the covariance matrix are needed to be said. We define ${\sigma }_{ij}^{\mathrm{corr}}$ to be the systematic error of i-th measurement correlated with j-th measurement. It is convenient to define the matrix $||{\upsilon }_{ij}||$ with ${\upsilon }_{ij}={\sigma }_{ij}^{\mathrm{corr}}$ for $i\ne j$ and ${\upsilon }_{ii}={\sigma }^{\exp }$. Then,

$$Cov(R_i^{\exp },R_j^{\exp })={10}^{-4}{\upsilon }_{ij}{\upsilon }_{ji}{R}_i^{\exp }{R}_j^{\exp }.$$

For each groups but 4 and 13, ${\sigma }_{ij}^{\mathrm{corr}}={\sigma }_{ik}^{\mathrm{corr}}$ for any $i\ne j\ne k$ within the group. In this way, the corresponding sub-matrices are unambiguously constructed from the data listed in Table A1, for example, for group 2, one gets:

$$||{\upsilon }_{ij}{||}_2\equiv {\mathit{V}}_2=\left(\begin{array}{ccccc}6.00& 4.84& 4.84& 4.84& 4.84\\ 4.76& 20.4& 4.76& 4.76& 4.76\\ 3.40& 3.40& 4.22& 3.40& 3.40\\ 3.40& 3.40& 3.40& 17.75& 3.40\\ 0& 0& 0& 0& 3.00\end{array}\right).$$

The sub-matrices for groups 4 and 13 are exceptional and have the form:

$$||{\upsilon }_{ij}{||}_4\equiv {\mathit{V}}_4=\left(\begin{array}{ccc}2.80& 2.80& 1.80\\ 2.80& 3.80& 2.80\\ 0.84& 0.84& 1.40\end{array}\right),$$

(22)

$$||{\upsilon }_{ij}{||}_{13}\equiv {\mathit{V}}_{13}=\left(\begin{array}{cccccccc}0.698& 0.673& 0.471& 0.471& 0.522& 0.521& 0.521& 0.519\\ 0.673& 0.699& 0.470& 0.471& 0.521& 0.521& 0.520& 0.519\\ 0.471& 0.470& 0.632& 0.604& 0.562& 0.561& 0.561& 0.559\\ 0.471& 0.471& 0.604& 0.637& 0.563& 0.562& 0.562& 0.560\\ 0.522& 0.521& 0.562& 0.563& 0.682& 0.560& 0.559& 0.558\\ 0.521& 0.521& 0.561& 0.562& 0.560& 0.680& 0.559& 0.557\\ 0.521& 0.520& 0.561& 0.562& 0.559& 0.559& 0.680& 0.556\\ 0.519& 0.519& 0.559& 0.560& 0.558& 0.557& 0.556& 0.693\end{array}\right).$$

(23)

Contribution from the relative data are taken as

$${\chi }_{\mathrm{rel}}^2=\sum _{\mathrm{g}}{\left({\mathit{R}}_{\mathrm{g}}^{\exp }-N_0{N}_{\mathrm{g}}{\mathit{R}}_{\mathrm{g}}^{\mathrm{th}}\right)}^T{\mathit{W}}_{\mathrm{g}}^{-1}\left({\mathit{R}}_{\mathrm{g}}^{\mathrm{th}}-N_0{N}_{\mathrm{g}}{\mathit{R}}_{\mathrm{g}}^{\mathrm{th}}\right),$$

where the sum is over the groups g = 11, 13, 17–19. The matrices ${\mathit{W}}_{\mathrm{g}}$ do not include the flux uncertainty term. The additional normalization factors $N_{\mathrm{g}}$ are calculated as

$$N_{\mathrm{g}}={\displaystyle \frac{{\left({\mathit{R}}_{\mathrm{g}}^{\exp }\right)}^T{\mathit{W}}_{\mathrm{g}}^{-1}{\mathit{R}}_{\mathrm{g}}^{\mathrm{th}}}{{\left({\mathit{R}}_{\mathrm{g}}^{\mathrm{th}}\right)}^T{\mathit{W}}_{\mathrm{g}}^{-1}{\mathit{R}}_{\mathrm{g}}^{\mathrm{th}}}}$$

(solution of the minimization condition $\partial {\chi }_{\mathrm{rel}}^2/\partial N_{\mathrm{g}}=0$).

The last term in Equation (21) accounts for the data on $R_{12}$ and is calculated as

$${\chi }_{12}^2=\sum _i\frac{{\left(R_{12,i}^{\exp }-R_{12,i}^{\mathrm{th}}\right)}^2}{{\sigma }_i^2},R_{12,i}^{\mathrm{th}}={\displaystyle \frac{1-L_0^2/L_{1i}^2}{1-L_0^2/L_{2i}^2}},L_{2i}>L_{1i}.$$

Here, ${\sigma }_i$ is the full uncertainty of the measurement i and sum is over the points # 1, 3, 8–10, 14, and 15 listed in Table A2. Note that the oscillation effect is almost negligible at $L\lesssim 250$ m and canceled if $L_{2i}\approx L_{1i}$ (as for the point # 9); in any case, the impact of these contributions is marginal.

All the fits are performed with the CERN function minimization and error analysis package MINUIT [206,207]. The errors of the output parameters quoted below correspond to one standard deviation.

6.5. Comparison with the Main Dataset

To study the robustness of the fitted parameters $N_0$ and $L_0$ (and their errors) with respect to different datasets, we carried out numerous fits which, in particular, involve

• full dataset,
• data at baselines $L>10$ m,
• data at baselines $L<200$ m,
• data at baselines $L=10-200$ m.

In addition, we consistently excluded from the analysis those groups of data that significantly affect the result. The main results of these analyses are summarized in

Table A3. A summary of the results of one- and two-parameter fits performed for different data subsets and ${\overline{\nu }}_e$ spectrum models; p is the right-tail p-value, $\Delta$ is the relative change (in percent) of the ${\chi }^2/$NDF values for each pair of these fits, and $N_{\sigma }$ is the roughly estimated number of standard deviations of $L_0$ from zero.

	One-Parameter Fit			Two-Parameter Fit
	${\mathit{N}}_{\mathbf{0}}$	${\chi }^2/\mathrm{NDF}$	$\mathit{p}$	${\mathit{N}}_{\mathbf{0}}$	${\mathit{L}}_{\mathbf{0}}$(m)	${\chi }^2/\mathrm{NDF}$	$\mathit{p}$	$\Delta$(%)	${\mathit{N}}_{\sigma }$
1. All data, NDF = 130 for one- and 129 for two-parameter fit
Mueller	$0.{953}_{-0.026}^{+0.027}$	$0.97$	$0.58$	$0.{957}_{-0.039}^{+0.042}$	$1.{38}_{-0.75}^{+0.46}$	$0.94$	$0.67$	$2.2$	$2.2$
Mueller + FRMII	$0.{957}_{-0.026}^{+0.027}$	$0.97$	$0.58$	$0.{960}_{-0.039}^{+0.042}$	$1.{40}_{-0.73}^{+0.45}$	$0.94$	$0.67$	$2.4$	$2.2$
Huber–Mueller	$0.{943}_{-0.026}^{+0.027}$	$0.97$	$0.58$	$0.{947}_{-0.038}^{+0.042}$	$1.{38}_{-0.74}^{+0.46}$	$0.94$	$0.67$	$2.2$	$2.2$
Huber + FRMII	$0.{946}_{-0.026}^{+0.027}$	$0.97$	$0.58$	$0.{950}_{-0.039}^{+0.042}$	$1.{40}_{-0.73}^{+0.45}$	$0.94$	$0.67$	$2.4$	$2.2$
Fallot	$1.{001}_{-0.027}^{+0.029}$	$0.97$	$0.58$	$1.{004}_{-0.041}^{+0.044}$	$1.{33}_{-0.82}^{+0.47}$	$0.96$	$0.61$	$1.8$	$1.5$
Fallot + FRMII	$1.{002}_{-0.027}^{+0.029}$	$0.97$	$0.58$	$1.{005}_{-0.041}^{+0.044}$	$1.{34}_{-0.81}^{+0.47}$	$0.96$	$0.61$	$1.8$	$1.5$
Silaeva–Sinev	$0.{971}_{-0.026}^{+0.028}$	$1.01$	$0.45$	$0.{974}_{-0.040}^{+0.043}$	$1.{26}_{-0.98}^{+0.49}$	$0.99$	$0.52$	$1.2$	$1.9$
Silaeva–Sinev + FRMII	$0.{980}_{-0.027}^{+0.028}$	$0.98$	$0.55$	$0.{983}_{-0.040}^{+0.043}$	$1.{30}_{-0.87}^{+0.48}$	$0.97$	$0.58$	$1.5$	$1.5$
RONS	$0.{973}_{-0.026}^{+0.028}$	$0.96$	$0.61$	$0.{977}_{-0.040}^{+0.043}$	$1.{38}_{-0.74}^{+0.46}$	$0.94$	$0.67$	$2.3$	$1.9$
BILL + FRMII	$0.{981}_{-0.027}^{+0.028}$	$0.97$	$0.58$	$0.{985}_{-0.040}^{+0.043}$	$1.{39}_{-0.74}^{+0.46}$	$0.94$	$0.67$	$2.3$	$2.2$
2.   $L>10$ m, NDF = 54 for one- and 53 for two-parameter fit
Mueller	$0.{954}_{-0.026}^{+0.027}$	$0.65$	$0.98$	$0.{957}_{-0.039}^{+0.042}$	$1.{44}_{-0.78}^{+0.48}$	$0.60$	$0.99$	$9.7$	$1.8$
Mueller + FRMII	$0.{957}_{-0.026}^{+0.027}$	$0.65$	$0.98$	$0.{961}_{-0.039}^{+0.042}$	$1.{45}_{-0.76}^{+0.47}$	$0.59$	$0.99$	$10.1$	$2.0$
Huber–Mueller	$0.{943}_{-0.026}^{+0.027}$	$0.65$	$0.98$	$0.{947}_{-0.039}^{+0.042}$	$1.{45}_{-0.77}^{+0.48}$	$0.59$	$0.99$	$9.8$	$2.0$
Huber + FRMII	$0.{947}_{-0.026}^{+0.027}$	$0.65$	$0.98$	$0.{950}_{-0.039}^{+0.042}$	$1.{46}_{-0.76}^{+0.47}$	$0.59$	$0.99$	$10.2$	$2.0$
Fallot	$1.{001}_{-0.027}^{+0.029}$	$0.69$	$0.96$	$1.{005}_{-0.041}^{+0.044}$	$1.{41}_{-0.82}^{+0.49}$	$0.64$	$0.98$	$7.9$	$1.8$
Fallot + FRMII	$1.{002}_{-0.027}^{+0.029}$	$0.68$	$0.97$	$1.{005}_{-0.041}^{+0.044}$	$1.{41}_{-0.82}^{+0.49}$	$0.63$	$0.98$	$8.0$	$1.8$
Silaeva–Sinev	$0.{971}_{-0.026}^{+0.028}$	$0.79$	$0.87$	$0.{974}_{-0.040}^{+0.043}$	$1.{36}_{-0.91}^{+0.50}$	$0.75$	$0.91$	$5.4$	$1.7$
Silaeva–Sinev + FRMII	$0.{980}_{-0.027}^{+0.028}$	$0.72$	$0.94$	$0.{983}_{-0.040}^{+0.043}$	$1.{39}_{-0.85}^{+0.49}$	$0.68$	$0.96$	$6.8$	$1.7$
RONS	$0.{974}_{-0.026}^{+0.028}$	$0.64$	$0.98$	$0.{977}_{-0.040}^{+0.043}$	$1.{45}_{-0.77}^{+0.48}$	$0.59$	$0.99$	$10.0$	$1.8$
BILL + FRMII	$0.{981}_{-0.027}^{+0.028}$	$0.65$	$0.98$	$0.{985}_{-0.040}^{+0.043}$	$1.{45}_{-0.77}^{+0.48}$	$0.59$	$0.99$	$9.8$	$2.0$
3.   $L<200$ m, NDF = 115 for one- and 114 for two-parameter fit
Mueller	$0.{947}_{-0.026}^{+0.028}$	$1.01$	$0.45$	$0.{953}_{-0.040}^{+0.043}$	$1.{32}_{-0.89}^{+0.49}$	$0.99$	$0.51$	$1.6$	$1.8$
Mueller + FRMII	$0.{950}_{-0.026}^{+0.028}$	$1.00$	$0.48$	$0.{956}_{-0.040}^{+0.043}$	$1.{32}_{-0.89}^{+0.49}$	$0.99$	$0.51$	$1.7$	$1.5$
Huber–Mueller	$0.{937}_{-0.026}^{+0.027}$	$1.01$	$0.45$	$0.{943}_{-0.040}^{+0.043}$	$1.{32}_{-0.89}^{+0.49}$	$0.99$	$0.51$	$1.7$	$1.8$
Huber + FRMII	$0.{939}_{-0.026}^{+0.027}$	$1.00$	$0.48$	$0.{945}_{-0.040}^{+0.043}$	$1.{33}_{-0.89}^{+0.49}$	$0.99$	$0.51$	$1.7$	$1.5$
Fallot	$0.{997}_{-0.028}^{+0.029}$	$1.02$	$0.42$	$1.{003}_{-0.042}^{+0.045}$	$1.{32}_{-0.90}^{+0.49}$	$1.01$	$0.45$	$1.6$	$1.5$
Fallot + FRMII	$0.{998}_{-0.028}^{+0.029}$	$1.02$	$0.42$	$1.{004}_{-0.042}^{+0.045}$	$1.{32}_{-0.90}^{+0.49}$	$1.00$	$0.48$	$1.6$	$1.8$
Silaeva–Sinev	$0.{971}_{-0.027}^{+0.028}$	$1.06$	$0.31$	$0.{977}_{-0.041}^{+0.044}$	$1.{31}_{-0.92}^{+0.49}$	$1.05$	$0.34$	$1.4$	$1.5$
Silaeva–Sinev + FRMII	$0.{978}_{-0.027}^{+0.029}$	$1.04$	$0.37$	$0.{984}_{-0.041}^{+0.045}$	$1.{32}_{-0.91}^{+0.49}$	$1.02$	$0.42$	$1.5$	$1.8$
RONS	$0.{967}_{-0.027}^{+0.028}$	$1.01$	$0.45$	$0.{973}_{-0.041}^{+0.044}$	$1.{32}_{-0.89}^{+0.49}$	$0.99$	$0.51$	$1.7$	$1.8$
BILL + FRMII	$0.{974}_{-0.027}^{+0.028}$	$1.01$	$0.45$	$0.{981}_{-0.041}^{+0.044}$	$1.{32}_{-0.89}^{+0.49}$	$0.99$	$0.51$	$1.7$	$1.8$
4.   $L=10-200$ m, NDF = 39 for one- and 38 for two-parameter fit
Mueller	$0.{948}_{-0.026}^{+0.028}$	$0.66$	$0.95$	$0.{954}_{-0.040}^{+0.043}$	$1.{39}_{-0.92}^{+0.51}$	$0.60$	$0.98$	$9.9$	$1.7$
Mueller + FRMII	$0.{950}_{-0.026}^{+0.028}$	$0.64$	$0.96$	$0.{956}_{-0.040}^{+0.043}$	$1.{38}_{-0.93}^{+0.51}$	$0.58$	$0.98$	$10.1$	$1.7$
Huber–Mueller	$0.{937}_{-0.026}^{+0.027}$	$0.65$	$0.96$	$0.{944}_{-0.040}^{+0.043}$	$1.{39}_{-0.93}^{+0.51}$	$0.59$	$0.98$	$10.0$	$1.7$
Huber + FRMII	$0.{940}_{-0.026}^{+0.027}$	$0.64$	$0.96$	$0.{946}_{-0.040}^{+0.043}$	$1.{38}_{-0.93}^{+0.51}$	$0.58$	$0.98$	$10.1$	$1.7$
Fallot	$0.{997}_{-0.028}^{+0.029}$	$0.72$	$0.90$	$1.{004}_{-0.042}^{+0.046}$	$1.{40}_{-0.90}^{+0.51}$	$0.66$	$0.95$	$9.1$	$1.7$
Fallot + FRMII	$0.{998}_{-0.028}^{+0.029}$	$0.72$	$0.90$	$1.{005}_{-0.042}^{+0.046}$	$1.{40}_{-0.90}^{+0.51}$	$0.66$	$0.95$	$9.2$	$1.7$
Silaeva–Sinev	$0.{971}_{-0.027}^{+0.028}$	$0.87$	$0.70$	$0.{978}_{-0.041}^{+0.044}$	$1.{42}_{-0.87}^{+0.50}$	$0.81$	$0.79$	$7.5$	$1.8$
Silaeva–Sinev + FRMII	$0.{978}_{-0.027}^{+0.029}$	$0.78$	$0.84$	$0.{984}_{-0.041}^{+0.045}$	$1.{41}_{-0.88}^{+0.50}$	$0.72$	$0.90$	$8.5$	$1.8$
RONS	$0.{967}_{-0.027}^{+0.028}$	$0.66$	$0.95$	$0.{974}_{-0.041}^{+0.044}$	$1.{39}_{-0.92}^{+0.51}$	$0.60$	$0.98$	$10.0$	$1.7$
BILL + FRMII	$0.{975}_{-0.027}^{+0.029}$	$0.65$	$0.96$	$0.{981}_{-0.041}^{+0.045}$	$1.{39}_{-0.93}^{+0.51}$	$0.59$	$0.98$	$10.0$	$1.7$
5.   All data without DANSS, NDF = 122 for one- and 121 for two-parameter fit
Mueller	$0.{953}_{-0.026}^{+0.027}$	$0.95$	$0.64$	$0.{957}_{-0.039}^{+0.042}$	$1.{45}_{-1.45}^{+0.65}$	$0.94$	$0.67$	$0.7$	$1.5$
Mueller + FRMII	$0.{957}_{-0.026}^{+0.027}$	$0.95$	$0.64$	$0.{961}_{-0.039}^{+0.042}$	$1.{48}_{-1.48}^{+0.64}$	$0.94$	$0.67$	$0.9$	$1.5$
Huber–Mueller	$0.{943}_{-0.026}^{+0.027}$	$0.95$	$0.64$	$0.{947}_{-0.039}^{+0.042}$	$1.{45}_{-1.45}^{+0.65}$	$0.94$	$0.67$	$0.8$	$1.5$
Huber + FRMII	$0.{946}_{-0.026}^{+0.027}$	$0.95$	$0.64$	$0.{950}_{-0.039}^{+0.042}$	$1.{49}_{-1.49}^{+0.64}$	$0.94$	$0.67$	$0.9$	$1.5$
Fallot	$1.{001}_{-0.027}^{+0.029}$	$0.96$	$0.61$	$1.{004}_{-0.041}^{+0.044}$	$1.{33}_{-1.33}^{+0.69}$	$0.96$	$0.61$	$0.3$	$1.0$
Fallot + FRMII	$1.{002}_{-0.027}^{+0.029}$	$0.96$	$0.61$	$1.{005}_{-0.041}^{+0.044}$	$1.{34}_{-1.34}^{+0.69}$	$0.96$	$0.61$	$0.3$	$1.0$
Silaeva–Sinev	$0.{971}_{-0.026}^{+0.028}$	$1.00$	$0.48$	$0.{973}_{-0.040}^{+0.043}$	$1.{15}_{-1.15}^{+0.76}$	$1.00$	$0.48$	$-0.2$	$1.0$
Silaeva–Sinev + FRMII	$0.{980}_{-0.027}^{+0.028}$	$0.97$	$0.58$	$0.{983}_{-0.040}^{+0.043}$	$1.{26}_{-1.26}^{+0.72}$	$0.97$	$0.58$	$0.0$	$1.0$
RONS	$0.{973}_{-0.026}^{+0.028}$	$0.95$	$0.64$	$0.{977}_{-0.040}^{+0.043}$	$1.{46}_{-1.46}^{+0.65}$	$0.94$	$0.67$	$0.8$	$1.5$
BILL + FRMII	$0.{981}_{-0.027}^{+0.028}$	$0.95$	$0.64$	$0.{985}_{-0.040}^{+0.044}$	$1.{46}_{-1.46}^{+0.65}$	$0.94$	$0.67$	$0.8$	$1.5$
6.   All data without PROSPECT, NDF = 61 for one- and 60 for two-parameter fit
Mueller	$0.{953}_{-0.026}^{+0.027}$	$0.75$	$0.93$	$0.{958}_{-0.039}^{+0.042}$	$1.{49}_{-0.67}^{+0.44}$	$0.68$	$0.97$	$9.7$	$2.2$
Mueller + FRMII	$0.{957}_{-0.026}^{+0.027}$	$0.75$	$0.93$	$0.{961}_{-0.039}^{+0.042}$	$1.{51}_{-0.66}^{+0.44}$	$0.68$	$0.97$	$10.2$	$2.2$
Huber–Mueller	$0.{943}_{-0.026}^{+0.027}$	$0.75$	$0.93$	$0.{947}_{-0.039}^{+0.042}$	$1.{50}_{-0.67}^{+0.44}$	$0.68$	$0.97$	$9.8$	$2.2$
Huber + FRMII	$0.{946}_{-0.026}^{+0.027}$	$0.75$	$0.93$	$0.{951}_{-0.039}^{+0.042}$	$1.{51}_{-0.66}^{+0.44}$	$0.68$	$0.97$	$10.3$	$2.2$
Fallot	$1.{001}_{-0.027}^{+0.029}$	$0.76$	$0.92$	$1.{005}_{-0.041}^{+0.044}$	$1.{45}_{-0.71}^{+0.46}$	$0.71$	$0.96$	$7.9$	$1.9$
Fallot + FRMII	$1.{002}_{-0.027}^{+0.029}$	$0.76$	$0.92$	$1.{006}_{-0.041}^{+0.044}$	$1.{45}_{-0.71}^{+0.46}$	$0.71$	$0.96$	$8.1$	$1.9$
Silaeva–Sinev	$0.{971}_{-0.026}^{+0.028}$	$0.83$	$0.83$	$0.{974}_{-0.040}^{+0.043}$	$1.{38}_{-0.80}^{+0.47}$	$0.79$	$0.88$	$5.4$	$1.8$
Silaeva–Sinev + FRMII	$0.{980}_{-0.027}^{+0.028}$	$0.79$	$0.88$	$0.{983}_{-0.040}^{+0.043}$	$1.{42}_{-0.75}^{+0.46}$	$0.74$	$0.93$	$6.8$	$2.0$
RONS	$0.{973}_{-0.026}^{+0.028}$	$0.74$	$0.94$	$0.{978}_{-0.040}^{+0.043}$	$1.{50}_{-0.67}^{+0.44}$	$0.67$	$0.98$	$10.0$	$2.2$
BILL + FRMII	$0.{981}_{-0.027}^{+0.028}$	$0.75$	$0.93$	$0.{985}_{-0.040}^{+0.043}$	$1.{50}_{-0.67}^{+0.44}$	$0.68$	$0.97$	$9.8$	$2.2$
7.   All data without Neutrino-4, NDF = 127 for one- and 126 for two-parameter fit
Mueller	$0.{953}_{-0.026}^{+0.027}$	$0.97$	$0.58$	$0.{957}_{-0.039}^{+0.042}$	$1.{40}_{-0.73}^{+0.46}$	$0.95$	$0.64$	$2.4$	$1.9$
Mueller + FRMII	$0.{957}_{-0.026}^{+0.027}$	$0.97$	$0.58$	$0.{960}_{-0.039}^{+0.042}$	$1.{42}_{-0.72}^{+0.45}$	$0.95$	$0.64$	$2.5$	$1.9$
Huber–Mueller	$0.{943}_{-0.026}^{+0.027}$	$0.97$	$0.58$	$0.{947}_{-0.039}^{+0.042}$	$1.{41}_{-0.73}^{+0.46}$	$0.95$	$0.64$	$2.4$	$1.9$
Huber + FRMII	$0.{946}_{-0.026}^{+0.027}$	$0.97$	$0.58$	$0.{950}_{-0.039}^{+0.042}$	$1.{42}_{-0.72}^{+0.45}$	$0.95$	$0.64$	$2.6$	$1.9$
Fallot	$1.{001}_{-0.027}^{+0.029}$	$0.98$	$0.55$	$1.{004}_{-0.041}^{+0.044}$	$1.{36}_{-0.80}^{+0.47}$	$0.96$	$0.61$	$1.9$	$1.9$
Fallot + FRMII	$1.{002}_{-0.027}^{+0.029}$	$0.98$	$0.55$	$1.{005}_{-0.041}^{+0.044}$	$1.{36}_{-0.79}^{+0.47}$	$0.96$	$0.61$	$2.0$	$1.9$
Silaeva–Sinev	$0.{971}_{-0.026}^{+0.028}$	$1.01$	$0.45$	$0.{974}_{-0.040}^{+0.043}$	$1.{28}_{-0.94}^{+0.49}$	$1.00$	$0.48$	$1.3$	$1.5$
Silaeva–Sinev + FRMII	$0.{980}_{-0.027}^{+0.028}$	$0.99$	$0.52$	$0.{983}_{-0.040}^{+0.043}$	$1.{33}_{-0.85}^{+0.48}$	$0.97$	$0.58$	$1.7$	$1.9$
RONS	$0.{973}_{-0.026}^{+0.028}$	$0.97$	$0.58$	$0.{977}_{-0.040}^{+0.043}$	$1.{41}_{-0.73}^{+0.46}$	$0.94$	$0.67$	$2.4$	$2.2$
BILL + FRMII	$0.{981}_{-0.027}^{+0.028}$	$0.97$	$0.58$	$0.{985}_{-0.040}^{+0.043}$	$1.{41}_{-0.73}^{+0.46}$	$0.95$	$0.64$	$2.4$	$1.9$
8.   All data with full Neutrino-4 dataset, NDF = 150 for one- and 149 for two-parameter fit
Mueller	$0.{953}_{-0.026}^{+0.027}$	$0.94$	$0.69$	$0.{957}_{-0.039}^{+0.042}$	$1.{40}_{-0.73}^{+0.45}$	$0.92$	$0.75$	$2.1$	$2.0$
Mueller + FRMII	$0.{957}_{-0.026}^{+0.027}$	$0.94$	$0.69$	$0.{960}_{-0.039}^{+0.042}$	$1.{41}_{-0.71}^{+0.45}$	$0.92$	$0.75$	$2.2$	$2.0$
Huber–Mueller	$0.{943}_{-0.026}^{+0.027}$	$0.94$	$0.69$	$0.{947}_{-0.039}^{+0.042}$	$1.{40}_{-0.73}^{+0.45}$	$0.92$	$0.75$	$2.1$	$2.0$
Huber + FRMII	$0.{946}_{-0.026}^{+0.027}$	$0.94$	$0.69$	$0.{950}_{-0.039}^{+0.042}$	$1.{42}_{-0.71}^{+0.45}$	$0.92$	$0.75$	$2.3$	$2.0$
Fallot	$1.{001}_{-0.027}^{+0.029}$	$0.95$	$0.66$	$1.{004}_{-0.041}^{+0.044}$	$1.{35}_{-0.79}^{+0.47}$	$0.93$	$0.72$	$1.7$	$2.0$
Fallot + FRMII	$1.{002}_{-0.027}^{+0.029}$	$0.95$	$0.66$	$1.{005}_{-0.041}^{+0.044}$	$1.{35}_{-0.79}^{+0.47}$	$0.93$	$0.72$	$1.7$	$2.0$
Silaeva–Sinev	$0.{971}_{-0.026}^{+0.028}$	$0.98$	$0.55$	$0.{974}_{-0.040}^{+0.043}$	$1.{28}_{-0.93}^{+0.49}$	$0.97$	$0.59$	$1.2$	$1.6$
Silaeva–Sinev + FRMII	$0.{980}_{-0.027}^{+0.028}$	$0.96$	$0.62$	$0.{983}_{-0.040}^{+0.043}$	$1.{32}_{-0.84}^{+0.48}$	$0.94$	$0.69$	$1.5$	$2.0$
RONS	$0.{973}_{-0.026}^{+0.028}$	$0.94$	$0.69$	$0.{977}_{-0.040}^{+0.043}$	$1.{40}_{-0.72}^{+0.45}$	$0.92$	$0.75$	$2.2$	$2.0$
BILL + FRMII	$0.{981}_{-0.027}^{+0.028}$	$0.94$	$0.69$	$0.{985}_{-0.040}^{+0.043}$	$1.{40}_{-0.72}^{+0.45}$	$0.92$	$0.75$	$2.1$	$2.0$

, which includes one- and two-parameter fits to eight datasets realized with ten models of the ${\overline{\nu }}_e$ flux. The table shows the ${\chi }^2/$NDF values, p-values, and also relative change of the reduced ${\chi }^2$s for each pair of the fits and the number of standard deviations of the best-fit $L_0$ from zero, roughly estimated as $\sqrt{{\chi }_1^2-{\chi }_2^2}$; this estimation is based on the fact of the closeness of the normalization parameters obtained in the one- and two-parameter fits. Selected results of the fits are illustrated in Figure 10, Figure 11, Figure 12, Figure 13, Figure 14, Figure 15, Figure 16, Figure 17, Figure 18 and Figure 19. Some outputs in the table and in the legends of Figure 11, Figure 12, Figure 13, Figure 14, Figure 15 and Figure 16 are written with excessive precision (only two digits in the mantissa for $N_0$ are significant). This is only needed to distinguish sometimes small differences between the results of the fits made with different data samples and spectrum models. The points related to the early experiments are explained in Figure 1. As mentioned above, each group of relative data (see Table A1) is always normalized to the best-fit curve. In all figures which involve the data points, the filled (open) symbols mark the absolute (relative) measurements. The curves and respective bands of uncertainty in the figures correspond to a reactor with pure ${}^{235}$U fuel. Nevertheless, in all calculations, we explicitly take into account the fuel compositions specific to each experiment, as they have a small, but not entirely insignificant, effect on the results. This is illustrated in Figure 10, which depicts the double ratio

$$R\left(\mathrm{real}\mathrm{fuel}\mathrm{composition}\right)/R\left(\mathrm{pure}{}^{235}\mathrm{U}\right)$$

(24)

calculated with four spectrum models (Mueller, Huber–Mueller, Silaeva–Sinev, and Fallot) for reactors with the same fuel composition as in those used in the experiments CHOOZ (Double Chooz), Palo Verde, KamLAND, Daya Bay, and RENO. In the calculations, we use the NuFIT 5.0 central values for the best-fit parameters Equation (18). It can be seen that, although the double ratio depends to some extent on the spectrum model, it differs very little from 1 (within ≲1.5%). At the baselines $L\lesssim$ 10 km, the double ratio Equation (24) is completely insensitive to variations of the fuel composition, and the sensitivity becomes noticeable at L∼60 km (slightly to the left of the first local minimum in the spectrum-weighted survival probability $\langle {\mathcal{P}}_{ee}\rangle$), at L∼200 km (to the right of the first local maximum in $\langle {\mathcal{P}}_{ee}\rangle$), and at L∼100 km (transition region between the minimum and maximum). Thus, for almost all the data except KamLAND, the calculation made for a Uranium-235 reactor is fine for visual comparison of the fits with the data. Given some uncertainty in determination of the flux-weighted baseline and in the averaged fuel composition of the reactors for the KamLAND point, even this exception is inessential. The aforesaid applies even more to the uncertainty associated with the time evolution of the isotopic fuel content in the industrial reactors: obviously, this is a subleading effect which can be safely ignored. The double ratios, calculated in the “3 + 1” framework with varying parameters $\Delta m_{41}^2$ and ${\theta }_{14}$, show practically the same shapes, which justifies comparisons shown in Figure 2 and Figure 3.

6.5.1. One-Parameter Fits (Flux Renormalization)

One-parameter fits for finding the ${\overline{\nu }}_e$ flux normalization of several model predictions are interesting in that they show the integral difference between these models and, if desired, can be easily transformed into the estimates of the limits to the mixing angle ${\theta }_{14}$. For our purposes, it is most important to study the dependency of the extracted normalization factors from the choice of the datasets and then to understand how much these factors will change compared to those obtained in the two-parameter fits (which include $L_0$) because small changes will signal that the ISLV effect is likely unrelated to some uncertainties in the ${\overline{\nu }}_e$ flux calculations but is an effect additional to the flux renormalization.

Table A3 in Appendix B collects the detailed results of the one-parameter fits for ten spectrum models and eight datasets. Figure 11, Figure 12 and Figure 13 illustrate a small part of this analysis for the full data set, which use, respectively, the Huber–Mueller, Silaeva–Sinev, and Fallot predictions for the ${\overline{\nu }}_e$ spectrum. These three models are chosen as predicting substantially different ${\overline{\nu }}_e$ energy spectra and requiring the largest, medium, and smallest renormalization, respectively. A total of 129 data points are plotted, some of which overlap due to close baselines. The obtained normalization factors, $N_0$, and corresponding values of ${\chi }^2/\mathrm{NDF}$ are shown in the legends in the top of each figure. As can be seen from the table and figures, the overall renormalization leads to quite reasonable agreement with most of the reactor data, giving wholly satisfactory ${\chi }^2$s for all spectrum models. In the case of the full dataset, all spectrum models yield similar ${\chi }^2$s, but rather different values of $N_0$ with similar errors (0.26–0.29): the Huber–Mueller model requires the largest ($5.7\%$) renormalization, while the Fallot and Fallot + FRMII models do not require any renormalization at all.

Removing the data with $L<10$ m from the fitted sample essentially improves consistency with the rest of the data for all models (from ${\chi }^2/\mathrm{NDF}$ = 0.96–1.01 to 0.64–0.79) without a change in the normalization factors and associated errors. The fits for the data subset “$L<200$ m” (i.e. without MBL and LBL data points) somewhat increase the reduced ${\chi }^2$s (to 1.00–1.06), and the fits for the subset “L = 10–200 m” give the results almost similar to those for $L>10$ m (with ${\chi }^2/\mathrm{NDF}$ = 0.64–0.87) and without a significant change of the normalization factor (for each flux model) in both cases. The observed stability of the normalization factors, in particular, means that the earlier SBL and MBL data are in broad agreement with the new precision measurements of Double Chooz, Daya Bay, and RENO. The required renormalization is expectedly very weakly sensitive to the relative SBL measurements (Neutrino-4, PROSPECT, STEREO, and DANSS) and to the “lonely extremes” (like Nucifer and ILL data points). However, as is seen from parts 2 and 6 of Table A3, the PROSPECT data noticeably affect the reduced ${\chi }^2$s due to the huge number and large scatter of the data points.

6.5.2. Two-Parameter Fits (Testing ISLV Hypothesis)

Figure 14, Figure 15 and Figure 16 show the comparison similar to the one shown in the previous figures, but for the two-parameter fits which include the ISL violating length $L_0$. Accordingly, the legends represent the obtained values of $N_0$ and $L_0$, together with the resulting values of ${\chi }^2/\mathrm{NDF}$. Bottom panels represent zooms of the very short-baseline region (below 13 m). The data from early experiments are explained in Figure 1. The detailed results for ten spectrum models and eight datasets are given in Table A3. The first thing that can be seen from the figures and Table A3 is that the reduced chi-squared values in all two-parameter fits are in fact always smaller than these from the respective one-parameter fits using the same inputs, and the best-fit values of the normalization factors are therewith almost the same within at least $0.6\%$ (and $<0.3\%$ for the four main datasets) that is well within the $1\sigma$ uncertainties (which are ≲2.9% for one- and ≲4.6% for two-parameter fits).

The main outputs for the ISLV scale $L_0$ are summed up in

Table 1. A summary of fits performed for four data subsets. Shown are the range of the best-fit values of $L_0$, formal average $\langle L_0\rangle$, 68% C.L. range of $L_0$, and error asymmetry; all the averaging are over the ten ${\overline{\nu }}_e$ spectrum models. The ranges and averagings have no statistical meaning.

Data Subset	Best-Fit ${\mathit{L}}_0$ (m)	$\langle {\mathit{L}}_0\rangle$ (m)	68% C.L. Range (m)	Asymmetry
Full dataset	1.26–1.40	$1.36$	0.57–1.82	−(0.23–0.33)
$L>10$ m	1.36–1.46	$1.43$	0.63–1.91	−(0.23–0.29)
$L<200$ m	1.31–1.33	$1.32$	0.42–1.81	−(0.29–0.31)
L = 10–200 m	1.38–1.42	$1.40$	0.48–1.90	−(0.27–0.29)

for the four datasets. The most salient result is the astonishing stability of the $L_0$ value, which is very weakly sensitive to the spectrum model and datasets, although the relative uncertainties are large and asymmetric: the lower error, ${\delta }_{-}$, ranges from $52.1\%$ to $77.8\%$ (73–98 cm), the upper error, ${\delta }_{+}$,—from $32.1\%$ to $38.9\%$ (44–51 cm), and the error asymmetry, $({\delta }_{+}-{\delta }_{-})/({\delta }_{+}+{\delta }_{-})$,—from $-0.33\%$ to $-0.23\%$ (roughly speaking, $L_0$ is easier to decrease than to increase). The best-fit values of $L_0$ are the same well within their uncertainties, for all spectrum models and investigated datasets. This indicates that the reactor data are consistent with the ISLV effect, and the latter is essentially independent of the spectrum models. To confirm this conclusion and better understand the results, we did several additional tests.

Figure 17 and Figure A1 in Appendix B show in different representations the confidence contours in the $(N_0,L_0)$ plane calculated for the four datasets and ten ${\overline{\nu }}_e$ spectrum models. The sequence of panels (left to right and top to bottom) in Figure A1 corresponds to increasing the value of the normalization factor $N_0$ for the fits of the full dataset. The tendency is about the same for other fitting datasets. As is seen from these figures, the extracted value of $L_0$ is very stable to variations of the input data and spectrum models, while the normalization factors, $N_0$, is expectedly very sensitive to the spectrum models and, to a much lesser extent, to the fitted dataset; at the same time, even the “extreme” contours (related to the smallest and largest $N_0$) intersect. One would expect that an increase in $N_0$ could partially offset the ISLV effect, i.e., the values of the parameters $L_0$ and $N_0$ can be correlated. Such a correlation does exist but is not statistically significant. All the above means that the ISLV effect is essentially model independent.

For the subsets “$L>10$ m” and “L = 10–200 m”, which exclude the highly fluctuating PROSPECT and Neutrino-4 data, as well as the “extreme” ILL and Nucifer data points, the extracted values of $L_0$ are even less sensitive to the spectrum model and, as a result, weaker correlate with the normalization factor. Why is that? The contributions to ${\chi }^2$ from numerous fluctuating (even within ∼2$\sigma$) data are rather large and the minimum of ${\chi }^2$ is largely determined by the normalization and therefore strongly depends on the spectrum model. When one excludes such contributions, the dependence on the spectral features practically disappears.

There is another reason to be more confident about the fits that exclude the data from very short baselines ($L\lesssim 10$ m). Since, as we see, the parameter $L_0$ can be as large as $1.9$ m (at $1\sigma$ level), the ISLV correction may exceed ≈3.6% at $L<10$ m and ≈7.4% at $L<7$ m. Therefore, one may expect that, at the very short distances, the next-to-leading order (NLO) ISLV correction ($\propto 1/L^4$) could become significant. However, the introduction of an additional adjustable parameter is now impractical given the large uncertainty of the VSBL data. Moreover, the VSBL region may turn out to be a transition region to another space–time behavior of the antineutrino flux, the so-called “off-shell regime” [208], which is still poorly understood and will not be discussed here. We only note that the seemingly inconsistent VSBL data may reflect an as-yet unexplored mode of the wave packet interactions. However, today this remark is purely speculative and cannot be further concretized.

Another (small) indetermination in the outputs from the VSBL data are related with the reactor core and detector dimensions, inhomogeneity of the core, and time evolution of its isotopic composition. We remind readers that the associated uncertainties are currently not important for all new VSBL experiments, but perhaps may start to deserve consideration for the DANSS experiment with further increase of statistics and improvement of control of systematics.

It is quite expected that narrowing the dataset used in the analysis increases the uncertainty of the extracted parameters, but it is more important that all the contours in Figure 17 and Figure A1 always overlap (even for spectrum models requiring the smallest and largest renormalization) and the lower bounds for $L_0$ are always positive. In this regard, it makes sense to recall that, in the absence of new SBL data from PROSPECT, STEREO, Neutrino-4, and DANSS, the lower limit on $L_0$ (but not the ISLV effect itself!) disappeared even at the $1\sigma$ level, when the ILL point was dropped from the analysis [63].

To clarify this point, we looked at cases where certain groups of data are excluded from the analysis. Most interesting examples are presented in Figure 18, where the first three panels demonstrate the cases in which data groups 16–18 are excluded, and the last panel shows the case where the narrow-binning Neutrino-4 dataset is used instead of the wide-binning one. All calculations in the figure are done with the three ${\overline{\nu }}_e$ flux models. The outputs for the ten spectrum models are collected in parts 5–8 of Table A3. The exclusion of the DANSS data leads to the most significant changes: in this case, the lower bound for $L_0$ cannot be determined and, in addition, the $ch{i}^2$ values become insignificantly different from those obtained from the one-parameter fits for all spectrum models (see part 5 of Table A3). Thus, the DANSS data are critical for testing the ISLV hypothesis. The reason for such high sensitivity to these data are obvious: it has the smallest statistical uncertainties, and excluding these data from the fit increases the statistical weight of the remaining SBL data. It is important to underline that the best-fit $L_0$ values (and upper bounds on $L_0$) remain practically unchanged even without the DANSS data. Excluding the PROSPECT and/or Neutrino-4 data, as well as replacing the four Neutrino-4 data points with the complete Neutrino-4 dataset, do not significantly affect the result, only slightly improving the lower boundaries of $L_0$. The same is also true for STEREO data, with the only caveat that these data alone are less consistent with the ISLV than with the standard ISL (see Section 6.6). It is notable that the PROSPECT data with a wider binning (14 bins instead of 70) reported in the previous publication of the PROSPECT collaboration [52] showed a bit more clearer hint of the ISLV effect. In our opinion, it would be useful to optimize the binning, which currently appears to be random. The formal ISLV fit to the full Neutrino-4 dataset alone yields $L_0=1.06$ m but with the same ${\chi }^2$ (≈18) as for the standard ISL fit.

The main inference that can be drawn from the considered exercises is that the confidence contours shown in Figure 18 perfectly overlap with each other and with the contours shown in Figure 17 and Figure A1. We can therefore conclude that the ISLV effect is not an accidental or artificial result based on a specific dataset, although the DANSS dataset most significantly affects the goodness of the fits.

As is seen, the minimax over all models obtained by using the fits to the most robust subset “$L_0>10$ m” yields the following $68\%$ confidence intervals:

$$0.63\mathrm{m}<L_0<1.91\mathrm{m},0.947<N_0<1.005.$$

(25)

From the first interval, we can roughly estimate the parameter ${\sigma }_{\mathrm{eff}}\sim \langle {\Sigma }_{\mathrm{eff}}\rangle$ representing the (multiply-averaged over the reactors and detectors) space–time overlaps of the wave packets of nuclei and particles:

$${\sigma }_{\mathrm{eff}}\sim (0.66-1.18)\mathrm{eV}⟺1/{\sigma }_{\mathrm{eff}}\sim (0.17-0.30)\mu .$$

(26)

In this estimation, we used the averaged ${\overline{\nu }}_e$ energy

$${\overline{E}}_{\nu }^{\mathrm{reactor}}=\sum _k{f}_k{\overline{E}}_k=4.2-4.4\mathrm{MeV},$$

(27)

where the mean energies ${\overline{E}}_k$ for the individual isotopes were borrowed from three model calculations reported in [86], and a typical fuel composition ($f_k=0.58$, $0.07$, $0.30$, $0.05$ for $k={}^{235}\mathrm{U}$, ${}^{238}\mathrm{U}$, ${}^{239}\mathrm{Pu}$, ${}^{241}\mathrm{Pu}$, respectively) was assumed. Although the results of [86] were derived for a specific configuration (a 1 MW${}_{\mathrm{th}}$ reactor in a 1 t detector at a standoff of 10 m measuring for 1 year for each isotope, assuming that only this isotope is fissioning), the estimation Equation (26) is robust, considering that the mean ${\overline{\nu }}_e$ energies, ${\overline{E}}_k$, are close to each other, within no more than $7\%$.

Remembering that ${\sigma }_{\mathrm{eff}}$ is by no means an antineutrino momentum dispersion, but is an effective parameter (that is a multiply averaged quantity, very nonlinearly dependent on the time–time (00) space–time (03) and transverse spatial (11, 12, 22) components of the inverse overlap tensors ${\tilde{\Re }}_s$ and ${\tilde{\Re }}_d$ defined by the WP states of all nuclei and particles participating in the ${\overline{\nu }}_e$ production and detection processes) that characterizes the ISLV effect, we should not be too surprised at the energy–momentum and spatial (“mesoscopic”) scales Equation (26) unusual in nuclear physics (see Section 4).

One more methodical test has been made to study sensitivity of the above results and conclusions to the uncertain value of the relative error of the theoretical flux prediction, ${\delta }_{\mathrm{flux}}$. Our calculations confirm that the parameters $N_0$ and $L_0$, as well as the corresponding ${\chi }^2$ values in both one- and two-parameter fits, do not depend on this parameter. The only quantity dependent of ${\delta }_{\mathrm{flux}}$ is the uncertainty in the extracted values of the normalization factor. Figure 19 illustrates dependencies of the asymmetric ($1\sigma$) errors of $N_0$ on ${\delta }_{\mathrm{flux}}$ evaluated for the three ${\overline{\nu }}_e$ spectrum models (Huber–Mueller [4,5], Silaeva–Sinev [173], and Fallot et al. [166]). It can be seen that both the upper and lower errors grow almost linearly with ${\delta }_{\mathrm{flux}}$.

There are many sub-leading uncertainties in the initial data (e.g., different detection thresholds and other kinematic cuts) that could potentially affect the results discussed, but our test calculations suggest that possible changes are well within the main uncertainty associated with the ${\overline{\nu }}_e$ flux predictions.

6.6. Comparison with Data on $R_{12}$

According to our hypothesis, we expect that, in the absence of the ordinary neutrino oscillations, the ratio $R_{12}$ behaves as

$$R_{12}={\displaystyle \frac{1-L_0^2/L_1^2}{1-L_0^2/L_2^2}}.$$

(28)

Although this ratio is a function of two independent variables, for visualization, it is convenient to represent it approximately, as a function of the single variable

$$x=\sqrt{1/L_1^2-1/L_2^2}.$$

Indeed, as is seen from Table A2, for all available data $L_2^2>L_1^2\gg L_0^2$ and thus Equation (28) can be in a very good approximation rewritten in the form

$$R_{12}\approx 1-{\left(L_{0}x\right)}^2.$$

(29)

Variable x maps closely spaced pairs $(L_1,L_2)$ and/or the pairs with $L_{1,2}\gg L_0$ to the narrow region $x\ll 1/L_0$, while the pairs with $L_1\ll L_2$ and with $L_{1,2}\ll L_0$ are mapped to the wide range $x\gg 1/L_0$ where the ISL violation effect is expected to be maximal. A comparison of Equation (29) with the data listed in Table A2 is shown in Figure 20. In addition, shown are the data form DANSS and STEREO (bins # 69–76 and # 148–152, respectively, as listed in Table A1) normalized, in each case, to the value in the bin with the shortest baseline (# 68 and 147, respectively); the errors were calculated assuming no correlations between the data. We do not show the data, which can be recalculated from the measurements of Neutrino-4 and PROSPECT (see Table A1), since the resulting values of $R_{12}$ are found to be very sensitive to the choice of the “reference bin” (the point at the distance $L_1$); in addition, the associated uncertainties are large even compared to those in the earlier measurements. The ISLV prediction is calculated using $L_0=1.{43}_{-0.95}^{+0.54}$, where the mean value is obtained from the global two-parameter fits to the reactor data with baselines $L>10$ m by averaged over the ten ${\overline{\nu }}_e$ spectrum models, and the errors represent the maximum and minimum values of the corresponding $1\sigma$ errors (see Table A3, part 2).

It is seen that the data on $R_{12}$, considered alone, are not convincing of the ISLV effect due to the large experimental uncertainties, given the small expected magnitude of the effect itself. One can only infer that the predictions based on the best fit values of $L_0$ do not contradict to the most of these data. The corresponding ${\chi }^2$s formally evaluated using different spectrum models and most reliable data from Table A2 are almost independent of the model and are comparable in magnitude for the ISLV- and ISL-based expectations. For the STEREO and DANSS data, we obtain

$$\begin{array}{cc}\hfill {\left(\frac{{\chi }_{12}^2}{\mathrm{NDF}}\right)}_{\mathrm{STEREO}}\approx & \left\{\begin{array}{cc}\hfill 1.81/5\approx 0.36& \mathrm{for}\mathrm{ISL},\hfill \\ \hfill 2.20/5\approx 0.44& \mathrm{for}\mathrm{ISLV},\hfill \end{array}\right.\hfill \\ \hfill {\left(\frac{{\chi }_{12}^2}{\mathrm{NDF}}\right)}_{\mathrm{DANSS}}\approx & \left\{\begin{array}{cc}\hfill 14.1/8\approx 1.76& \mathrm{for}\mathrm{ISL},\hfill \\ \hfill 6.05/8\approx 0.76& \mathrm{for}\mathrm{ISLV},\hfill \end{array}\right.\hfill \end{array}$$

where uncertainties of $L_0$ are not taken into account, and so the numbers for the ISLV are rather upper estimates of ${\chi }_{12}^2$/NDF. With this caveat, the STEREO data are better consistent with the classical ISL, but do not exclude the ISLV. The data from DANSS, on the contrary, are in better agreement with the ISLV, although future improvements in knowledge of the detection efficiencies and account for correlations between the data points may make a difference. The situation with the PROSPECT and Neutrino-4 data are less certain for the reasons given above. Formal calculation, which involves the data from STEREO and DANSS, and independent measurements listed in Table A2 (except for obviously outdated ones), yields

$${\left(\frac{{\chi }_{12}^2}{\mathrm{NDF}}\right)}_{\mathrm{all}}\approx \left\{\begin{array}{cc}\hfill 27.3/26\approx 1.05& \mathrm{for}\mathrm{ISL},\hfill \\ \hfill 20.9/26\approx 0.80& \mathrm{for}\mathrm{ISLV}.\hfill \end{array}\right.$$

In this calculation, the DANSS measurements evidently provide the dominant contribution. We can infer that the ISLV hypothesis describes the data on $R_{12}$ better than the “classical forecast”, although the latter is also in agreement with the data and, let us repeat, the trend is critically dependent on the only experiment. As is usually said in such cases, this result alone is indicative but inconclusive. However, the potential of the new VSBL reactor experiments in testing the ISLV is evident. The effect of the data on $R_{12}$ in the global fits which use the full reactor dataset is small, but not fully negligible as it leads to about $2\%$ decrease in the best-fit value of $L_0$. It is notable that adding the data on $R_{12}$ to the earlier dataset (shown in Figure 1) results in about $20\%$ reduction of $L_0$. The new relative data (NUCIFER, PROSPECT, STEREO, Neutrino-4, DANSS) lead to further reduction of $L_0$, which however “stabilizes”, i.e., as we saw above, stops changing significantly with variations of the datasets involved in the analysis.

7. Gallium Neutrino Anomaly

The gallium neutrino anomaly (GAA) is a deficit in the number of events caused by electron neutrinos from intense artificial ${}^{51}$Cr and ${}^{37}$Ar sources and measured by the gallium-based solar neutrino detectors GALLEX [209,210,211] and SAGE [212,213,214]. The aims of these measurements were to calibrate the detection efficiency of the detectors and to test the experimental procedures, including chemical extraction, counting, and analysis techniques. The results of the measurements are summarized in

Table 2. Ratios $R^{\mathrm{Ga}}$ of measured and predicted ${}^{71}$Ge production rates in the two GALLEX radioactive source experiments with ${}^{51}$Cr [209] (${}^{51}$Cr-1) and [210] (${}^{51}$Cr-2), and in the SAGE experiments with ${}^{51}$Cr [212,213] and ${}^{37}$Ar [214], as reported in [214]. Statistical and systematic uncertainties have been combined in quadrature.

Source	${}^{51}$Cr-1	${}^{51}$Cr-2		${}^{51}$Cr	${}^{37}$Ar
$R^{\mathrm{Ga}}=\frac{\mathrm{measured}\mathrm{rate}}{\mathrm{predicted}\mathrm{rate}}$	$1.{00}_{-0.10}^{+0.11}$	$0.{81}_{-0.10}^{+0.11}$		$0.{95}_{-0.12}^{+0.12}$	$0.{79}_{+0.09}^{-0.10}$
$R_{\mathrm{combined}}^{\mathrm{Ga}}$	$0.90\pm 0.07$			$0.86\pm 0.08$
$R_{\mathrm{weighted}\mathrm{average}}^{\mathrm{Ga}}$	$0.88\pm 0.05$

.

As explained in [214], the ratios of the measured-to-predicted rates, $R^{\mathrm{Ga}}$, given in Table 2 for the two GALLEX experiments have been revised due to changes in counter efficiencies, the solar neutrino subtraction, and the ${}^{222}$Rn cut inefficiency subtraction (see [215] for more detail). The mean ratio formally combined from the four measurements is $R^{\mathrm{Ga}}=0.88\pm 0.05$. The quality of the fit to the average value is characterized by ${\chi }^2/\mathrm{NDF}\approx 0.6$ and the goodness-of-fit $\mathrm{GoF}\approx 0.6$. The ratios shown in Table 2 have been calculated with respect to the rates estimated using the best-fit values of the cross section of the detection process Equation (30) calculated by Bahcall [216]. Haxton’s calculation [217] reduces the ratios by about 10%. In recent years, these results have been reconsidered by many authors (see, e.g., [218,219,220,221,222,223] and references therein). In

Table 3. Ratios of measured and expected ${}^{71}$Ge event rates in the four radioactive source experiments calculated with several models for the gallium cross sections for ${}^{51}$Cr and ${}^{37}$Ar neutrinos (borrowed from [222,224]).

		GALLEX		SAGE
		51Cr(G1)	51Cr(G2)	51Cr(S1)	51Ar(S2)
Bahcall	[216]	$0.95\pm 0.11$	$0.81\pm 0.11$	$0.95\pm 0.12$	$0.79\pm 0.08$
Haxton	[217]	$0.86\pm 0.13$	$0.74\pm 0.12$	$0.86\pm 0.14$	$0.72\pm 0.10$
Frekers et al.	[220]	$0.93\pm 0.11$	$0.79\pm 0.11$	$0.93\pm 0.12$	$0.77\pm 0.08$
Kostensalo et al.	[222]	$0.97\pm 0.11$	$0.83\pm 0.11$	$0.97\pm 0.12$	$0.81\pm 0.08$
Barinov et al.	[224]	$0.93\pm 0.11$	$0.80\pm 0.11$	$0.93\pm 0.12$	$0.{77}_{-0.08}^{+0.09}$

, we show the results of five such reanalyses [216,217,220,222,224].

Isotope ${}^{51}$Cr decays by electron capture to the ground state of ${}^{51}$V (90.1%) and to the 320 keV excited level (9.9%), and ${}^{37}$Ar decays 100% by electron capture to the ground state of ${}^{37}$Cl. Taking into account the atomic levels to which the transitions can occur, the chromium (argon) neutrino energy spectrum consists of four (two) lines. Detection of neutrinos is achieved through the charged-current neutrino capture reaction

$${\nu }_e+{}^{71}{\mathrm{Ga}}_{\mathrm{ground}\mathrm{state}}\to {}^{71}{\mathrm{Ge}}_{\mathrm{lowest}\mathrm{states}}+e^{-},$$

(30)

which produces a ${}^{71}\mathrm{Ge}$ nucleus in one of four (for ${}^{51}$Cr neutrinos) or five (for ${}^{37}$Ar neutrinos) lowest excited states, including the ground state. The neutrino spectra along with the pertinent fractions and gallium cross sections are listed in

Table 4. Neutrino energies ($E_{\nu }$), branching ratios ($f=f\left(E_{\nu }\right)$), and capture cross sections on gallium ($\sigma =\sigma \left(E_{\nu }\right)$) for the two artificial sources (borrowed from [224]).

Isotope	${\mathit{E}}_{\mathit{\nu }}$ (keV)	f (%)	$\mathit{\sigma }\left({10}^{-46}{\mathbf{cm}}^2\right)$	${\mathit{E}}_{\mathit{\nu }}$ (keV)	f (%)	$\mathit{\sigma }\left({10}^{-46}{\mathbf{cm}}^2\right)$
${}^{51}$Cr	752	$8.49$	$63.22\pm 1.40$	432	$0.93$	$27.14\pm 0.52$
	747	$81.63$	$62.58\pm 1.39$	427	$8.95$	$26.72\pm 0.51$
${}^{37}$Ar	813	$9.80$	$71.63\pm 1.62$	811	$90.20$	$71.35\pm 1.61$

taken from [224]. The cross sections are model dependent (see, e.g., [222]).

If the ISLV mechanism really works for reactor’s ${\overline{\nu }}_e$s, it must also work for the Cr and Ar ${\nu }_e$s albeit on a different spatial scale (considering the differences in the neutrino energies and production/detection processes). From the extended GS theorem, it follows that the corresponding scale is of the order of

$$L_0^{\mathrm{Ga}}=〈E_{\nu }{\left({\Sigma }_{\mathrm{eff}}\right)}^{-2}〉\sim {\overline{E}}_{\nu }^{\mathrm{Ga}}{\left({\sigma }_{\mathrm{eff}}^{\mathrm{Ga}}\right)}^{-2},$$

where ${\sigma }_{\mathrm{eff}}^{\mathrm{Ga}}\sim \langle {\Sigma }_{\mathrm{eff}}\rangle$ is an effective parameter similar to that for the reactor data and

$${\overline{E}}_{\nu }^{\mathrm{Ga}}=\frac{\langle E_{\nu }f\left(E_{\nu }\right)\sigma \left(E_{\nu }\right)\rangle }{\langle f\left(E_{\nu }\right)\sigma \left(E_{\nu }\right)\rangle }$$

(31)

is the mean neutrino energy (defined through the neutrino energy lines, decay fractions, and capture cross sections listed in Table 4). The value of ${\sigma }_{\mathrm{eff}}^{\mathrm{Ga}}$ does not have to be the same (even in order of magnitude) as the value extracted from the reactor data. Moreover, it can be in general different for the chromium and argon ${\nu }_e$ sources. However, in order to check, at least at a rough qualitative level, the possibility that the GNA can be relevant to the ISLV, let us speculate that, for both ${\nu }_e$ sources, ${\sigma }_{\mathrm{eff}}^{\mathrm{Ga}}$ is of the same order of magnitude as the reactor ${\sigma }_{\mathrm{eff}}$. Then, using Equations (27) and (31), the effective length can be crudely estimated to be

$$L_0^{\mathrm{Ga}}\sim \frac{L_0^{\mathrm{reactor}}{\overline{E}}_{\nu }^{\mathrm{Ga}}}{{\overline{E}}_{\nu }^{\mathrm{reactor}}}=\left\{\begin{array}{cc}\hfill {24}_{-16}^{+10}\mathrm{cm}& \mathrm{for}{}^{51}\mathrm{Cr},\hfill \\ \hfill {27}_{-18}^{+11}\mathrm{cm}& \mathrm{for}{}^{37}\mathrm{Ar},\hfill \end{array}\right.$$

where the errors are estimated at $68\%$ C.L. Surprisingly, these ranges are comparable in order of magnitude with the average neutrino paths in the gallium detectors. Therefore, it makes sense to undertake a bit more detailed analysis.

Figure 21 schematically shows a sectional view of the detectors GALLEX and SAGE (both are cylindrically symmetrical) and their dimensions, including the dimensions of the cylindrical openings with radioactive sources. Obviously the approximation Equation (13) is not applicable for the zone of the fiducial volume, $S_{\mathrm{cut}}$, close to the source (the L0 ISLV correction is clearly insufficient for $L\sim L_0$ and the neutrino path is not much larger than the size of the neutrino source). Moreover, the general expression for the count rate Equation (5) becomes incorrect since the trivial but important condition (i) (Section 4) does not hold in these zones. Therefore, a bit more general formulas and much more technically complex analysis are required to study the gallium anomaly in terms of the ISLV effect, as well as (it is pertinent to emphasize here) to quantitatively explain the anomaly in terms of the sterile neutrino hypothesis. The latter circumstance is simply ignored (without any reservations) in the related studies based on the standard QM approach. Below, we restrict ourselves to very crude estimates, remaining within the same framework as in the case of the reactor anomaly.

Let us roughly represent $S_{\mathrm{cut}}$ as a sphere of radius $r_{\mathrm{cut}}$ which is determined by the same effective spatial scale $L_0^{\mathrm{Ga}}\sim$ 20–30 cm responsible for the ISLV effect. Thus, although we cannot estimate the ISLV effect in the whole fiducial volume V of GALLEX or SAGE, we can do it for the volume $V_{\mathrm{cut}}$ with the cut-out sphere $S_{\mathrm{cut}}$ surrounding the source region. Clearly, the larger the radius of the sphere $S_{\mathrm{cut}}$, the more accurate the estimate. While this will not be a complete solution to the problem (especially for the small detector SAGE), it will allow us to check the consistency of the assumptions made.

In this formulation and assuming no mixing with sterile neutrinos, we may expect the following measured-to-predicted ratio due to the ISLV effect:

$$R_{\mathrm{cut}}^{\mathrm{Ga}}={\int }_{V_{\mathrm{cut}}}\frac{d\mathit{r}}{{|\mathit{r}|}^2}\left[1-{\left(\frac{L_0^{\mathrm{Ga}}}{|\mathit{r}|}\right)}^2\right]{\left({\int }_{V_{\mathrm{cut}}}\frac{d\mathit{r}}{{|\mathit{r}|}^2}\right)}^{-1}.$$

(32)

Here, the origin is positioned at the center of the source; the shape of the volume $V_{\mathrm{cut}}$ in the GALLEX and SAGE detectors is clear from Figure 21. Let us parametrize the radius of the sphere $S_{\mathrm{cut}}$ as follows:

$$r_{\mathrm{cut}}=(1-\alpha ){\overline{E}}_{\nu }/{\sigma }_{min}^2.$$

(33)

In the calculations, we use ${\sigma }_{min}=0.6$ eV (this is the lower bound for the reactor ${\sigma }_{\mathrm{eff}}$) and $0\le \alpha \le 0.5$ (this is just an artificial parameter that makes it easier to visualize the dependency of the ratio Equation (32) on the size of the cut volume).

Figure 22 shows the comparison of the ratios $R_{\mathrm{cut}}^{\mathrm{Ga}}$ as functions of ${\sigma }_{\mathrm{eff}}^{\mathrm{Ga}}$, evaluated with several values of the radius $r_{\mathrm{cut}}$. An increase of the radius leads to a decrease of the volume $V_{\mathrm{cut}}$ and, as a consequence, to an increase of the neutrino mean path and the ratio $R_{\mathrm{cut}}^{\mathrm{Ga}}$. The correction becomes ∼1 (which is meaningless for the first term of the asymptotic series) at the inner edge of the remaining fiducial volume for the given range of ${\sigma }_{\mathrm{eff}}^{\mathrm{Ga}}$. Thus, the calculations extending around and below these borders are very rough extrapolations of the LO calculations. The curves for the G1 and G2 experiments are indistinguishable in the physical regions since the only difference in the calculations for these experiments is in the positions of the chromium source ($h_1$ and $h_2$ in Figure 21). Since the SAGE tank is much smaller compared to the GALLEX one, the mean neutrino path is shorter and therefore the ISLV effect is larger. On the other hand, the volume $V_{\mathrm{cut}}$ in SAGE is also very small, which makes comparison with the SAGE data more uncertain.

The SAGE and GALLEX data are represented in Figure 22 by the wide shaded bands; the upper and lower boundaries on the measured ratios are model-dependent and are defined here as, respectively, $max(R+\Delta R)$ and $min(R-\Delta R)$, where the maximum and minimum are taken over the $1\sigma$ uncertainties of the values of $R=R_{\mathrm{cut}}^{\mathrm{Ga}}$ listed in Table 3; the lower bound is given by the Haxton model [217] and the upper bound—by the model of Kostensalo et al. [222]. Semi-transparent rectangles in the figure enclose the allowed range of ${\sigma }_{\mathrm{eff}}$ obtained from our global analysis of the reactor data (see Equation (26), Section 6). It is indicative that it is in this range that an apparent correlation between the predicted curves and gallium data is observed. Considering that the ISLV effect is very sensitive to variation of ${\sigma }_{\mathrm{eff}}$, this observation a posteriori validates our assumption that this parameter may be of the same order of magnitude in the reactor and gallium experiments. However, it is quite possible that we are dealing with just a coincidence here.

Sadly, the forced roughness of our estimations does not allow us to draw more definite conclusions. It is hoped that the new-generation experiments with artificial neutrino sources, like BEST/BEST-2 [225,226], SOX/CeSOX [227,228,229,230], and CeLAND [231,232], as well as a more rigorous analysis of the data, will make it possible to test the ISLV hypothesis more directly and conclusively. Although the primary aim of these experiments is to clarify the sterile neutrino puzzle [64], they are quite suitable for this purpose as well. However, of course, a dedicated experimental setup (large sectioned detector) would be highly desirable.

8. Conclusions

The QFT approach to neutrino oscillations predicts that the classical inverse-square law can be violated at short, but possibly macroscopic distances from the neutrino source. The numerical value of the spatial scale at which the deviation becomes essential cannot be predicted from the present-day theory, but it can be extracted, under reasonable assumptions, from the data of the past and current reactor antineutrino experiments. Our statistical analysis of the IBD event rates is based on the conventional $3\nu$ mixing scenario with the oscillation parameters adjusted to the current solar, atmospheric, accelerator, and reactor data, and supplemented with the ISLV factor and proper ${\overline{\nu }}_e$ flux renormalization. The latter has a double meaning in the analysis because it can be associated either with the uncertainties in the ${\overline{\nu }}_e$ flux calculations or (to a certain degree of precision) with the effect of light (eV-mass scale) sterile neutrinos. The results of the analysis demonstrate that the averaged over the reactor antineutrino spectrum value of the ISLV spatial scale ($L_0$) ranges from ≈0.6 to ≈1.9 m (68% C.L.) that approximately corresponds to 0.7–1.2 eV for the spectrum and reactor/detector averaged effective parameter ${\sigma }_{\mathrm{eff}}$, reciprocal of which, $1/{\sigma }_{\mathrm{eff}}\sim$ (0.2–0.3) $\mathsf{\mu }$, characterizes, in order of magnitude, the spatial-temporal overlap between the wave packets of the particles and nuclei involved into the antineutrino production and detection processes. These results are roughly consistent, within uncertainties, with our earlier estimates [62], which did not include the data from the new reactor experiments, the most influential of which is DANSS.

In addition, the best-fit value of $L_0$ is very stable with respect to the choice of the reactor data subsets and ${\overline{\nu }}_e$ spectrum models. To check the latter, we performed detailed calculations with the models based on different methods, as well as with combinations of these models and the cumulative ${\overline{\nu }}_e$ spectrum from ${}^{238}$U fission measured with the scientific neutron source FRM II in Garching. We conclude from these exercises that the best-fit value of $L_0$ is almost insensitive (well within the $1\sigma$ errors) to the features of the ${\overline{\nu }}_e$ spectra. It is needless to say that this is not the case for the flux normalization parameter $N_0$, required to match the reactor data, which varies within relatively wide range (0.95–1.01), reflecting the uncertainties in the modern calculations of the reactor ${\overline{\nu }}_e$ flux or hinting at new physics. This parameter is weakly sensitive to exclusion of different data subsets, which indicates the consistency of the data obtained in the majority of the experiments. The Huber–Mueller model (most popular in the field) requires the most renormalization, while the flux predicted by Fallot et al. does not need it at all and therefore does not support the sterile neutrino hypothesis. However, both models (like all other tested models) support the ISLV hypothesis, yielding very similar best-fit values and allowed ranges for the parameter $L_0$. The current reactor data cannot definitely confirm or exclude the light sterile neutrino hypotheses; they also do not provide unambiguous evidence in favor of the ISLV effect, but the latter is in fairly good agreement with most of the data. To put this another way, the sterile neutrino hypothesis is not the only possible explanation of the RAA, and the ISLV can be regarded as an alternative or additional effect, fully, or at least in part responsible for the observed anomaly.

Our very preliminary analysis provides a hint that the ISLV effect can be also relevant to the gallium neutrino anomaly, therewith with approximately the same allowed range of values for the parameter ${\sigma }_{\mathrm{eff}}$ as in the reactor experiments. Although the latter can be considered an accidental coincidence, it would be very interesting to confirm or disconfirm the ISLV in the next-generation experiments with artificial neutrino sources, since the current level of the gallium data does not allow more definite conclusions.