Precision measurement of reactor antineutrino oscillation at kilometer-scale baselines by Daya Bay

We present a new determination of the smallest neutrino mixing angle ${\theta}_{13}$ and the mass-squared difference ${\Delta}{\rm m}^{2}_{32}$ using a final sample of $5.55 \times 10^{6}$ inverse beta-decay (IBD) candidates with the final-state neutron captured on gadolinium. This sample was selected from the complete data set obtained by the Daya Bay reactor neutrino experiment in 3158 days of operation. Compared to the previous Daya Bay results, selection of IBD candidates has been optimized, energy calibration refined, and treatment of backgrounds further improved. The resulting oscillation parameters are ${\rm sin}^{2}2{\theta}_{13} = 0.0851 \pm 0.0024$, ${\Delta}{\rm m}^{2}_{32} = (2.466 \pm 0.060) \times 10^{-3}{\rm eV}^{2}$ for the normal mass ordering or ${\Delta}{\rm m}^{2}_{32} = -(2.571 \pm 0.060) \times 10^{-3} {\rm eV}^{2}$ for the inverted mass ordering.

Neutrino oscillation has been firmly established by multiple observations since its discovery in 1998 [1]. As this phenomenon is not required by the Standard Model, it offers opportunities to search for new interactions and physical principles. The three-neutrino paradigm of neutrino oscillation can be parametrized by three mixing angles, two mass-squared differences, and a CP phase [2]. This framework has been very successful in explaining most of the observations made with accelerator, atmospheric, reactor and solar neutrinos. Our knowledge of the smallest neutrino mixing angle θ 13 has been steadily improving since the first definitive determination in 2012 [3]. Besides being the best-measured neutrino mixing angle at present, precise knowledge of θ 13 is important for testing the three-neutrino paradigm of neutrino mixing and as an invaluable input to modelbuilding and to other experiments, most notably in resolving the neutrino mass hierarchy [4] and the search for CP violation in neutrino oscillation [5,6].
Nuclear reactors produce low-energy electron antineutrinos, ν e s, that are ideal for determining θ 13 and the mass-squared difference ∆m 2 32 through the study of ν e disappearance. This is best accomplished by comparing the energy spectra obtained with identically designed detectors positioned at different distances from the reactors. This relative approach cancels the uncertainties in the absolute detection efficiency that are correlated between detectors and heavily suppresses the effect of the uncertainty in the reactorν e flux determination, thus enabling precision measurement of the oscillation parameters. The ν e s are detected via the inverse beta-decay re-action (IBD), ν e +p → e + +n, with the kinetic-energy loss and annihilation of the positron giving rise to a promptenergy (E p ) signal, and the subsequent neutron capture to a delayed-energy (E d ) signal. The energy of the ν e , E ν , central to measurements of neutrino oscillation, is inferred from E p with E ν ≈ E p + 0.78 MeV.
In this Letter we report a new measurement of sin 2 2θ 13 and ∆m 2 32 using a final sample of 5.55 × 10 6 IBD candidates with the final-state neutron captured on gadolinium (n-Gd) acquired by the Daya Bay reactor neutrino experiment in 3158 days of operation.
We utilized up to eight antineutrino detectors (ADs) to detect ν e s emitted from three pairs of 2.9-GW th reactors at the Daya Bay-Ling Ao nuclear power facility in Shenzhen, China. The ADs were installed in three underground experimental halls, EH1, EH2, and EH3, having a flux-averaged baseline of about 500 m, 500 m, and 1650 m from the reactors, respectively. To suppress ambient radiation, the ADs were submerged in water pools. Each pool was optically divided to function as inner (IWS) and outer (OWS) water Cherenkov detectors for detecting cosmic-ray muons. Four layers of Resistive Plate Chambers (RPCs) covering the top of each water pool provided another independent muon detector. IBD events were detected with 20 tonnes of liquid scintillator doped with 0.1% gadolinium by weight (GdLS) in each AD [7][8][9]. The GdLS was contained in a 3-m-diameter acrylic cylinder enclosed inside a 4-mdiameter acrylic cylinder filled with 22 tonnes of undoped liquid scintillator (LS). Optical photons generated in the scintillator were detected with 192 photomultiplier tubes (PMTs) covering the barrel surface of the AD [10]. The PMTs were arranged in 8 horizontal rings and 24 vertical columns. Highly reflective disks sandwiching the 4-m acrylic vessel were used to enhance the detection efficiency of scintillation photons. Radioactive sources and LEDs were stored in three automatic calibration units (ACUs) on top of each AD [11,12]. Detailed information of the experiment can be found in Refs. [13,14]. For each AD, a cylindrical coordinate system with the vertical zaxis being the symmetry axis and z = 0 at the AD center was used.
The Daya Bay experiment was operated with three different configurations of ADs in the three EHs. From 24 December 2011 to 28 July 2012 (217 days), the experiment ran in an initial six-AD configuration with 2 ADs in EH1, 1 AD in EH2 and 3 ADs in EH3 that resulted in the first observation of ν e disappearance at O(1 km) baselines [3]. An AD was added to both EH2 and EH3 during the summer of 2012 and this eight-AD configuration was operated from 19 October 2012 until 20 December 2016 (1524 days). Seven-AD operation occurred from 26 January 2017 until 12 December 2020 (1417 days) with one AD in EH1 re-purposed for liquid scintillator R&D for the JUNO experiment [15].
The results presented in this Letter are based on the data collected in the three configurations. Throughout the entire data analysis process, multiple groups within the collaboration provided validation and cross-checks.
Details of the analysis process and techniques can be found in Refs. [16,17]. In this Letter we focus on the improvements to the analysis techniques.
Accurate and precise measurement of the prompt energy E p is essential for extracting the oscillation parameters from the spectra. After the gain of each PMT was calibrated with the single-photoelectron peak from dark noise, a correction for the non-linear response of the electronics was applied to each channel. This correction was derived from the waveform output from a flash-ADC readout system running in parallel with the default ADC system of EH1-AD1 in 2016 [18]. The observed charge profile was then used to reconstruct the position of the event using Reconstruction B in Ref. [16].
To obtain the reconstructed energy (E rec ), an additional correction to the non-uniform detector response was applied to account for a few non-functional PMTs toward the end of data collection. We used the energy deposited by spallation neutron capture on Gd in the GdLS and delayed α-particles from correlated decays of natural radioactivity, 214 Bi→ 214 Po→ 210 Pb, in the LS to determine this additional position-dependent correction. The active volume of each AD was divided into 100 voxels in z and r 2 , where r is the radial distance from the z-axis. For each voxel, the correction was defined as the ratio of the reconstructed energy to the reconstructed energy averaged over the entire GdLS volume. The temporal dependence of this correction was accommodated by two calibration periods, before and after 31 March 2017. The largest additional per-voxel correction was about 3%.
In this study, the prompt energy was obtained by directly correcting E rec for the non-linear response of the LS which was determined from calibration [19]. Weekly calibration was performed by remotely lowering the calibration sources into the ADs from the ACUs. Specialized calibration runs were taken during the re-configuration periods [20]. The positron response model of Ref. [19] was updated, taking into account the measured responses of γ-rays from various sources and electrons from β-decay of cosmogenic 12 B of the full dataset as inputs. The best-fit model had a Birks' coefficient k B = 0.0143 g/cm 2 /MeV for the quenching effect and k C = 0.023 for the contribution of Cherenkov radiation to the non-linearity; both parameters agreed well with the previous result [19]. The improved energy response model for the positron achieved a precision of < 0.5% for E p > 2 MeV.
IBD candidates were selected with the following criteria. Events caused by spontaneous light emission of the PMTs, so-called flashers, were removed. Candidates must have a prompt-like signal with 0.7 MeV < E p < 12 MeV separated by 1 to 200 µs from a delayed-like signal with 6 MeV < E d < 12 MeV. Candidate pairs were vetoed if their delayed-like events occurred (i) within a (-202 µs, 600 µs) time-window with respect to an IWS or OWS trigger with a PMT-hit multiplicity (nHit) > 12, or (ii) within a (-202 µs, 410 µs) time-window with respect to an IWS trigger with 6 < nHit ≤ 12, or (iii) within a (-202 µs, 1400 µs) time-window with respect to triggers in the same AD with energy between 20 MeV and 2 GeV or (iv) within a (-202 µs, 0.4 s) time-window with respect to triggers in the same AD with energy higher than 2 GeV. This targeted muon veto efficiently removed spurious triggers that followed a muon as well as most muon-induced spallation products and muon decays. To remove any ambiguity in the candidate-pair selection, no additional AD triggers with energy between 0.7 MeV and 20 MeV were allowed within (-400 µs, 200 µs) of the delayed candidate.
A new source of flashers was observed in the 7-AD operation period that were not suppressed by the previous criteria [16]. Additional selection criteria targeting the characteristic charge pattern and temporal distribution of these new flashers were devised that rejected over 99% of this instrumental background with an IBD selection efficiency over 99.99%.
The selected IBD candidates consisted of genuine IBD and background events. The background comprised uncorrelated accidental pairs, and correlated prompt-anddelayed signals coming from fast neutrons, β-n decays of spallation 9 Li/ 8 He, neutrons leaking from the 241 Am-13 C calibration sources and 13 C(α, n) 16 O with the α coming from natural radioactivity. The latter two correlated backgrounds and the accidental background, detailed in Ref. [16], did not require any improved treatment in this analysis. The muon detection efficiency of the IWS and OWS dropped with time due to the gradual loss of functional PMTs near the top of the water pools, particularly in the 7-AD period. With this loss of detection efficiency, a new background, dubbed "muon-x" (described below), became apparent.
The largest correlated background is β-n decay of cosmogenic radio-isotopes 9 Li and 8 He. To determine this background, muons were paired with all IBD candidates within ±2 seconds. To improve discrimination of 9 Li/ 8 He from other processes, candidate events were separated into several samples based on the visible energy deposited by the muon (E µ ) in the AD and the distance between the prompt and delayed signals, ∆r. The rates and energy spectra of the dominant cosmogenic radio-isotopes were extracted with a simultaneous fit to 12 two-dimensional histograms defined by the different muon samples in the three experimental halls for the two ∆r regions with a probability density function φ(E p , ∆t), where ∆t is the time difference of the prompt-energy signal and the muon [21]. The distribution in ∆t is described by a sum over radio-isotopes, taking into account the known isotope lifetimes, and a term for uncorrelated muon-IBD pairs that is well-constrained by the pairs in the region ∆t < 0 s. Since the lifetimes of 9 Li and 8 He are comparable, we simply measured the sum of these two radioisotopes. The E p distributions of the radio-isotopes were determined from the fit. This method provides higher statistics and a better determination of the low-energy part of the β spectrum of 9 Li/ 8 He than the previous determination [16] while reducing the rate uncertainty to less than 25%.
We determined the combined contribution of the fastneutron and muon-x processes to the background.
Energetic neutrons generated by cosmic-ray muon interactions in the vicinity of the water pool can enter the active volume of an AD. Proton recoil from neutron scattering in the LS and the subsequent neutron capture on gadolinium constituted the prompt and delayed signals of this fast-neutron background and dominates the correlated events with E p > 12 MeV. The energy spectrum was determined from prompt signals in coincidence with muons detected only by the OWS or RPC within 0.5 µs of a prompt candidate with 0.7 MeV < E p < 250 MeV. The measured spectra were similar in shape among the three experimental halls and stable throughout all three operational periods.
The muon-x background was caused by low-E µ muons that passed through the IWS undetected. These events typically consisted of the muon as the prompt signal and a Michel electron from muon decay or a spallation neutron as the delayed signal. Due to the decrease in efficiency of the IWS that occurred mainly during the 7-AD period, this background could not be removed with the IWS nHit > 12 criterion, but was efficiently (> 80%) removed by vetoing events with a delayed signal less than 410 µs after a muon identified with an IWS nHit satisfying 6 < nHit ≤ 12, which led to a loss in livetime of less than 0.1%. The muon-x prompt-and delayed-energy spectra of the remaining muon-x background were approximated by the muon-x sample obtained with IWS nHit = 7.
To determine the rates of these two backgrounds, the prompt-energy spectrum above 12 MeV and the delayedenergy spectrum of the IBD-candidate sample with the prompt energy extended to 250 MeV were fitted to the spectra of the previously described OWS-tagged fastneutron sample and the muon-x sample with IWS nHit = 7. Their rates in the range of 0.7 MeV < E p < 12 MeV were estimated by extrapolation. We found no muon-x background in the 6-AD period. In the 8-AD period, the rate of muon-x background was about 0.038 of that of the fast neutron in EH1, 0.0055 in EH2, and 0.023 in EH3. The ratios increased to 0.26, 0.15, and 0.53, respectively, in the 7-AD period, consistent with the reduction in the number of functional PMTs in the three IWSs. The combined systematic uncertainty of these two backgrounds was estimated to be about 20%, which is dominated by the fast neutron background. As validation, the muon-x background was deduced by comparing the prompt-and delayed-energy spectra of data samples before the 7-AD period with and without masking the PMTs that failed subsequently. Consistent results were obtained. Table I summarizes the IBD candidates and backgrounds for the final n-Gd sample. We obtained a total of 4.8 million IBD candidates at the near halls and 0.76 million at the far hall with less than 2% background.
The ν e flux without oscillation at each AD was predicted by using the thermal-power data and fission fractions of each fuel cycle, provided by the power plant operator, as a function of burn-up. The power data carried an uncorrelated uncertainty of 0.5% per core, while a 0.6% uncorrelated uncertainty per core in the ν e yield was introduced by the uncertainties of the fission fractions. Due to the nature of the near-far relative measurement, 95% of the uncorrelated uncertainty of each core cancelled and extraction of the oscillation parameters was insensitive to the spectral shape of the no-oscillation prediction.
The detector-related uncertainties have been presented in Ref. [16]. Detection efficiency uncertainties that are correlated between detectors did not contribute to this near-far relative measurement. The total uncorrelated uncertainty in the detection efficiency remained at 0.13%. The largest contribution of 0.10% coming from the fraction of neutrons captured on gadolinium was obtained by comparing the capture-time distributions of the ADs. The next largest uncorrelated uncertainty of 0.08% in the delayed-energy selection criterion was due to a 0.2% spread in the relative energy scale among the ADs. The relative detection efficiency estimate was validated by comparing the ν e rates of neighboring ADs in each EH for each data-taking period [22]. The rates were consis-TABLE I. Summary of IBD signal and background. Rates are corrected for the muon veto and multiplicity selection efficiencies εµ × εm. The sum of the fast neutron and muon-x background rates is reported as "Fast n + muon-x". The AD numbering scheme reflects the time order of AD fabrication and deployment.

EH1
EH2 EH3  AD1  AD2  AD3  AD8  AD4  AD5  AD6  AD7  νe candidates  794335  1442475  1328301  1216593  194949  195369  193334  tent with the predictions that took the tiny variations in the baseline and number of protons into account. Furthermore, no significant deviation in the spectral distributions among the ADs in the same experimental hall was found. We extracted the oscillation parameters using the survival probability of three-flavor oscillation given by P = 1 − cos 4 θ 13 sin 2 2θ 12 sin 2 ∆ 21 −sin 2 2θ 13 cos 2 θ 12 sin 2 ∆ 31 + sin 2 θ 12 sin 2 ∆ 32 (1) where ∆ ij = 1.267∆m 2 ij L/E with ∆m 2 ij in eV 2 , L is the baseline in meters between an AD and a reactor core and E is the energy of the ν e in MeV. We used sin 2 θ 12 = 0.307 ± 0.013 and ∆m 2 21 = (7.53 ± 0.18) × 10 −5 eV 2 [2]. Alternatively, for short baselines of a few kilometers, the survival probability can be parametrized as P = 1 − cos 4 θ 13 sin 2 2θ 12 sin 2 ∆ 21 − sin 2 2θ 13 sin 2 ∆ ee . (2) Here, the effective mass-squared difference ∆m 2 ee is related to the wavelength of the oscillation observed at Daya Bay, and is independent of the choice of neutrino mass ordering as well as the value and uncertainty of the mixing angle θ 12 [16]. We adopted fitting Method B reported in Ref. [16] to extract the oscillation parameters. The fit minimized a χ 2 function defined as [21]: where χ 2 stat is the standard statistical term that compares all the measured background-subtracted prompt-energy spectra with the predictions. For each period of operation, the spectrum of each AD was divided into 26 bins. The predictions were derived from the calculated reactor ν e flux, survival probability, IBD cross section [23] and detector response obtained with a detailed Geant4-based simulation [24][25][26]. The term χ 2 syst (ν) contains the detector and background systematic uncertainties as pulls of the nuisance parameters expressed as a vector ν. Figure 1 shows the covariance contours in the ∆m 2 eesin 2 2θ 13 space. The best-fit point with χ 2 /ndf = 559/517 yields sin 2 2θ 13 = 0.0851 ± 0.0024, and ∆m 2 32 = (2.466 ± 0.060) × 10 −3 eV 2 for the normal mass hierarchy or ∆m 2 32 = −(2.571 ± 0.060) × 10 −3 eV 2 for the inverted mass hierarchy. Using Eq. 2, we obtained sin 2 2θ 13 = 0.0852 ± 0.0024 and ∆m 2 ee = (2.519 ± 0.060) × 10 −3 eV 2 with the same reduced-χ 2 value. Results determined with the other fitting methods described in Ref. [16] were consistent to < 0.2 standard deviations. Error ellipses in the ∆m 2 ee -sin 2 2θ13 space with the best-fit point indicated. The error bars display the onedimensional one-standard deviation confidence intervals. The colored contours correspond to one, two, and three standard deviations. The ∆χ 2 distributions are also shown. These onedimensional distributions were obtained by determining the smallest ∆χ 2 value after scanning through ∆m 2 ee (sin 2 2θ13 ) for a given sin 2 2θ13 (∆m 2 ee ).
The best-fit prompt-energy distribution is in excellent agreement with the observed spectra in each experimental hall, as shown in Fig. 2. Figure 3 depicts the normalized signal rate of the three halls as a function of L ef f / E νe with the best-fit curve superimposed, where L ef f and E νe are the effective baseline and average ν e energy, respectively [16]. The oscillation pattern related to θ 13 is unambiguous. The present improved result in sin 2 2θ 13 is consistent with our previous determinations [3,16,17] and agrees with other measurements of reactor ν e disappearance by RENO [27] and Double Chooz [28,29] as well as electron neutrino and antineutrino appearance measurements by T2K [6]. Daya Bay's measured ∆m 2 32 is consistent with the results of NOvA [5], T2K [6], MINOS/MINOS+ [30], IceCube [31] and SuperK [32] that were obtained with muon (anti)neutrino disappearance. The agreement in sin 2 2θ 13 and ∆m 2 32 between Daya Bay measurements using ν e and the muon neutrino and antineutrino determinations provides strong support of the three-neutrino paradigm.
To conclude, we have presented a new determination of sin 2 2θ 13 with a precision of 2.8% and the mass-squared differences reaching a precision of about 2.4%. The reported sin 2 2θ 13 will likely remain the most precise measurement of θ 13 in the foreseeable future and be crucial to the investigation of the mass hierarchy and CP violation in neutrino oscillation.