Search for a sub-eV sterile neutrino using Daya Bay’s full dataset

,

The existence of neutrino oscillation has been firmly demonstrated by many experiments, ruling out the original Standard Model's postulate that all neutrino states are massless.Although most of the findings agree with the hypothesis that there are three neutrino mass states, there exist a few anomalies [1,2] and indications [3] that may be explained by the existence of extra neutrino mass states [4].Precision measurements of the Z-boson width are consistent with three light neutrino species that participate in the weak interaction [5] so any additional neutrino species must be "sterile", that is, not subject to the weak interaction.
As for the mass of an extra neutrino state, in theory any value is possible.A mass as large as 10 15 GeV is considered by the seesaw mechanism, which can both generate the very light active-neutrino masses and produce the baryon asymmetry of the universe [6][7][8][9][10].In contrast, a keV-range sterile neutrino is a possible candidate for warm dark matter [11].The most stringent limits on the mass of a light relativistic sterile neutrino again come from cosmology, which is sensitive to the sum of neutrino masses [4].The current limit of m ν < 0.12 eV at 95% C.L. [12] remains consistent with the existence of a sub-eV sterile neutrino; this mass constraint is loosened to a few MeV when sterile neutrino self-interactions are allowed [13].
In the minimal "3+1" extension of the three-neutrino model, considering one sterile neutrino in addition to the three active neutrinos, the flavor eigenstates ν α (α = e, µ, τ, s) are related to the four mass eigenstates ν i as: where U is a unitary 4 × 4 mixing matrix, analogous to the Pontecorvo-Maki-Nakagawa-Sakata (PMNS) matrix from the three-neutrino scenario.The matrix is in general parameterized [14,15] by six mixing angles θ ij and three CP-violating phases δ i .The survival probability of electron antineutrinos (P νe νe ) is a function of the neutrino energy E and the distance traveled (i.e.baseline) L as: where , and ∆m 2 ji = m 2 j −m 2 i is the mass-squared difference of the mass eigenstates ν j and ν i .The survival probability of ν e depends only on three mixing angles θ 12 , θ 13 and θ 14 and six masssquared differences, only three of which are independent.According to Eq. ( 2) the survival probability oscillates with wavelength proportional to E/∆m 2 .Assuming m 1 to be the lightest and ∆m 2  21 ≪ ∆m2 41 , the survival probability may be approximated for baselines below hundreds of meters and for MeV-scale energies as1 : P νe νe ≈ 1 − sin 2 2θ 14 (cos 2 θ 13 sin 2 ∆ 41 + sin 2 θ 13 sin 2 ∆ 43 ) − cos 4 θ 14 sin 2 2θ 13 sin 2 ∆ 32 .
(3) One possible hint of a sterile neutrino would thus be a deficit of the electron-antineutrino event rate, accompanied by the dependence of the energy-spectrum distortion pattern on the distance between source and detector.
Such a deficit would be visible in the Daya Bay Reactor Neutrino Experiment with baselines spanning from 360 m to 1900 m [16] that can provide world-leading sensitivity to a sterile neutrino with ∆m 2 41 < 0.2 eV 2 .This Letter reports the search for such a state with the Daya Bay experiment's full data sample of inverse beta decay (IBD) interactions identified by subsequent neutron-capture on gadolinium.
From 2011 to 2020, the Daya Bay experiment operated with up to eight identically designed antineutrino detectors (ADs) near the Daya Bay and Ling Ao nuclear power plants in southern China [17,18].The ADs were distributed among two Near experimental halls (EH1 and EH2, each with up to 2 ADs) and one Far hall (EH3, with up to 4 ADs).Each AD contained a target mass of ∼20 tons of gadolinium-doped liquid scintillator, observed by 192 8-inch photomultiplier tubes (PMTs).
As usual for reactor-based experiments, the Daya Bay experiment detects electron antineutrinos via the IBD interaction: νe + p → e + + n.The outgoing positron ionizes the scintillator and annihilates with an electron, producing a signal corresponding to a "prompt" energy with E p ≈ E ν − 0.8 MeV.The neutron is captured with an average delay of 28 µs, predominantly by gadolinium, which deexcites via a cascade of gamma rays totaling ∼8 MeV, forming the "delayed" signal.In this channel, the close temporal coincidence of the two signals, together with the high energy of the delayed signal, allow for an average background-to-signal ratio of ∼ 1.5%.
The Daya Bay experiment operated with different detector configurations in three consecutive periods, which we label the 6-AD, 8-AD, and 7-AD periods, based on the number of active ADs.In this analysis, we adopt the same IBD data sample used in the most recent three-flavor neutrino oscillation analysis [18].This full dataset, with a total of ∼5.55 million IBD candidates, has twice the statistics of the sample used in the previous sterile neutrino search [19].The daily IBD rates and the estimated backgrounds in the three halls are summarized in Table I for the three periods.In addition to the increased statistics, this analysis benefits from improvements to four systematic uncertainties: i) The 9 Li/ 8 He background was estimated using a new multi-dimensional fitting method [18]; ii) The effect of spent nuclear fuel was derived from detailed history on reactor operation [20,21]; iii) The channel-wise electronics nonlinearity was recalibrated using an FADC readout system [21,22]; iv) The energy response model was constrained with new calibration data [21,22].
The background rate is dominated by accidental coincidences of uncorrelated signals satisfying the IBD selection criteria.To mitigate this background, in the summer of 2012, the 241 Am-13 C neutron sources were removed from two of the automated calibration units on each far AD, halving the rate of delayed-like uncorrelated signals in the far hall [16].Although these accidentals dominate the background rate, the uncertainty of the total background rate is dominated by the contribution from correlated pairs induced primarily by cosmic-ray muons.In particular, cosmogenic 9 Li/ 8 He is a well-known background in liquid scintillation detectors used for reactor neutrino experiments, and in the Daya Bay experiment it is the leading contributor to the background uncertainty.In this analysis, the relative uncertainty of the estimated 9 Li/ 8 He background rate has been reduced from ∼35% [16] to <25% by taking into account the correlated temporal and energy information of the IBD candidates [18].
The antineutrino flux is predicted from thermal-power data and from the calculated fission fractions of each fuel cycle.Uncertainties in the thermal-power data lead to a core-to-core uncorrelated flux uncertainty of 0.5%, while a 0.6% uncorrelated uncertainty per core in the νe yield was introduced by the uncertainties of the fission fractions.The spent nuclear fuel in the cooling pools adjacent to each core contributes 0.3 ± 0.1% to the predicted neutrino flux [20,23].The correlated flux uncertainty includes contributions from the 0.2% uncertainties of the mean energy released per fission [24].However, the primary contributors to the correlated uncertainty are the theoretical uncertainties on forbidden decays and the missing information in published nuclear data tables [25], which lead to a conservative total of 5% for the core-tocore correlated flux uncertainty.The size of the reactor cores and ADs has negligible impact on the sterile neutrino sensitivity due to the relatively long baselines.
The nominal reactor antineutrino spectra from U 235 , Pu were predicted using the models of Huber [26] and Mueller et al. [27].Multiple methods [16] were used to account for the known disagreements between these models and existing measurements [28].One method applied enlarged uncertainties, ranging from 10% to 40% as a function of energy, to the neutrino spectra.An alternative method used the observations of the near detectors to predict the observations of the far detectors, while another used free parameters to modify the predicted antineutrino spectra in the fit.As an additional cross-check, the analysis was also repeated using summation spectra [29].All of the aforementioned methods produce consistent results.
To predict the IBD rate and reconstructed promptenergy spectrum at each AD, the reactor spectra were multiplied by the neutrino oscillation probability (Eq.( 2)), the IBD differential cross section, the number of target protons, the detection efficiency, and the detector's energy response, which maps positron energy to reconstructed prompt energy.The energy response model considers energy resolution and nonlinearity and possible energy loss in the inner acrylic vessel.
To further improve the energy nonlinearity model, in December 2015 a full flash ADC (FADC) readout system was installed in EH1-AD1, recording PMT waveforms at 1 GHz and 10-bit resolution [22].A deconvolution method was applied to the waveforms to minimize any dependence on the single-photoelectron pulse shape (in particular the overshoot) and to extract the integrated charge with minimal bias.The residual nonlinearity in the corrected charge from a single waveform was estimated to be less than 1%, resulting in a 0.2% nonlinearity in the total charge measurement for each event [22].In addition, a special calibration campaign in January 2017 improved the knowledge of optical shadowing by the radioactive source enclosures, reducing the energy calibration uncertainties of γ rays from 1% to 0.5% [22].This combination of improvements, when applied to calibration data, led to a significant improvement in the energy response model for positrons.The uncertainty improved from ∼ 1% [16,19] to < 0.5% [21,22] for E p > 2 MeV.
In the region of ∆m 2 41 where the Daya Bay experiment provides world-leading sensitivity on sin 2 2θ 14 , the bulk of the sensitivity comes from relative measurements between the experimental halls.Figure 1 shows the double ratio of the measured to three-neutrino-predicted spectra of EH2 and EH3 to that of EH1.The data are compared to two four-neutrino predictions.The data are well contained in the uncertainty band, indicating that the data are consistent with the three-neutrino predic-

tion.
To quantify our measurement, several different statistical methods have been used in this analysis, yielding consistent results.One method used a χ 2 statistic based on the binned maximum likelihood ratio with systematic uncertainties treated via Gaussian nuisance terms.An alternative method built a χ 2 function using a covariance matrix generated by a toy Monte Carlo that includes fluctuations due to systematic uncertainties.Hybrids of the two approaches have also been tested.The free parameters are sin 2 2θ 13 , sin 2 2θ 14 and ∆m 2 41 .We used ∆m 2 32 = (2.453± 0.034) × 10 −3 eV 2 [14] (normal mass ordering).
To test the consistency between the Daya Bay experimental data and the three-(3ν) or four-neutrino (4ν) hypotheses, we first calculated χ 2 (η) at each point in the ≡ sin 2 2θ 14 , ∆m 2 41 parameter space by profiling over sin 2 2θ 13 and the nuisance parameters.These parameters are those which minimize χ 2 at each point in η.We then defined a test statistic ∆χ 2 = χ 2 (η (0,0) ) − χ 2 (η bf ) with 2 degrees of freedom sin 2 2θ 14 , ∆m 2 41 , where η (0,0) represents the 3ν oscillation assumption and η bf represents the global best fit under the assumption of 4ν oscillation.The Daya Bay experimental data gave ∆χ 2 = 2.3, corresponding to a p-value [14] of 0.86, obtained from the ∆χ 2 distribution generated by Monte Carlo samples under the 3ν oscillation hypothesis including statistical and systematic variations.This indicates that no significant signal of a sterile neutrino was observed.
The Feldman-Cousins (FC) method [30] was used to set confidence intervals in the η parameter space.For each point η, a distribution of ∆χ 2 = χ 2 (η) − χ 2 (η bf ) was generated from 1000 pseudoexperiments with both statistical and systematic variations considered.Based on the ∆χ 2 distribution and the ∆χ 2 observed with the data, a p-value for each η point was calculated.The 1 − α confidence interval boundary was set where p-value=α.
An alternative method for determining limits is the Gaussian CL s method [31] based on a two-hypothesis test, comparing the null hypothesis (3ν) and the alternative hypothesis (4ν) for each point η.Using the resulting test statistic ∆χ 2 = χ 2 (η) − χ 2 (η (0,0) ), the corresponding CL s value was calculated [14]: where p 4ν and p 3ν are the p-values of the two hypothesis.
The condition of CL s < α was used to set the (1 − α) CL s exclusion region [14].
The results of the application of the FC and the CL s methods are shown in Figure 2. The acronym "95% C.L." represents both the 95% confidence interval for the FC method and the 95% exclusion region for the CL s method.Both contours show consistent features.The CL s method provides more stringent limits due to the different definition of the test statistic ∆χ 2 .The decrease of sensitivity in the region of ∆m 2 41 ≈ ∆m 2 32 is related to the fact that the oscillations to the sterile neutrino state are not distinguishable from 3ν oscillations.At the baseline of the Daya Bay experiment, the choice of neutrino mass ordering has a marginal impact on the results.The best sensitivity to sin 2 2θ 14 is achieved in the region 10 −2 eV 2 ≲ ∆m 2 41 ≲ 0.1 eV 2 where the measurement relies on the relative spectral differences between the detectors of the Near and Far experimental halls.For the higher mass-squared difference, the sensitivity decreases as the detectors become insensitive to the shape of sterileactive oscillations.The contour tends to a constant in the region |∆m 2 41 | ≳ 0.5 eV 2 where the prediction becomes limited by the uncertainty in the reactor antineutrino flux.To demonstrate the effects of different types of uncertainties [16], the sensitivity with the CL s method under various scenarios is shown in Figure 3.The main effect of the reactor antineutrino flux uncertainties falls in the region of |∆m 2 41 | ≳ 4 × 10 −3 eV 2 , where the sensitivity is dominated by the relative spectral difference between the two Near EHs [15].The sensitivity for ∆m 2  41 ≲ 3 × 10 −2 eV 2 is affected by the uncertainties of the detector energy response model, where the relative difference between EH3 and the two Near EHs plays the most important role [15].Background uncertainties have a negligible effect on the sensitivity contour due to both the low background level and the accurate estimation of the background.
A comparison with the results of other short and medium baseline reactor neutrino experiments is shown in Figure 4.The Day Bay experiment is able to set the most stringent upper limits on light sterile-active neutrino mixing for 2 × 10 −4 eV 2 ≲ ∆m 2 41 ≲ 0.2 eV 2 due to its high statistics and well-controlled systematic uncertainties.Currently, short baseline (≲ 100 m) experiments like Bugey-3 and NEOS give more stringent limits than Daya Bay in the region of ∆m 2 41 ≳ 0.2 eV 2 .In the future, the JUNO experiment will dominate in the region of ∆m 2 41 ≲ 10 −3 eV 2 , because of its long baseline (∼ 53 km) and superb energy resolution [32].
In summary, we report the results of a light sterile neutrino search using the full dataset of the Daya Bay

ACKNOWLEDGMENT
The Daya Bay experiment is supported in part by the Ministry of Science and Technology of China, the U.S. Department of Energy, the Chinese Academy of Sciences (CAS), the National Natural Science Foundation of China, the New Cornerstone Science Foundation, the Guangdong provincial government, the Shenzhen municipal government, the China General Nuclear Power Group, the Research Grants Council of the Hong Kong Special Administrative Region of China, the Ministry of Education in Taiwan, the U.S. National Science Foundation, the Ministry of Education, Youth, and Sports of the Czech Republic, the Charles University Research Centre (UNCE), the Joint Institute of Nuclear Research in Dubna, Russia, and the National Commission of Scientific and Technological Research of Chile.We acknowledge Yellow River Engineering Consulting Co., Ltd., and China Railway 15th Bureau Group Co., Ltd., for building the underground laboratory.We are grateful for the with those from other reactor experiments.The methods used to generate these contours are not identical.The Daya Bay experiment (DYB) used the CLs method, while RENO used traditional method based on Wilk's theorem [33,34], and Double Chooz (D-Chooz) [35], Bugey-3 [36] and NEOS [37] all used the raster scan (RS) [38] method.The contours from applying the RS and traditional method to the Daya Bay data are shown in the supplemental material, and are consistent with the contour from the CLs method.
ongoing cooperation from the China Guangdong Nuclear Power Group and China Light & Power Company.

2 FIG. 1 .
FIG.1.Impact of sterile neutrino oscillation on the double ratios (R) of the spectra, obtained at EH2(3), to EH1 as R = (M EH2(3) /P 3ν EH2(3) )/(MEH1/P 3ν EH1 ), where M is the measured prompt energy spectrum, and P 3ν is the threeneutrino prediction with the best-fit parameters.The error bands represent the total uncertainties (statistical and systematic) of the three-neutrino prediction.The error bars represent only the statistical uncertainties.For the predictions, sin 2 2θ14 = 0.1 and two ∆m 2 41 values are shown.

FIG. 2 .
FIG.2.Exclusion contours at 95% C.L., obtained by the FC and CLs methods.For the FC method the 1σ and 2σ bands are also shown, which account for the statistical and all systematic uncertainties.The parameter space to the right of the contours is excluded.

FIG. 4 .
FIG. 4. Comparison of the Daya Bay experiment's resultswith those from other reactor experiments.The methods used to generate these contours are not identical.The Daya Bay experiment (DYB) used the CLs method, while RENO used traditional method based on Wilk's theorem[33,34], and Double Chooz (D-Chooz)[35], Bugey-3[36] and NEOS[37] all used the raster scan (RS)[38] method.The contours from applying the RS and traditional method to the Daya Bay data are shown in the supplemental material, and are consistent with the contour from the CLs method.

TABLE I .
Summary of IBD signal and background for the full dataset for 3 detector configurations and 3 experimental halls (EHs).The errors include statistical and systematic uncertainties.The effective live time is the product of live time, the efficiency of the muon veto and the efficiency of the multiplicity selection.AD/day] 12.36 ± 0.58 9.96 ± 0.59 3.16 ± 0.08 10.87 ± 0.57 8.26 ± 0.42 1.11 ± 0.03 10.61 ± 0.81 7.75 ± 0.41 0.94 ± 0.03