Improved Measurement of the Reactor Antineutrino Flux at Daya Bay

This work reports a precise measurement of the reactor antineutrino flux using 2.2 million inverse beta decay (IBD) events collected with the Daya Bay near detectors in 1230 days. The dominant uncertainty on the neutron detection efficiency is reduced by 56% with respect to the previous measurement through a comprehensive neutron calibration and detailed data and simulation analysis. The new average IBD yield is determined to be $(5.91\pm0.09)\times10^{-43}~\rm{cm}^2/\rm{fission}$ with total uncertainty improved by 29%. The corresponding mean fission fractions from the four main fission isotopes $^{235}$U, $^{238}$U, $^{239}$Pu, and $^{241}$Pu are 0.564, 0.076, 0.304, and 0.056, respectively. The ratio of measured to predicted antineutrino yield is found to be $0.952\pm0.014\pm0.023$ ($1.001\pm0.015\pm0.027$) for the Huber-Mueller (ILL-Vogel) model, where the first and second uncertainty are experimental and theoretical model uncertainty, respectively. This measurement confirms the discrepancy between the world average of reactor antineutrino flux and the Huber-Mueller model.


I. INTRODUCTION
Nuclear reactors are an intense man-made source of electron antineutrinos and were used for the first observation of the neutrino [1]. Electron antineutrinos can be detected through inverse beta decay (IBD) on target protons, where a prompt positron and a delayed neutron capture signals are measured in time coincidence. Since the early 2000s, the energy and baseline (the distance between source and detector) dependent neutrino disappearance at nuclear reactors [2][3][4] has provided strong evidence of neutrino oscillation [5][6][7]. However, a recent re-evaluation of the theoretical prediction (referred to as Huber-Mueller model [8,9]) of the reactor neutrino flux resulted in a ∼6% deficit in measured flux from short-baseline experiments [10] and the previous ILL-Vogel model [11][12][13][14]. The difference between the data and Huber-Mueller prediction, i.e. the so-called "reactor antineutrino anomaly" (RAA), could be interpreted as active-to-sterile neutrino oscillation with a mass-squared splitting (∆m 2 ) around 1 eV 2 . It is also shown in Refs. [15][16][17] that the allowed parameter space is compatible with earlier anomalies from LSND [18,19], MiniBooNE [20], GALLEX [21], and SAGE [22]. On the other hand, a number of authors [23][24][25][26] have argued that the RAA may be due to the theoretical uncertainties in the flux calculations. Recent antineutrino flux evolution results from Daya Bay are in tension with the sterile-neutrino-only explanation of RAA [27].
The uncertainty of the reactor antineutrino flux in our previous measurement [27] is dominated by the uncertainty of neutron detection efficiency. The neutron detection efficiency was determined to be ε n = (81.83 ± 1.38)% [28,29], and the ratio with respect to the total uncertainty is σ 2 εn /σ 2 total = 65%. To further elucidate the RAA situation, this work presents an updated flux measurement from Daya Bay using the same 1230-day data set, but with a more precise determination of the neutron detection efficiency. Key improvements include an elaborated neutron calibration campaign covering a wide range of neutron energy and positions, an improved simulation with different physics models, and a data-driven correction to the neutron efficiency. This paper is organized as follows. In Sec. II, we explain the general method to measure reactor neutrino yield, and highlight our approach here to improve its estimate. Sec. III discusses the neutron calibration campaign and the analysis of calibration and simulation data. In Sec. IV we present an improved reactor antineutrino flux measurement and a comparison with the world data and theoretical models.

II. METHOD
A. Overview of procedure The Daya Bay experiment has four near and four far identically designed antineutrino detectors (ADs), located at different baselines (360 m-1900 m) [30] measuring the electron antineutrino flux from six reactor cores. The structure of the detector is shown in Figure 1. Each AD consists of a cylindrical target volume with 20 tons of 0.1% gadolinium loaded liquid scintillator (GdLS, 3.1 m in diameter and 3.1 m in height), surrounded by a layer of 42-cm thick liquid scintillator (LS) to enclose the gammas or electrons escaped from the central GdLS region. The GdLS and LS are separated by a 1-cm thick acrylic vessel. An energy deposit in the GdLS and LS regions is detected by photomultiplier tubes (PMTs). The origin of the coordinate system is set at the geometrical center of the GdLS cylinder, with the z-axis pointing up. IBD neutrons are detected by delayed capture either on hydrogen emitting one 2.2 MeV gamma or on gadolinium emitting several gammas with total energy of about 8 MeV. The kinetic energy of the IBD neutrons is less than 50 keV. The average capture time in the GdLS region is about 28.5 µs and 216 µs in the LS [31].
The reactor antineutrino IBD candidates are selected with the same criteria as in Ref. [32], which are also described here. 1) Removal of events caused by PMT light emission. 2) The time between the prompt and delayed signal is in the range of [1,200] µs. 3) Prompt signal must have a reconstructed energy, E, between 0.7 and 12 MeV. 4) Delayed signal must have E between 6 and 12 MeV to select gadolinium captures. 5) Muon anti-coincidence. 6) Multiplicity cut to remove events with E > 0.7 MeV in the interval 200 µs before the prompt signal, 200 µs after the delayed signal, or between the prompt and delayed signals.
The dominant backgrounds are accidental coincident events and cosmic-ray muon induced 9 Li/ 8 He, which are less than 2% of the signal IBD rate for the four near ADs. After statistical subtraction of background, the total number of IBD signals, N IBD , is 2.201 × 10 6 for the four near detectors.
To compare to the theoretical predictions, the reactor antineutrino yield σ f , defined as the number of antineutrinos times IBD cross-section per fission, can be calculated by solving the following equation: where the index d is for four near detectors, index r is for the six reactor cores, N P d is the number of target protons of detector d, ε IBD is the IBD detection efficiency, P rd sur is the mean neutrino survival probability from the reactor r to detector d, N f r is the predicted number of fissions of the rth reactor core, L rd is the distance from reactor r to detector d, and c SNF is a correction term for spent nuclear fuel. P rd sur is calculated by integrating the cross-section-weighted oscillation survival probability over theν e energy spectrum, using sin 2 2θ 13 and |∆m 2 ee | determined from the same data [32]. The average oscillation correction for near detectors is 1.5%.
The IBD detection efficiency is divided into two factors: where ε n is the neutron selection efficiency due to the [6,12] MeV cut and ε other is for the PMT light emission, prompt energy, and coincident time cuts. The predicted number of fissions of the rth reactor core is where W r is the thermal power of the rth core, E iso is the mean energy released per fission for each isotope, and f iso r is the average fission fraction of the rth core for each isotope, and the ratio is integrated over the live time of the detectors. The original thermal power and fission fuel composition data are provided by the power plant. c SNF was estimated to be (0.3 ± 0.3)% previously [29].
Different components of relative uncertainties for the antineutrino yield measurement from previous work [28,29], including ε n and ε other , are summarized in Table I. Clearly ε n dominates the uncertainty, and is the target of improvement in this paper.

B. Principle of improvement
The neutron detection efficiency, ε n , is composed of three individual factors.
• The Gd capture fraction is the fraction of neutrons produced by IBD in the GdLS target that are captured on Gd. The capture fraction is lower at the edge of GdLS volume because neutrons may drift into un-doped LS volume (spill-out effect).
• The nGd gamma detection efficiency is the fraction of neutron Gd capture signals with detected gamma energy above 6 MeV.
• Spill-in: efficiency increase due to IBD events produced in the LS and acrylic but with neutron capture on Gd.
We note that the estimation of Gd capture fraction and spill-in effects are strongly correlated since they are both driven by the modeling of neutron propagation including neutron scattering in materials and the subsequent nuclear capture. The estimation of the nGd gamma detection efficiency depends on the modeling of gamma emission including the multiplicity and energy spectrum of the emitted gammas.
In the previous study [29], we attempted to use different neutron calibration data to estimate these individual effects. The main difficulty was that no data can cleanly separate their uncertainties. In this paper, instead, we evaluate ε n and its uncertainty directly using a new neutron calibration data set and a data-simulation comparison. This approach is data driven and allows a significant reduction of the uncertainty.

A. Neutron calibration campaign
An extensive neutron calibration campaign was carried out in Daya Bay at the end of 2016. Two types of custom sources were fabricated, 241 Am-13 C (AmC, neutron rate ∼ 100 Hz) [33] and 241 Am-9 Be (AmBe, neutron rate ∼ 30 Hz). They produce neutrons through 13 C(α, n) 16 O or 9 Be(α, n) 12 C reactions with the final nucleus either in the ground state (GS) or excited state (ES). The kinetic energy of the neutrons from AmC (AmBe) in the GS and ES are [3,7] MeV ([6, 10] MeV) and <1 MeV ([2, 6] MeV), respectively. Calibration events are formed from the prompt energy of the proton recoil and deexcitation gammas, if 16 O * or 12 C * is created, and the delayed neutron capture. The high neutron rate and delayed-time-coincidence present a high signal-to-background ratio for the calibration study.
The sources, sealed in a small stainless steel cylinder (8 mm in both diameter and height), were enclosed in a highly reflective PTFE (Polytetrafluoroethylene) shell. These sources were deployed vertically into a near-site AD using the automated calibration units (ACUs) [34] along the central axis (ACU-A), an edge axis of GdLS at a radius of 1.35 m (ACU-B), and through a middle axis of the LS layer at a radius of 1.77 m (ACU-C). During deployment, the absolute precision of the source z location is 7 mm [34]. All calibration positions are illustrated in Figure 1. In total, data in 59 different source (final nucleus states) and location points (SLPs), were collected.
Delayed coincidence events for the calibration sample are selected with a time requirement of 1 µs < ∆t < 1200 µs for all events with E between [0. 3,12] MeV. The 1200 µs selection cut is set to efficiently include neutron captures in the LS and acrylic region. Two example distributions of the prompt-delay energies of AmC (ACU-B z=0 m) and AmBe (ACU-B z=0 m) samples are shown in Figure 2 and Figure 3, respectively. The data in different channels are selected using the following prompt energy cut: [0. 3,4] MeV for AmC GS, [5.5, 7] MeV for AmC ES, [1,4] MeV for AmBe GS, and [4.2, 7] MeV for AmBe ES. The accidental background contributes 0.1-20% of the neutron candidates depending on the SLP, and is estimated by randomly paired single events [31]. The reactor antineutrino and cosmogenic backgrounds are estimated by applying the same selection cuts on the data acquired immediately before and after the calibration campaign. They contribute <0.1% to the neutron source signals. All of these backgrounds are statistically subtracted.

B. Neutron and gamma modeling in simulation
The neutron calibration data were compared to the model predictions obtained using the Geant4-based [35] (v4.9.2) Daya Bay Monte Carlo (MC) simulation framework NuWa [3,30] with improvements to the calibration pipe geometry, detector energy response, and neutron transport modelling [36]. These modifications improve the agreement between the neutron calibration data [32] and simulation.
Neutrons lose energy through various scattering processes before capture on a nucleus. There are no scattering models in Geant4 for the Daya Bay scintillator (average hydrogen-to-carbon ratio CH 1.61∼1.64 ) or acrylic , and a polyethylene model (CH 2 , "poly") [35]. The latter two models are built based on the ENDF database [37] and are quite different from the free gas model. The total scattering cross-section as a function of energy for the three models is shown in Figure 4. To approximate the Daya Bay (scintillator, acrylic) material pair, five combinations (Table II)  including a) (water, free gas) b) (water, poly), c) (poly, poly), d) (poly, free gas), and e) (free gas, free gas).
For the neutron capture gamma energy and multiplicity distributions, four different models (Table II) were selected, including 1) a native Geant4 model, 2) a Geant4 model with the photon evaporation process, 3) a model based on the Nuclear Data Sheets by L. Groshev et al. [38], and 4) a model based on the measured single gamma distribution of nGd capture at Caltech [29]. The energy spectra of the deexcitation gammas of gadolinium-155 and 157 are shown in Figure 5(a) and Figure 5(b), respectively, for these models. The gamma model-3 has the hardest gamma spectra.
The 20 available combinations provided by the five neutron scattering model combinations (a-e) and the four gamma models (1-4) are used to estimate ε n . Model a-1 was used in the previous analyses [29,32].

C. Data and simulation comparison
For each calibration SLP, a ratio F is calculated.
where N ( [6,12] MeV) and N ([1.5, 12] MeV) are the numbers of events with reconstructed delayed energy in the range of [6,12] MeV and [1.5, 12] MeV, respectively. 1.5 MeV is chosen to include the hydrogen capture peak. F is very sensitive to the relative strength of H vs. Gd capture peaks and the containment of the 8 MeV of gamma energy from nGd capture, and therefore provides a crucial benchmark for the neutron and gamma simulation models. For the 59 calibration SLPs, a χ 2 is constructed to measure the overall difference between data and MC predictions,  where (F data − F MC ) is a vector with 59 elements of the difference of F between the data and MC, and V is the covariance matrix. For most of the calibration points, the statistical uncertainty is dominant, but for the points near a GdLS boundary (ACU-A top, ACU-A bottom, ACU-B top, or ACU-B bottom), the distance to the acrylic and LS volume is comparable to the neutron drift distance, so they share a large common uncertainty due to the source location z uncertainty. The χ 2 and ε n values for all models are shown in Table II. The eight combinations with either neutron model-e or gamma model-3 are discrepant (χ 2 > 300), and therefore are excluded. The remaining twelve (4×3) models agree with the data reasonably well with χ 2 in the range of 52.1-156 ("reasonable models"). The best model with minimum χ 2 is model b-1. In Figure 6, the delayed energy spectra at two boundary calibration locations from data are compared to models b-1 and e-1, where model b-1 shows a better agreement with data for F . The data and best MC F values and their differences are shown explicitly in Figure 7 for all sources and locations. The systematic variations among the twelve reasonable models are overlaid, where the full spread among them, maximum minus minimum, are plotted as the gray bars. The variation in F from 1 to 85% for the 59 data points is due to the differences in the local geometry and neutron kinetic energy and is well reproduced by simulation. For most points, the best MC model, b-1, reaches an agreement with data at the sub-percent level, and the residual difference is mostly smaller than the model spread. Another quantity, F = N ([3, 4.5] MeV)/N ( [3,12] MeV), is also constructed, which is also sensitive to the gamma model and energy leakage. The same data-model comparison procedure confirms that gamma model 3 should be rejected and that gamma model 1 is reasonable.
To further investigate potential effects due to the discreteness of the calibration sources and the energy difference between neutron sources and IBD neutrons, we exploited a large sample of IBDs from data as a special SLP to be compared to the model prediction. Due to the resolution in position reconstruction, selection of pure GdLS IBDs is impossible. Instead, a GdLS+LS IBD sample from all four near site ADs were selected using cuts identical to those used on the neutron calibration data (Sec. III-A), except that the prompt energy cut was adjusted to be greater than 3.5 MeV to suppress accidental background. About 2 million GdLS+LS IBD events in total were selected. The measured ratio F is consistent AD to AD with an average of 47.1%±0.1%. The ratio from model b-1 is 47.0%, and the full model spread is from 46.7% to 47.5%.

D. Neutron detection efficiency determination
Each model can give a prediction on ε n for IBDs (see Table II). For the twelve reasonable models, ε n ranges Delayed energy spectra of neutron sources at two calibration locations for data and simulation and the corresponding residual plots (MC-data). Two MC models, b-1 (best fit) and e-1 (rejected), are overlaid. Normalization is determined using the integral between 1.5 and 12 MeV. The difference of relative gadolinium/hydrogen capture ratio between sub-figure (a) and (b) is due to the relative position to the GdLS volume. In (b), a weak signal of neutron capture on carbon can be seen. In (b), a mismatch of the energy scales of data and MC at nH peak is observed, but our selection cut efficiency is not sensitive to the difference. from 81.61% (model c-1) to 82.55% (model a-4), and that from model b-1 is 81.75%. Instead of taking the prediction as is, one can translate the data and MC difference in F to a correction to ε n , since the two are intrinsically correlated (linear to the lowest order) through the neutron and gamma models mentioned above. In mathematical form, for the ith SLP, we have where ε MC,best is the neutron detection efficiency given by the best MC model. c i characterizes the linear correlation between F i and ε n , and can be estimated through a linear regression (fit) using predicted values of ε n and F i from all 20 MC models. This procedure is illustrated in Figure 8. The eight rejected models were also included here by default as larger variations in ε n and F are allowed. Excluding them from the fit does not change the result (See Figure 8 for an example). In addition to ε n determined from 59 individual SLPs, a multiple regression procedure was also applied to model the linear relation between ε n and F from a set of SLPs. Taking all values of corrected ε n into account, a shift of −0.27% with a standard deviation of 0.47% was obtained, relative to that from the best model (81.75%). The standard deviation of 0.47% is consistent with the model spread on ε n . Aside from the model uncertainties, other systematic  effects (e.g. gadolinium abundance, source geometry, absolute energy scale, and material density variations) have been studied in the MC and found to be negligible.

IV. ANTINEUTRINO YIELD AND COMPARISON WITH PREDICTION
Using the new neutron detection efficiency ε n and Eq. 2, the IBD detection efficiency ε IBD is (80.25 ± 0.61)%. Using the procedure as in Eq. 1, the mean IBD reaction yield per nuclear fission is σ f = (5.91 ± 0.09) × 10 −43 cm 2 /fission, where the major uncertainties (Table I) are from the target proton fraction 0.92% (relative uncertainty), dominated by the hydrogen-to-carbon ratio due to instrumental uncertainty in the combustion measurements, and reactor-related uncertainty 0.90% (relative uncertainty) due to reactor power and fission fractions.
The ratio of the yield to the prediction of the Huber-Mueller (or ILL-Vogel) reactor model can be calculated. The effective fission fractions for four fission isotopes are defined as where iso refers to one of the four major fission isotopes, i.e. 235 U, 238 U, 239 Pu, and 241 Pu, N f,iso r is the predicted number of fissions contributed by the iso th isotope in the rth reactor core, and other symbols are defined in Eq. 1. In the analyzed data, the effective fission fractions for the four fission isotopes ( 235 U, 238 U, 239 Pu, and 241 Pu) are determined to be (0.564, 0.076, 0.304, and 0.056), respectively. The predicted IBD yield is the sum due to all four isotopes, including corrections due to nonequilibrium effects, in which S iso (E ν ) is the predicted antineutrino spectrum for each isotope given by Huber-Mueller or ILL-Vogel model, σ IBD (E ν ) is the IBD cross section, and k NE iso (E ν ) corrects for the non-equilibrium long-lived isotopes. The calculation integrates over neutrino energy E ν and the non-equilibrium effect contributes +0.6% [29]. The ratio between the measured to predicted reactor antineutrino yield R is 0.952 ± 0.014 ± 0.023 (Huber-Mueller) and 1.001 ± 0.015 ± 0.027 (ILL-Vogel), where the first uncertainty is experimental and the second is due to the reactor models themselves. A breakdown of the experimental uncertainties can be seen in Table I (see also Ref. [29]). The uncertainties from power, spent fuel, and non-equilibrium are treated to be uncorrelated among different reactor cores in the oscillation analysis [32], and those from fission fraction, IBD cross section, and energy/fission are treated to be correlated. They are conservatively treated as fully correlated in this analysis, and the total reactor-related uncertainty is 0.9%. The total experimental uncertainty has been reduced to 1.5%, which is a relative 29% improvement on our previous study. The new flux measurement is consistent with the ILL-Vogel model, but differs by 1.8 standard deviations with respect to the Huber-Mueller model, with the uncertainty now dominated by the theoretical uncertainty. With the new result, a comparison with the other measurements is updated using the same method presented in Ref. [29]. A summary figure is shown in Figure 9. The Daya Bay new result on R is consistent with the world data. The new world average of R is 0.945 ± 0.007 (exp.) ± 0.023 (model) with respect to the Huber-Mueller model. This more precise measurement further indicates that the origin of RAA is unlikely to be due to detector effects.

V. SUMMARY
In summary, an improved antineutrino flux measurement is reported at Daya Bay with a 1230-day data set. The precision of the measured mean IBD yield is improved by 29% with a significantly improved neutron detection efficiency estimation. The new reactor antineutrino flux is σ f = (5.91±0.09)×10 −43 cm 2 /fission. The ratio with respect to predicted reactor antineutrino yield R is 0.952 ± 0.014 ± 0.023 (Huber-Mueller) and 1.001 ± 0.015 ± 0.027 (ILL-Vogel), where the first uncertainty is experimental and the second is due to the reactor models. This yield measurement is consistent with the world data, and further comfirms the discrepancy between the world reactor antineutrino flux and the Huber-Mueller model.