Improved Constraints on Sterile Neutrino Mixing from Disappearance Searches in the MINOS, MINOS+, Daya Bay, and Bugey-3 Experiments

Searches for electron antineutrino, muon neutrino, and muon antineutrino disappearance driven by sterile neutrino mixing have been carried out by the Daya Bay and MINOS+ collaborations. This Letter presents the combined results of these searches, along with exclusion results from the Bugey-3 reactor experiment, framed in a minimally extended four-neutrino scenario. Significantly improved constraints on the $\theta_{\mu e}$ mixing angle are derived that constitute the most stringent limits to date over five orders of magnitude in the sterile mass-squared splitting $\Delta m^2_{41}$, excluding the 90\% C.L. sterile-neutrino parameter space allowed by the LSND and MiniBooNE observations at 90\%~CL$_s$ for $\Delta m^2_{41}<5\,$eV$^2$.Furthermore, the LSND and MiniBooNE 99\%~C.L. allowed regions are excluded at 99\%~CL$_s$ for $\Delta m^2_{41}$~$<$~1.2 eV$^2$.

Searches for electron antineutrino, muon neutrino, and muon antineutrino disappearance driven by sterile neutrino mixing have been carried out by the Daya Bay and MINOS+ collaborations. This Letter presents the combined results of these searches, along with exclusion results from the Bugey-3 reactor experiment, framed in a minimally extended four-neutrino scenario. Significantly improved constraints on the θµe mixing angle are derived that constitute the most stringent limits to date over five orders of magnitude in the sterile masssquared splitting ∆m 2 41 , excluding the 90% C.L. sterile-neutrino parameter space allowed by the LSND and MiniBooNE observations at 90% CLs for ∆m 2 41 < 5 eV 2 . Furthermore, the LSND and MiniBooNE 99% C.L. allowed regions are excluded at 99% CLs for ∆m 2 41 < 1.2 eV 2 .
Since the discovery of neutrino oscillations two decades ago [1,2], the progress achieved allows more rigorous tests than ever to be conducted. The wealth of experimental data to date overwhelmingly demonstrates that the weak and mass eigenstates of the neutrino mix.
Most of the measurements made so far with solar, atmo-spheric, reactor and accelerator neutrinos [3][4][5][6][7][8][9][10][11][12][13][14][15] can be fully explained with three neutrino states that mix as described by the PMNS formalism [16][17][18]. There are, however, some experimental observations that cannot be accommodated in the three-neutrino mixing model, such as the excess of electronlike events in a muon (anti)neutrino beam observed over short baselines by the Liquid Scintillator Neutrino Detector (LSND) [19] and MiniBooNE [20,21] experiments. These observations may be explained by mixing with at least one additional fourth neutrino state with mass-squared splitting ∆m 2 41 |∆m 2 32 |, where the ∆m 2 ji = m 2 j − m 2 i represent neutrino mass-squared differences and m i is the mass of the i-th mass eigenstate. The addition of such states, a natural occurrence in many extensions of the Standard Model that incorporate neutrino masses, results in new neutrino states that are commonly deemed to be sterile in accordance with the tight constraints from precision electroweak measurements [22,23] on the number of neutrinos that couple to the Z boson. The far-reaching implications of sterile neutrinos in particle physics and cosmology makes their possible existence one of the key questions in physics.
Sterile neutrinos could be detected in oscillation experiments as a deviation from the standard three-neutrino oscillation behavior if they mix with the three active neutrinos. In 2016, the Daya Bay and MINOS experiments reported limits on active-to-sterile oscillations obtained by combining the results of their electron antineutrino and muon (anti)neutrino disappearance measurements, respectively [24], with those from the Bugey-3 experiment [25]. This Letter presents significantly improved limits obtained by utilizing a data set with roughly twice the exposure in the case of Daya Bay [26], and by adding 5.80 × 10 20 protons-on-target (POT) of MINOS+ data, recorded with the medium-energy configuration of the NuMI beam [27], to the full MINOS data sample [28]. Some key systematic uncertainties are reduced in the case of Daya Bay, and a new two-detector fit technique is employed for MINOS and MINOS+. The resulting limits provide leading constraints on possible mixing between active and sterile neutrinos, and can be used to examine the sterile neutrino interpretation of the appearance claims made by the LSND and MiniBooNE experiments in a way that is independent of CP violation and mass-ordering effects.
The results of the combined analysis presented in this Letter are interpreted within the framework of a 3+1 model, which includes one new mass eigenstate and one sterile weak eigenstate in addition to the three known mass eigenstates and active neutrino flavors. We parameterize the extended 4×4 unitary matrix U describing mixing between weak and mass eigenstates following Ref. [29], and the expressions for the elements of U that are relevant to this work become |U e3 | 2 = cos 2 θ 14 sin 2 θ 13 , |U µ4 | 2 = sin 2 θ 24 cos 2 θ 14 .
Under the assumption of neutrino-antineutrino invariance, in the ∆m 2 41 |∆m 2 31 | approximation for Daya Bay and Bugey-3 baselines, the survival probability of electron antineutrinos with energy E after traveling a distance L approximates to which yields the following sin 2 2θ 14 -dependent expression: Long-baseline experiments like MINOS and MINOS+ constrain sin 2 θ 24 by looking for muon neutrino and antineutrino disappearance, for which we can approximate the survival probability as In addition, long-baseline experiments can also look for deficits of neutral-current (NC) neutrino interactions between the Near and Far detectors, approximately described by Besides sensitivity to both θ 24 and ∆m 2 41 , the NC channel provides sensitivity to θ 14 , θ 34 and δ 24 . Sterile (anti)neutrinodriven muon to electron (anti)neutrino appearance at short baselines has been advanced as a possible explanation of the LSND and MiniBooNE results. This appearance probability is described by where Therefore, electron antineutrino disappearance constraints from reactors on sin 2 2θ 14 , combined with muon neutrino and antineutrino disappearance constraints from long-baseline experiments on sin 2 θ 24 , can place stringent constraints on the quadratically-suppressed electron neutrino or antineutrino appearance described by sin 2 2θ µe within the framework of the 3+1 model [30]. While Eqs. 3 and 4 show leading terms to illustrate the general behavior of the oscillation probabilities, exact formulae of the full survival probabilities are used in the analyses reported in this Letter. The Daya Bay reactor antineutrino experiment consists of eight identically-designed antineutrino detectors (ADs) placed in three underground experimental halls (EHs) at different distances from three pairs of 2.9 GW th nuclear reactors in the southeast of China. The two near halls, EH1 and EH2, house two ADs each and have flux-averaged baselines on the order of 550 m. The far hall, EH3, houses four ADs and has a flux-averaged baseline around 1600 m. The overburdens of EH1, EH2 and EH3 are 250, 265, and 860 meters-waterequivalent, respectively. Electron antineutrinos are detected via the inverse beta decay (IBD) reaction,ν e + p → e + + n, whose two products are visible in the ADs. Further details about the Daya Bay experiment can be found in Ref. [31].
Daya Bay's unique configuration with multiple baselines makes it well suited to search for sterile neutrino mixing. A relative comparison of the flux and spectral shape of reactor antineutrinos observed in the EHs at different baselines provides most of the sensitivity to sterile neutrino oscillations in the 10 −3 eV 2 |∆m 2 41 | 0.3 eV 2 region. For |∆m 2 41 | 0.3 eV 2 , the oscillations are too fast to be resolved by the detectors, and the sensitivity arises primarily from comparing the measured flux with the expectation. The uncertainty in the expected reactor antineutrino flux is conservatively set to 5% as motivated by recent re-evaluations in light of the so-called reactor antineutrino anomaly [32,33].
A new search for light sterile neutrino mixing was performed at Daya Bay with a data set acquired over 1230 days. This represents a factor of ∼2 increase in exposure over the previous result [34]. The analysis of this data set incorporates other improvements, such as a more precise background assessment, the inclusion of a time-dependent correction for spatial nonuniformity within each AD, and a reduction in the relative detection efficiency uncertainty to 0.2%, which is the dominant source of systematic error. The IBD selection, background rejection, and assessment of systematic uncertainties for this data set are described in detail in Ref. [26]. The normal mass ordering is assumed for ∆m 2 31 and ∆m 2 41 . The results reported here are largely insensitive to this choice.
The same two complementary methods applied in previous sterile neutrino searches at Daya Bay [34,35] are used to set the exclusion limits in the (∆m 2 41 , sin 2 2θ 14 ) parameter space. The first one is based on a purely relative comparison between the near and the far data, and relies on the frequentist approach proposed by Feldman and Cousins to determine the exclusion limits [36]. The second one uses the predicted antineutrino spectra to simultaneously fit the observations in the three halls, and uses the CL s statistical method [37,38] to set the limits.
The CL s method is a two-hypothesis test, here used to discriminate between the three-neutrino (3ν) and four-neutrino (4ν) scenarios where each combination of (∆m 2 41 , sin 2 2θ 14 ) is treated as a separate 4ν scenario. We define the test statistic ∆χ 2 = χ 2 4ν − χ 2 3ν , where χ 2 3ν is the χ 2 resulting from a fit to the 3ν hypothesis (with free θ 13 ) and χ 2 4ν is the χ 2 from a fit to the 4ν hypothesis (with free θ 13 , and ∆m 2 41 and θ 14 set to the corresponding 4ν scenario under consideration). Other parameters, namely sin 2 2θ 12 , ∆m 2 21 and |∆m 2 32 |, are constrained using external data [22]. We produce a ∆χ 2 3ν distribution by fitting simulated pseudo-experiments with ∆m 2 41 = sin 2 2θ 14 = 0 and θ 13 fixed to the best-fit value in the data. The same is done to construct a ∆χ 2 4ν distribution for every point in the (∆m 2 41 , sin 2 2θ 14 ) parameter space. Since these distributions are normally distributed, we estimate their mean and variance from Asimov data sets [39], greatly reducing the amount of computation needed. For each point in (∆m 2 41 , sin 2 2θ 14 ) the observed ∆χ 2 obs is compared to the ∆χ 2 3ν and ∆χ 2 4ν distributions in order to obtain the corresponding p-values. The CL s statistic is defined by where p H is the p-value for hypothesis H. The 90% exclusion contour is obtained by requiring CL s ≤ 0.1.
As seen in Fig. 1, consistent results are obtained by the two methods. It has been shown that the CL s approach can yield more stringent contours than the Feldman-Cousins approach with null data sets [39]. Moreover, a study using a very large number of simulated experiments found that the purely relative near-far comparison method that is used to produce the Feldman-Cousins contours had slightly lower sensitivity in the ∆m 2 41 2 × 10 −3 eV 2 region than the method where the near and far observations are fit simultaneously. This study also found that the two methods can react slightly differently to statistical fluctuations. The small differences observed in Fig. 1 are thus well within expectation.
A CL s -based analysis is also applied to the published data from the Bugey-3 experiment [25]. This reactor experiment operated at shorter (<100 m) baselines, allowing it to provide valuable constraints on sterile neutrino mixing from electron antineutrino disappearance for higher values of ∆m 2 41 compared to Daya Bay. The same methodology detailed in Ref. [24] was followed to generate the exclusion contour for Bugey-3. The main adjustments made with respect to the original Bugey-3 analysis were: (i) the use of the Gaussian CL s method, instead of the raster scan technique; (ii) the use of an updated neutron lifetime in the IBD cross-section calculation; and (iii) the use of the Huber+Mueller [40,41] model instead of the original ILL+Vogel model [42,43] to make the flux prediction at the different baselines. The reproduced contour is very similar to the one published originally by the Bugey-3 collaboration, shown in Fig. 1.
The MINOS and MINOS+ experiments used two detectors placed on the NuMI beam axis, the Near Detector (ND), located 1.04 km downstream from the production target at Fermilab at a depth of 225 meters-water-equivalent, and the Far Detector (FD), located 734 km further downstream, in the Soudan Underground Laboratory in Minnesota at a depth of 2070 meters-water-equivalent. The detectors were functionally-identical magnetized, tracking, sampling calorimeters composed of steel-scintillator planes read out by multi-anode photomultiplier tubes [44]. The NuMI neutrino beam is produced by colliding 120 GeV protons accelerated by the Main Injector complex at Fermilab with a graphite target. The emerging secondary beam of mostly π and K mesons is focused by two parabolic electromagnetic horns and allowed to decay in a 675 m long helium-filled pipe, resulting in a neutrino beam composed predominantly of ν µ , with a 1.3% contamination of ν e [45]. The detectors accumulated a 10.56 × 10 20 POT beam exposure during the MINOS neutrino runs, with the observed neutrino energy spectrum peaked at 3 GeV. In the MINOS+ phase, the detectors sampled a higher-intensity NuMI beam, upgraded as part of the NOvA experiment [46], with the neutrino energy peaked at 7 GeV. The higher-energy neutrino beam, although less favorable for three-flavor oscillation measurements (for MINOS' baseline and three-neutrino standard oscillations, the muon neutrino disappearance maximum occurs at E ν ≈ 1.6 GeV), provides greater sensitivity to sterile-induced muon neutrino disappearance by increasing the statistics in regions of L/E ν where oscillations driven by large mass-squared splittings would occur. A new search for sterile neutrino mixing using an additional exposure of 5.80 × 10 20 POT of MINOS+ data has been recently published [28]. Unlike the previous MINOS analysis that was based on the ratio between the measured neutrino energy spectra in the two detectors (Far-over-Near ratio) [47][48][49][50] and that was limited by the statistical error of the lowerstatistics FD sample, the new analysis employs a two-detector fit method, simultaneously fitting the reconstructed neutrino energy spectra in both detectors [51]. The new technique exploits the full power of the large ND statistics for L/E ν regions probed by the ND baseline.
The analysis employs both the charged-current (CC) ν µ and the NC data samples from MINOS and MINOS+. The CC ν µ disappearance channel has sensitivity to θ 24 and ∆m 2 41 , in addition to the three-flavor oscillation parameters ∆m 2 32 and θ 23 . The NC sample adds nontrivial sensitivity to θ 34 , θ 24 and ∆m 2 41 , albeit with a worse energy resolution (due to the missing energy carried by the outgoing final-state neutrino) than in the CC case, as well as lower statistics due to the lower NC interaction cross section. As detailed in Refs. [28,51], the analysis is approximately independent of the angle θ 14 and the phases δ 13 , δ 14 , and δ 24 , so these parameters are all set to zero in the fit. The MINOS and MINOS+ combined search for sterile neutrinos places the most stringent limit to date on the mixing parameter sin 2 θ 24 for most values of the sterile neutrino mass-splitting ∆m 2 41 > 10 −4 eV 2 . Following the same approach used in the first joint analysis by MINOS and Daya Bay [24], the CL s contours for the new two-detector fit of MINOS and MINOS+ data are obtained using a similar prescription to the one used by Daya Bay, but where the test statistics ∆χ 2 3ν and ∆χ 2 4ν are approximated by MC simulations of pseudo-experiments without assuming they have Gaussian distributions. The consistency with the published Feldman-Cousins corrected limits is displayed in Fig. 2. The new MINOS and MINOS+ limits are combined with the Daya Bay and Bugey-3 limits described above to obtain a new improved limit on anomalous ν µ to ν e oscillations, as discussed below.
The disappearance measurements from the three experiments are combined using the same methodology as in Ref. [24]. For each fixed value of ∆m 2 41 , the ∆χ 2 obs value and the ∆χ 2 3ν and ∆χ 2 4ν distributions for each (sin 2 2θ 14 , ∆m 2 41 ) point from the Daya Bay and Bugey-3 combination are paired with those for each (sin 2 θ 24 , ∆m 2 41 ) point from the MINOS and MINOS+ experiments, resulting in specific (sin 2 2θ µe , ∆m 2 41 ) combinations according to Eq. 7. Since systematic uncertainties of accelerator and reactor experiments are largely uncorrelated, the combined values of ∆χ 2 obs are obtained by simply summing the corresponding values from the reactor and accelerator experiments. Similarly, the combined ∆χ 2 3ν and ∆χ 2 4ν distributions are calculated by random sampling the distributions from each experiment and summing. Since several different combinations of (sin 2 2θ 14 , sin 2 θ 24 ) can yield the same sin 2 2θ µe , the combination with the largest CL s value is conservatively selected to be used in the final result.
The new combined 90% and 99% CL s limits from searches for sterile neutrino mixing in MINOS, MINOS+, Daya Bay, and Bugey-3 in the 3+1 neutrino model are shown in Figs. 3 and 4, respectively. Constraints on the sin 2 2θ µe electron (anti)neutrino appearance parameter are provided over 7 or- ders of magnitude in the sterile mass-squared splitting ∆m 2 41 . These limits are the world's most stringent over 5 orders of magnitude, for ∆m 2 41 10 eV 2 . The new constraints exclude the entire 90% C.L. allowed regions from LSND and MiniBooNE for ∆m 2 41 < 5 eV 2 , with regions at higher values being excluded by NO-MAD [54]. Further, the 99% C.L. allowed regions from LSND and MiniBooNE are excluded for ∆m 2 41 < 1.2 eV 2 . The allowed region from a global fit to data from sterile neutrino probes, intentionally excluding MINOS, MINOS+, Daya Bay, and Bugey-3 contributions, computed by the authors of Refs. [55,56], is fully excluded at the 99% C.L. The allowed region resulting from a fit to all appearance data, updated by the authors of Ref. [57] to include the MiniBooNE 2018 results [21], is equally strongly excluded. The new limits presented here thus significantly increase the tension between pure sterile neutrino mixing explanations of appearance-based indications and the null results from disappearance searches. The sole consideration of additional sterile neutrino states cannot resolve this tension, which stems from the non-observation ofν e and (−) ν µ disappearance beyond what is expected from the three-neutrino mixing model. This inconsistency may be further quantified in additional detector exposures in the process of being analyzed, specifically the last year of MINOS+ data taking, representing an additional sample of similar size to the one used here, as well as over two more years of Daya Bay data.  Bugey-3 combined 90% CLs limit on sin 2 2θµe to the LSND and MiniBooNE 90% C.L. allowed regions. Regions of parameter space to the right of the red contour are excluded. The regions excluded at 90% C.L. by the KARMEN2 Collaboration [53] and the NOMAD Collaboration [54] are also shown. The combined limit also excludes the 90% C.L. region allowed by a fit to global data by Gariazzo et al. where MINOS, MINOS+, Daya Bay, and Bugey-3 are not included [55,56], and the 90% C.L. region allowed by a fit to all available appearance data by Dentler et al. [57] updated with the 2018 MiniBooNE appearance results [21].  Bugey-3 combined 99% CLs limit on sin 2 2θµe to the LSND and MiniBooNE 99% C.L. allowed regions. The limit also excludes the 99% C.L. region allowed by a fit to global data by Gariazzo et al. where MINOS, MINOS+, Daya Bay, and Bugey-3 are not included [55,56], and the 99% C.L. region allowed by a fit to all available appearance data by Dentler et al. [57] updated with the 2018 MiniBooNE appearance results [21].