Improved Measurement of the Evolution of the Reactor Antineutrino Flux and Spectrum at Daya Bay

Reactor neutrino experiments play a crucial role in advancing our knowledge of neutrinos. A precise measurement of reactor electron antineutrino flux and spectrum evolution can be key inputs in improving the knowledge of neutrino mass and mixing as well as reactor nuclear physics and searching for physics beyond the standard model. In this work, the evolution of the flux and spectrum as a function of the reactor isotopic content is reported in terms of the inverse-beta-decay yield at Daya Bay with 1958 days of data and improved systematic uncertainties. These measurements are compared with two signature model predictions: the Huber-Mueller model based on the conversion method and the SM2018 model based on the summation method. The measured average flux and spectrum, as well as their evolution with the 239Pu fraction, are inconsistent with the predictions of the Huber-Mueller model. In contrast, the SM2018 model is shown to agree with the average flux and its evolution but fails to describe the energy spectrum. Altering the predicted IBD spectrum from 239Pu does not improve the agreement with the measurement for either model. The models can be brought into better agreement with the measurements if either the predicted spectrum due to 235U is changed or the predicted 235U, 238U, 239Pu, and 241Pu spectra are changed in equal measure.

The detection of reactor electron antineutrinos with the inverse-beta-decay (IBD) process plays a crucial role in advancing our knowledge of neutrinos including the discovery of neutrinos [1], establishment of large mixing angle solution of neutrino oscillation [2], and the discovery of non-zero mixing angle θ 13 [3].Looking forward, the JUNO experiment requires an accurate knowledge of the reactor neutrino spectrum to determine the neutrino mass ordering [4].
For commercial reactors, uranium isotopes are introduced at beginning of a fueling cycle and plutonium isotopes are gradually generated.Four fission isotopes 235 U, 238 U, 239 Pu, and 241 Pu account for over the 99.7% of the antineutrino flux with energy above the IBD detection threshold [5].A reactor antineutrino prediction, the Huber-Mueller (HM) model [6,7], is determined by converting cumulative beta spectra to antineutrino spectra for 235 U, 239 Pu, and 241 Pu and by summing all involved beta decay branches in databases for 238 U.The average of reactor neutrino flux measurements is only 95%-96% of the HM prediction, known as the reactor antineutrino anomaly (RAA) [8][9][10][11].Another anomaly is about the spectrum.The measured neutrino spectrum is poorly described by the HM model, e.g. a notable "bump" around 5 MeV [12][13][14].
Another prediction approach is the summation method, which adds up all related decay branches from databases for all four isotopes.One such example, the SM2018 calculation [34], with the latest experimental inputs, predicted a uniformly lower flux from 235 U than the HM model.Kopeikin et al. [35] reported the measured ratio between cumulative β spectra from 235 U and 239 Pu that is also systematically lower than the HM prediction.Both SM2018 and Kopeikin imply a much smaller discrepancy with neutrino flux measurement than HM.
The most recent results from Daya Bay on the total flux in terms of IBD yield, i.e., the number of antineutrinos per fission multiplied by the IBD cross section [9] and evolution of the spectrum as a function of reactor burnup used a 1230-day data sample [36].These results showed that the 235 U yield is about 8% less than the HM prediction while the 239 Pu yield is consistent with the model.The latest total and energy differential yields from 235 U and 239 Pu with a 1958-day data sample are reported in Ref. [37].Evolution studies have been performed for the NEOS [38] and RENO [39] experiments.
In this Letter, using the 1958-day data sample taken from December 2011 to August 2017 with the Daya Bay experiment [40], we report the direct measurement of the total and energy differential IBD yields, σ and σ e , and their evolution with reactor status with improved systematic uncertainties.Compared to the unfolded spectra of 235 U and 239 Pu [37], the measurements in this work do not introduce extra uncertainties from the unfolding method and the theoretical uncertainty of 238 U and 241 Pu which allows a more powerful examination of the combined reactor flux and spectrum prediction of the HM and SM2018 models.
The Daya Bay experiment, equipped with eight antineutrino detectors (ADs), measures the electron antineutrinos from six commercial reactors [41][42][43].The results in this Letter are based on approximately 3.5 million IBD candidates detected with the four near-site ADs.
The IBD process, νe + p → e + + n, is identified by the prompt-delayed coincidence.The delayed signal corresponds to the neutron captured on gadolinium.The prompt energy, E p , including the kinetic energy of positron and its annihilation gammas, is related to the antineutrino energy E ν ≈ E p + 0.78 MeV.The true E p deposit is reconstructed as E rec .The reconstructed energy resolution is about 8% at 1 MeV and a detector response matrix M (E rec , E ν ) is constructed taking into account all detector effects [40].The measured energy spectrum is corrected for the spent-nuclear-fuel contribution and the nonequilibrium contribution [36,37] for each AD and week, instead of being treated as time independent in the previous analysis [36].
To measure the IBD yield, a quantity N dw i is calculated for the d th AD and w th week, and i is 5, 8, 9, and 1 for 235 U, 238 U, 239 Pu, and 241 Pu, respectively [5].It describes the number of fissions of an isotope detected by an AD, and the definition is where N Proton d is the number of target protons of the d th AD, P sur dwr is the average survival probability of reactor electron antineutrinos integrated over energy from the r th reactor to the d th AD calculated under 3-active-neutrino framework in the w th week, ε is the detection efficiency, L dr is the distance of the AD-reactor pair, W rw is the thermal power of the r th reactor for the w th week, which is provided by the reactor company, T dw is the running time of that AD in that week, f irw is the fission fraction of the i th isotope in the r th reactor and w th week, and e i is the energy per fission of the isotope [44].The effective fission fraction for the i th isotope, F i (F 5 , F 8 , F 9 , and F 1 ), for that AD and week, F dw i , is defined by Data are sorted into 13 groups according to their effective 239 Pu fission fraction F dw 9 , which represents the burnup status of reactors and is analogous to the use of F dw 5 [36].In this data set, F 9 ranges from approximately 0.22 to 0.36, and F 5 , correspondingly, ranges from 0.66 to 0.49.The first group corresponds to F 9 between 0.22 and 0.24, due to low statistics, with the additional 12 groups each having a 0.01 interval in F 9 from 0.24 to 0.36.The effective fission fraction of the g th group, F g i , is calculated as N dw , where the information in each AD and week are added together if their F dw 9 's belong the g th group.The effective fission fractions averaged over all detectors and time ( F5 , F8 , F9 , and F1 ) are (0.564, 0.076, 0.304, and 0.056).
The energy differential IBD yield is measured for six reconstructed energy regions: 0.7-2, 2-3, 3-4, 4-5, 5-6, and 6-8 MeV and the energy differential yield, σ eg , for the e th energy region and the g th fission group is calculated as [5,36] where the integral is over the energy region, S dw (E rec ) is the measured energy spectrum of the d th AD in the w th week, the divisor gives the total number of fissions for the energy region, and the calculation of N dwe is the same N dw , except that the neutrino survival probability in Eq. 1 is calculated for the e th E rec region only.The sum over e is the total yield, σ g = e σ eg , of that group.The evolution of total and energy differential yield with F g 9 are plotted in Fig. 1.The uncertainties in σ g have statistical, background and the following systematic components.For the IBD detection efficiency, the AD-correlated uncertainty is improved from 1.7% to 0.75% [9], and the AD-uncorrelated uncertainty is 0.11% [37].The uncertainty of the number of target protons is 0.92% and is AD-correlated [5].The reactor power measurement uncertainty is 0.5% and is assigned to be reactor-uncorrelated and time-correlated [5].The uncertainty of the energy per fission is taken into account [44].The fission fraction uncertainty for the each isotope and reactor is 5%, but the uncertainties of the four isotopes are further constrained with the normalization condition and the correlation matrix [5] and are assigned to be reactorand time-correlated.The spent nuclear fuel uncertainty is improved from 100% to 30% [40].The nonequilibrium effect uncertainty is 30% [5].The θ 13 -induced oscillation uncertainty is also included [40].The uncertainty of the energy differential yield of σ eg further includes all the energy spectrum uncertainties from the background shape and detector response [37], in which the uncertainties in the absolute energy scale is reduced to be less than 0.5% for E rec larger than 2 MeV.
The predicted total and energy differential yields of the i th isotope, (σ 5 , σ 9 , σ 1 , and σ 8 ) and (σ e 5 , σ e 9 , σ e 1 , and σ e 8 ) are obtained by convolving the product of model prediction and IBD cross section [5] with the detector response matrix.The total yield predictions is defined as where σ i are the yields per isotope.Likewise, using the energy differential predictions, σ e i , we define the predicted energy differential yields The evolution plots of σ Pred,g and σ Pred,eg with F g 9 are shown in Fig. 1.The differences between the measured and predicted total and energy differential yields are also plotted as a function of F g 9 in Fig. 1.The uncertainties of σ Pred,g and σ Pred,eg are from all sources involved in the effective fission fraction calculation as described in Eq. 1, 3 and 4. Model uncertainties are poorly defined and not included unless explicitly stated otherwise.
The total yield evolution is compared to the predictions with two characteristic variables, average yield σ and normalized evolution slope (dσ/dF 9 )/σ.The average yield of σ and slope of dσ/dF 9 are two direct observables in Fig. 1.The evolution of the predicted yield can be described as a linear function of F 9 for the observed range of F 9 .In addition, if the prediction in Eq. 3 is off by a normalization factor η, for example, induced by large-mass sterile neutrinos [8,45,46] or by a global uncertainty, e.g. from the detection efficiency, the prediction would be The comparison in the normalized evolution slope (dσ/dF 9 )/σ is free of any normalization issue.
The total yield measurements in the 13 fission groups are fitted to the following linear function, with the χ 2 function, to extract σ and (dσ/dF 9 )/σ, where V is a 13×13 covariance matrix determined by randomly sampling all the related uncertainty sources described above.The best-fit results are σ = (5.89± 0.07) × 10 −43 cm 2 /fission and [(dσ/dF 9 )/σ] = −0.300± 0.024 with the χ 2 over the number of degrees of freedom, χ 2 /NDF, of 9.6/11.The dominant uncertainty of σ is from the IBD detection efficiency and number of target protons.The dominant uncertainty of (dσ/dF 9 )/σ is from statistics.The uncertainties from the effective fission fraction calculation are not significant for them.The best-fit line is shown in Fig. 1, and the results and 68% confidence level contour are shown in Fig. 2. For predictions, σ and (dσ/dF 9 )/σ can be directly calculated for a set of known fission fractions at Daya Bay.A joint distribution of σ and (dσ/dF 9 )/σ is obtained by randomly sampling the effective fission fractions according to their covariance matrix.The mean values and uncertainties of σPred and [(dσ/dF 9 )/σ] Pred are obtained with the distribution.The results for the HM are σHM = (6.18±0.04)× 10 −43 cm 2 /fission and [(dσ/dF 9 )/σ] HM = −0.387± 0.016 ((6.18 ± 0.16) × 10 −43 cm 2 /fission and −0.387 ± 0.018 if including the model uncertainties [6,7]).The HM prediction in σ and (dσ/dF 9 )/σ are rejected at 3.6 and 3.0 standard deviations.For SM2018, the results are consistent with the Daya Bay measurements.These results are shown in Fig. 2 and the best-determined lines are plotted in Fig. 1.
The energy differential yield evolution is compared to models with the average yields and normalized evolution slopes in six reconstructed energy regions.The data are simultaneously fitted to six linear functions, σ Lin,eg = σe {1 + [(dσ/dF 9 )/σ] e (F g 9 − F9 )}, with the χ 2 function, to extract six pairs of parameters of σe and [(dσ/dF 9 )/σ] e , where U is a 78 × 78 covariance matrix with a combined row (column) index of eg (e ′ g ′ ) for the e th (e ′th ) reconstructed energy region and g th (g ′th ) fission fraction group.U is also determined by a random sampling method of all the related uncertainty sources described earlier.The best-fit χ 2 /NDF is 76/66.The fit also gives the 12 × 12 covariance matrix of σe and [(dσ/dF 9 )/σ] e , which includes the 6 × 6 covariance matrix, A, for the six σe and the 6 × 6 covariance matrix, B, for [(dσ/dF 9 )/σ] e .The six σe results are strongly correlated because their dominant uncertainties are from the IBD detection efficiency and number of target protons, and the matrix A deviates strongly from a diagonal matrix.The six [(dσ/dF 9 )/σ] e results are all limited by data statistics and largely uncorrelated, and B is close to diagonal.The correlation between σe and [(dσ/dF 9 )/σ] e is insignificant.
For the predictions, a joint 12-dimension distribution of σe and [(dσ/dF 9 )/σ] e is obtained by randomly sampling the effective fission fractions as for the study of the predicted total yield and its normalized evolution rate.The mean values of σPred,e and [(dσ/dF 9 )/σ] Pred,e are obtained with the distribution as well as the covariance matrix for σPred,e , A Pred , and the covariance matrix for [(dσ/dF 9 )/σ] Pred,e , B Pred .The difference of σPred,e with the measurement and [(dσ/dF 9 )/σ] Pred,e results are plotted in Fig. 3.The uncertainty associated with prediction is much smaller than that from measurement.
The average IBD yields of six energy regions, σe are compared to the HM and SM2018 predictions σPred,e .Their difference is quantified as a χ 2 calculated with the difference of σe -σ Pred,e and their covariance matrix of A+A Pred .The resulting χ 2 /NDF and the corresponding extent of discrepancy in standard deviations are shown in Tab.I.The models do not agree with Daya Bay, and because of the deficit around 3 MeV and/or the bump around 5 MeV found in the measurement (Fig. 1) and the strong correlation among the measurements in different energy regions, their χ 2 /NDF's are rather large, and they correspond to 25 and 27 standard deviations for the HM and SM2018 models, respectively.The latter, due to the larger discrepancy in the 4-6 MeV region with the measurement, has a slightly worse χ 2 /NDF than HM.
The normalized evolution slopes of the six energy regions, [(dσ/dF 9 )/σ] e , are compared to HM and SM2018.Their difference is quantified with a χ 2 calculated with the difference of [(dσ/dF 9 )/σ] e -[(dσ/dF 9 )/σ] Pred,e and their covariance matrix of B+B Pred .The resulting χ 2 /NDF is shown in Tab.I.While the HM and SM2018 models poorly predict the spectral shape, their predicted relative changes with the fuel composition have much better agreement with the measurement.
To understand the difference between the Daya Bay differential IBD yield evolution and the predictions, three types of modified models with new free parameters are introduced on top of the HM and SM2018 predictions.
The first modification to each model is to alter only the 235 U energy differential yield prediction in each reconstructed energy region by the fraction f e 5 together with the global normalization factor η, as in Eq. 5, Depending on what the base model is, the modified models are further labelled as HM+ 235 U and SM2018+ 235 U.This is motivated by the fact that the majority of the the neutrino flux is due to 235 U.
In the second modification to each model, the prediction is where only the 239 Pu energy differential yield predictions in each reconstructed energy region is allowed to change by the fraction f e 9 together with the global normalization factor η. The modified models are labelled as HM+ 239 Pu and SM2018+ 239 Pu next.This is motivated given that 239 Pu is the second largest contributor to the neutrino flux.
The third modification to each model is to equally scale the predicted spectra of four isotopes in each reconstructed energy region by the fraction f e E , σ model,eg The motivation is that particular studies [30,47] have suggested that all four isotopes may have a common problem in predictions.They are labelled as HM+Equ and SM2018+Equ.
We fit the measured energy differential yields in the 6 energy regions and 13 fission fraction groups to the modified models with the following χ 2 where six f e 's and/or η are free parameters and Q is a 78 × 78 covariance matrix including all uncertainties for the measurement and predictions determined as V of Eq. 7 or U of Eq. 9. When using Eq. 10 or Eq.11, fits are also performed with η fixed to 1.The best-fit χ 2 /NDF, the corresponding extent of discrepancy in standard deviations, and best-fit η, when applicable, are shown in Tab.II.The best-fit f e 5 and f e 9 of Eq. 10 and Eq.11 with η fixed to 1, and f e E in Eq. 12 are shown in Fig. 4. The difference of the deduced σmodel,e with measurement and the deduced [(dσ/dF 9 )/σ] model,e are also shown in the figure, and where the first and third model modifications are preferred with respect to the second model.
Even when the 239 Pu energy spectra are modified, both the HM and SM2018 model predictions remain incompatible with the data at well over three standard deviations as shown in Tab.II.For both models, as seen in Fig. 4, the required changes of the 239 Pu spectrum in some regions are higher than 40%, which is far beyond the range of uncertainties caused by the various postulated mechanisms [25][26][27][28][29][30][31][32][33] and is unreasonable.This observation can be phenomenologically traced back to the features of Fig. 1.For example, the σe -σ HM,e in the 2-4 MeV region shows a positive slope and is not proportional to F 9 , which contradicts the assumption of pure 239 Pu-caused anomaly [38,39].
The attempts to adjust the predicted spectrum of 235 U or all spectra in equal measure all lead to good agreement with the data using this metric.As shown in Tab.II, their best-fit η results for 235 U-adjusting models are all consistent with 1.In summary, the total and energy differential IBD yield evolution as a function of fuel composition are measured and compared to the predictions of two signature models: the HM model based on the conversion method and the SM2018 model based on the summation method.While the measurement of the total IBD yield evolution is found to be incompatible with the HM model prediction, it is consistent with the SM2018 prediction.On the other hand, the predictions of spectrum evolution for both HM and SM2018 model show large discrepancies from the data.We exclude at high significance the hypothesis that the 239 Pu energy spectrum in HM or SM2018 models is responsible for the entire difference with the data, regardless of how the normalization of the Daya Bay data is treated.In contrast, good consistency with the data can be achieved either by altering the 235 U spectrum or all four isotopes' spectra in equal measure in the SM2018 model.For the HM model, the 235 U spectrum adjustment works slightly better than adjusting all spectra, as indicated by the total yield evolution measurement.Future enhancements to the models could prioritize 235 U-specific causes or factors common to the four isotopes.
Daya Bay is supported in part by the Ministry of Science and Technology of China, the U.S. Department of Energy, the Chinese Academy of Sciences, the National Natural Science Foundation of China, the Guangdong provincial government, the Shenzhen municipal government, the China General Nuclear Power Group, Key Laboratory of Particle and Radiation Imaging (Tsinghua University), the Ministry of Education, Key Laboratory of Particle

3 FIG. 1 .
FIG.1.The panels a.1 and b.1 show the total IBD yields in [0.7,8] MeV and energy differential yield in six reconstructed energy regions as a function of the effective fission fraction of 239 Pu, F9, respectively.The best-fit and best-determined lines for the measurements and predictions of the evolution of the total yield are shown in a.1, respectively.The difference between the measurement and the HM and SM2018 predictions for the total yield (a.2 and a.3) and energy differential yields (b.2 and b.3) are also shown.The error bars represent the statistical uncertainties.The units of all panels are 10 −43 cm 2 /fission.

FIG. 3 .
FIG.3.The upper panel shows the difference between the measured energy differential yields and predictions for six reconstructed energy bins, where the error bars are from the measurement.The lower panel shows the normalized evolution slopes for the measurement and predictions, where the uncertainties of measurement are shown.

FIG. 4 .
FIG. 4. The best-fit f e , i.e. f e 5 , f e 9 or f e E , values of the modified models of HM+ 235 U, SM2018+ 235 U (Eq. 10 with η fixed to 1), HM+ 239 Pu, SM2018+ 239 Pu (Eq.11 with η fixed to 1), HM+Equ, and SM2018+Equ (Eq.12) are shown in the upper panels, where the error bars are fit results.The deduced σmodel,e predictions with the corresponding f e values for each model are shown as the difference with the measurement in the middle panels and the error bars shown are from the measurement.The measured [(dσ/dF9)/σ] e and deduced [(dσ/dF9)/σ] model,e are shown in the lower panels and only the error bars of measurement are shown.

TABLE II .
For the six modified models in Eq. 10, 11, and 12 (the first column), the best-fit χ 2 /NDF when fitting to data and the corresponding number of standard deviations are shown in the second column and the determined normalization factor η in the third column.Trials are also done with η fixed to 1 for Eq. 10 and Eq.11.The deduced σmodel,e and [(dσ/dF 9 )/σ] model,e are consistent with the measurements as shown in Fig.4.HM+ 235 U works slightly better than HM+Equ model, as their best-fit χ 2 /NDF shown in Tab.II.But with the current precision of the Daya Bay data set, it is difficult to distinguish whether 235 U, by itself, or a mix of fission isotopes, are responsible for the flux and spectrum anomalies.