Antineutrino energy spectrum unfolding based on the Daya Bay measurement and its applications

Nuclear reactors are a powerful source of electron antineutrinos ($ \overline{\nu}_{e} $) and have played a significant role in neutrino physics, including the discovery of neutrinos [1], the measurement of the neutrino mixing angle $ \theta_{12} $ and the neutrino mass-squared splitting $ \Delta m^2_{21} $ [2], and the observation of the neutrino oscillation driven by $ \theta_{13} $ [3-5]. Many very short baseline experiments are producing exciting results in the search for sterile neutrinos and in precise measurements of the reactor antineutrino energy spectrum [6-10]. Looking forward, resolution of the neutrino mass ordering is the design goal of the JUNO reactor neutrino experiment [11] at a baseline of 53 km.

Commercial pressurized water reactors produce a large number of $ \overline{\nu}_{e} $ 's emitted from the beta decay chains of the fission products from four main isotopes ²³⁵U, ²³⁸U, ²³⁹Pu, and ²⁴¹Pu, while other isotopes contribute less than 0.3%. In general, about $ 2\times 10^{20} $ antineutrinos per second are released per GW thermal power. The understanding of the antineutrino spectra is a key issue for reactor antineutrino experiments. At present there are two methods to obtain the predicted antineutrino spectra from these four fission isotopes.

The first method is the conversion method. It is based on the measured beta spectra from thermal-neutron induced fission of ²³⁵U, ²³⁹Pu, and ²⁴¹Pu performed at the ILL High Flux Reactor [12-14]. The electron spectra of fission isotopes were fitted with a set of virtual beta decay branches based on allowed beta decay transitions. In 2011, Mueller et al. and Huber re-evaluated the flux and spectra including an improved beta spectrum calculation [15-17] (the Huber-Mueller model). The uncertainties were estimated from detailed studies of corrections to the allowed $ \beta $-spectrum shape, the inversion errors based on synthetic data sets, and the reliability of nuclear structure data. The re-evaluated flux was found to be about 5% higher than past measurements [18]. The discrepancy is commonly referred to as the "Reactor Antineutrino Anomaly" (RAA). Precise measurements of the reactor antineutrino spectral shape indicated another discrepancy (spectral distortion) compared to the Huber-Mueller model prediction [19-24]. A number of different hypotheses have been proposed to explain this discrepancy, such as improper treatment of the shape corrections for forbidden transitions [25].

The second method is the summation method. The total $ \overline{\nu}_{e} $ energy spectra of all known fission decay channels are calculated based on the fission yields of the fission products, Q values and decay branching fractions in the nuclear data libraries. This method generally has unknown uncertainties because the correlation among different sets of nuclear data uncertainties have not been properly cataloged. In addition, this method suffers from unknown uncertainties due to missing or incomplete information in the nuclear data libraries. It was recently pointed out that the bias of the pandemonium effect (the inability to accurately measure the complicated beta spectra for nuclei when the energy available for beta decay is large [26]) on nuclear structure measurements also impacts antineutrino summation predictions [27]. Including pandemonium-free data of the major contributors can decrease the antineutrino flux from all fission isotopes [28], thus improving agreement with the measured flux from the fuel evolution study at Daya Bay [29]. Nonetheless, a recent comprehensive analysis with updated reactor antineutrino flux models showed the rate anomaly and spectral anomaly still persist [30].

The above studies illustrate the difficulty of constructing a reliable model because of the complexity of beta decay theory and the imperfect information in nuclear data libraries. To minimize the uncertainty due to the prediction of the antineutrino spectrum, a relative measurement comparing antineutrino events at near and far site detectors is utilized in the experiments aiming at the measurement of $ \theta_{13} $ [3-5]. In addition, the relative measurement technique is also utilized in experiments with segmented detectors for short-baseline sterile neutrino searches [7, 8]. The antineutrino flux and spectrum with little influence from neutrino oscillation are well measured in the near-site detectors. The uncertainty of the total measured antineutrino flux and spectrum is better than that of model predictions for commercial reactors [31, 32]. These absolute spectrum measurements can be used as an additional resource for validating standard nuclear databases [33]. Moreover, these measurements can also provide a reference spectrum for other reactor antineutrino experiments, especially with single unsegmented detectors. These measurements have been utilized in the studies of reactor antineutrino experiments via the inverse beta decay (IBD) reaction [11, 21]. These precise measurements can also be used as an input for reactor antineutrino experiments utilizing coherent elastic neutrino-nucleus scattering (${\rm{CE}} \nu {\rm{NS}}$) [34, 35].

Reactor antineutrinos are measured via the IBD reaction at Daya Bay. The positron carries most of the antineutrino energy and forms a prompt signal in the detector. In Ref. [32], the measured prompt energy spectrum based on 1958 days of data acquisition was provided in bins of 0.25 MeV. Moreover, the prompt energy spectra of ²³⁵U and ²³⁹Pu were extracted for the first time based on the commercial reactor data. To provide a data-driven prediction for other experiments with different fission fractions than Daya Bay, the correlation of the total measured prompt energy spectrum and the extracted prompt energy spectra for specific fission isotopes is obtained(the term "isotopic energy spectrum" will be used to refer to an energy spectrum for a specific fission isotope in this manuscript). Furthermore, the corresponding antineutrino energy spectra are presented and utilized in this paper, by removing the contribution from detector response with the Wiener-SVD unfolding method [36]. Finally, a new method is utilized to predict the antineutrino spectrum for other experiments at reactors with arbitrary fission fractions.

This paper is organized as follows: Sec. II introduces the antineutrino detectors, the previous studies, the event selection, and the detector response at Daya Bay. Sec. III provides the correlation of the total measurement and the extracted isotopic spectra. Sec. IV presents the generic $ \overline{\nu}_{e} $ energy spectra of the IBD reaction. Application of generic $ \overline{\nu}_{e} $ energy spectra is explained in detail in Sec. V. A short summary is given in Sec. VI.

The Daya Bay experiment studies the flux and spectrum of $ \overline{\nu}_{e} $'s from six $ 2.9\; {\rm{GW}} $ thermal power commercial pressurized-water reactor cores at the Daya Bay nuclear power complex. Eight identically designed antineutrino detectors (ADs) are distributed at two near-site experiment halls (EH1 and EH2) and one far-site experiment hall (EH3). The near site halls are used to monitor the $ \overline{\nu}_{e} $ flux and spectrum from the reactor cores (with flux-weighted baselines of 560 and 600 m for EH1 and EH2, respectively), while the far-site experiment hall is used to measure the oscillated spectrum and flux deficit driven by $ \theta_{13} $ (with flux-weighted baselines of 1640 m). Each AD consists of three nested cylindrical vessels. The inner acrylic vessel (IAV) contains 20 tons of 0.1% gadolinium-loaded liquid scintillator (GdLS) and serves as the central $ \overline{\nu}_{e} $ target. The outer acrylic vessel surrounding the target creates a 42 cm thick pure liquid scintillator (LS) region to improve the collection of gammas escaping from the GdLS region to reduce the energy leakage. Scintillation light is detected by 192 8-inch PMTs (Hamamatsu R5912), which are positioned on the outermost cylinder of stainless steel. A mineral oil buffer shields the GdLS and LS regions from gamma rays from natural radioactivity in the PMTs. The experimental setup is described in detail in Refs. [37, 38].

In a commercial reactor core, the chain reaction enabled by the neutrons produced through the fission of uranium and plutonium isotopes maintains the overall burn-up of the nuclear reactor fuel. The fissile isotopes ²³⁹Pu and ²⁴¹Pu are readily produced by neutron capture on ²³⁸U and ensuing reactions. The fraction of nuclear fissions attributed to a parent isotope, such as ²³⁹Pu, is called the fission fraction. The fission fractions of different isotopes evolve with burn-up. Figure 1 shows an example of the evolution of the fission fractions as a function of burn-up within a refueling cycle [31]. The fission fraction of ²³⁵U decreases with burn-up while the fission fractions of ²³⁹Pu and ²⁴¹Pu increase. The fission fraction of ²³⁸U is relatively small (7.6%) and stable over time. Over 99.7% of the $ \bar\nu_e $'s produced from thousands of beta decays of fission daughters are due to the four isotopes ²³⁵U, ²³⁹Pu, ²³⁸U and ²⁴¹Pu.

Figure 1.  (color online) Fission fractions of isotopes as a function of burn-up from a simulation of a reactor core at the Daya Bay nuclear power plant. Other isotopes contribute less than 0.3% in total.

The Daya Bay experiment has measured the total flux and spectrum of reactor antineutrinos [19, 20, 31, 39]. Based on the evolution with respect to the fission fractions, the $ \overline{\nu}_{e} $ flux and spectra of two primary fission isotopes, ²³⁵U and ²³⁹Pu, were extracted [29, 32]. These results have contributed significantly to the understanding of the RAA.

The first measurement of the antineutrino flux and the prompt energy spectrum with ~0.3 million IBD events were reported in Refs. [19, 20]. A ~5% flux deficit, and a ~10% spectral distortion in the 4-6 MeV prompt energy region were found, relative to the Huber-Mueller model predictions. The details of this analysis were published in Ref. [31] with an increased IBD sample of about 1.2 million events. Moreover, the antineutrino energy spectrum weighted by the IBD cross section was obtained. It has provided model-independent predictions for other reactor antineutrino experiments. To reduce the systematic uncertainties in the flux and spectrum measurements, the detection efficiency and the energy response were carefully determined [39, 40]. These efforts led to the improved flux measurement in Ref. [39], and the improved spectrum measurement in Ref. [32]. Two novel measurements were performed by introducing the information about reactor fuel burning. The antineutrino fluxes of ²³⁵U and ²³⁹Pu were extracted in Ref. [29] using 2.2 million IBD candidates. In Ref. [29], a 7.8% discrepancy between the measured and predicted ²³⁵U antineutrino yield was found, indicating that ²³⁵U is likely the primary contributor to the RAA. In addition, the measurement at Daya Bay disfavored the hypotheses of equal deficits for all fission isotopes and a ²³⁹Pu-only deficit as the reason for the discrepancy with the predicted flux at the $ 2.8\sigma $ and $3.2 \sigma $ confidence levels, respectively. The individual prompt energy spectra of ²³⁵U and ²³⁹Pu were extracted for the first time in commercial reactors [32] using 3.5 million IBD candidates. The work in this paper uses the same data set as Ref. [32].

This analysis follows the work in Ref. [32], which makes use of 3.5 million inverse beta decay candidates from the four near-site ADs. IBD candidates are selected following the same criteria as Ref. [41]. After applying these cuts, the estimated signal and background rates, as well as the efficiencies of the muon veto, $ \epsilon_{\mu} $, and multiplicity selection, $ \epsilon_{\rm m} $, were determined [37]. Backgrounds remaining in the IBD samples from accidental coincidences, fast neutrons, cosmogenic ⁸He/⁹Li production, AD-intrinsic alpha radioactivity, and AmC neutron calibration sources were estimated using a variety of techniques described in detail in a previous publication [42]. Additionally, spent nuclear fuel (SNF) present in the cooling pool adjacent to each reactor core contributes ~0.3% to the IBD rate [31]. The effect of the SNF was subtracted from the IBD spectrum. The relative uncertainty from SNF was reduced from 100% [41] to 30% [43] using the spent fuel inventory history provided by the nuclear power plant. The relative rate of background is less than 2% and contributes less than 0.15% to the uncertainty on the IBD rate and less than 0.02% to the uncertainty of IBD shape in the 1.5-6.0 MeV prompt energy region.

The reconstructed positron energy relative to the true interaction energy is nonlinear due to ionization quenching and the Cherenkov light emission in the scintillator, and the underestimate of the charge of the PMT signals in the readout system. Previously, the energy nonlinearity model was constructed with 1% uncertainty based on a semi-empirical analytic approach [41]. In December 2015, a full Flash-ADC (FADC) readout system was installed at EH1 AD1, which recorded the PMT waveforms simultaneously with the previous front end electronics (FEE) systems. Based on a deconvolution method, the integral charge was extracted with a minimum bias based on the PMT waveform recorded by the FADC [44]. The integral charge from the FADC was compared with the reconstructed charge from the FEE. The uncertainty on the electronics nonlinearity was reduced to 0.2% based on an event-by-event comparison of the total charge from these two readout systems. In addition, the uncertainty in the visible energy from $ \gamma $ rays was improved from 1% to 0.5% after a special calibration campaign in January 2017 that deployed ⁶⁰Co sources with different enclosures to quantify the optical shadowing effect. Finally, the statistics of cosmogenic ¹²B candidates detected in four near-site ADs increased to ~470000 which were used to further refine the nonlinearity model. Based on these improvements, the uncertainty on the energy nonlinearity was improved to <0.5% for prompt energy larger than 2 MeV. More details can be found in Refs. [37, 40].

Comparison of the mean reconstructed energy between antineutrino detectors was done with a variety of calibration references (neutron-capture on hydrogen and gadolinium, gammas from external ⁴⁰K and ²⁰⁸Tl decays, and $ \alpha $'s from ²¹²Po, ²¹⁵Po and ²¹⁹Po decays). The variations were less than 0.2% [41].

As mentioned above, the IAV contains the GdLS ($ \overline{\nu}_{e} $ target), but the IAV is non-scintillating material. Positrons and the annihilation $ \gamma $-rays from IBD reactions around or inside the acrylic may lose energy invisibly in the acrylic, and this phenomenon leads to energy leakage and a slight distortion of the prompt energy spectrum. The relative uncertainty on the measured prompt energy spectrum from this effect (called the "IAV effect" in this manuscript) was estimated to be 4% (0.1%) below (above) 1.25 MeV by using simulation and it is assumed to be correlated among detectors [31].

The energy resolution of the detectors was studied based on the measured reconstructed energy spectra of a variety of calibration sources deployed at the detector center, IBD and spallation neutrons, and alpha sources from radioactivity. The relative energy resolution of the ADs as a function of energy was modeled using the expression

where $ \sigma_E $ is the standard deviation of the reconstructed energy $ E_{\rm{rec}} $, and the parameters $ a $, $ b $, and $ c $ quantify the contribution to the resolution from detector energy nonuniformity, photoelectron statistics, and PMT dark noise, respectively. The best-fit parameters of the model are $ a $ = 0.016, $ b $ = 0.081 MeV^1/2 and $ c $ = 0.026 MeV. The relative degradation of the energy resolution is less than 6% after 6 years of operation (~8.5% to ~9.0% at $ E_{\rm{rec}} $ = 1 MeV), and the degradation has negligible affect on following analysis.

The response matrix mapping the $ \overline{\nu}_{e} $ energy to the reconstructed energy was constructed using a full-detector simulation based on Geant4 [45]. In this case, the response matrix included the IBD energy shift (the energy shift when $ \overline{\nu}_{e} $ energy is transferred to a positron and reconstructed as prompt energy), the IAV effect, and energy resolution. The reconstructed energy of IBD events is corrected for energy nonlinearity, which means the reconstructed energy of IBD events is the energy deposited in the LS. The response matrix was provided in the supplemental materials of Ref. [32]. For convenience, it is provided with the supplemental materials of this paper as well.

The total measurement of the prompt energy spectrum in this paper is based on the data from the four near-site ADs. For each AD, the detected $ \overline{\nu}_{e} $ 's are from 6 reactor cores with different baselines. Antineutrinos from each reactor come from four main isotopes with different relative contributions, hence the total prompt energy spectrum is the sum of the four isotopic energy spectra of prompt signal. To define the fractional contribution from different isotopes to the total number of IBD events in near-site ADs, the effective fission fraction of the total prompt energy spectrum is defined as:

Here $ f_r^{\rm{iso}} $ is the fission fraction of $ r $-th core for each isotope; $ L_{dr} $ is the distance between the $ d $-th detector and the $ r $-th reactor core; $ N_r^f $ is the predicted number of fissions for each isotope from the $ r $-th reactor core, and it is calculated based on:

where $ t $ represents the detector data acquisition (DAQ) live time, $ W_r $ represents the thermal power of $ r $-th core, $ E_{\rm{iso}} $ represents the mean energy release per fission for each isotope. The average effective fission fractions for ²³⁵U, ²³⁸U, ²³⁹Pu and ²⁴¹Pu in near-site ADs for the analyzed data set are $ f^{235}_{\rm{eff}}:f^{238}_{\rm{eff}}:f^{239}_{\rm{eff}}:f^{241}_{\rm{eff}} = 0.564: 0.076: $$ 0.304: 0.056 $. The detailed description of the measurement as well as the comparison with the Huber-Mueller model can be found in Ref. [32].

As the observed prompt energy spectrum evolves as a function of fission fractions, the isotopic prompt energy spectra from ²³⁵U and ²³⁹Pu have been extracted at Daya Bay [32]. Daya Bay data are less sensitive to individual spectra of ²³⁸U and ²⁴¹Pu due to their smaller fission fractions than ²³⁵U and ²³⁹Pu. Constraints on individual spectra of ²³⁸U and ²⁴¹Pu are needed to obtain the dominant isotopic prompt energy spectra from ²³⁵U and ²³⁹Pu. The constraints on the prompt energy spectra of ²³⁸U and ²⁴¹Pu are given in Ref. [32]. Since the fission fraction of ²⁴¹Pu evolves approximately proportional to ²³⁹Pu (the ratio is fitted to be 0.183), the spectra of ²³⁹Pu ($ s_{239} $) and ²⁴¹Pu ($ s_{241} $) can be treated as one component, defined as $ s_{\rm{combo}} = s_{239} + R^{\rm{Pu}} \times s_{241} $. Here $ R^{\rm{Pu}} $ should be a constant to make $ s_{\rm{combo}} $ an invariant and $ R^{\rm{Pu}} $ is chosen to be 0.183 to have the least dependence on the input of the ²⁴¹Pu spectrum, thus leading to the best precision of Pu combo spectrum. The residual contribution of ²⁴¹Pu spectrum on the data was corrected when the fission fraction ratios of ²⁴¹Pu to ²³⁹Pu deviate from 0.183. The extracted isotopic (²³⁵U and Pu combo) prompt energy spectra of the IBD reaction were published in Ref. [32], which provided uncertainties for both the flux and spectral shape.

The isotopic energy spectra have larger uncertainties than the total observed prompt energy spectrum since the individual components extracted from the total one have dependence on extra information, which brings in additional uncertainties. The total prompt energy spectrum can provide a data-driven prediction for other reactor $ \overline{\nu}_{e} $ experiments with similar fission fractions. The isotopic prompt energy spectra can be used to correct effects from the differences in fission fractions between experiments. Since the total and isotopic prompt energy spectra are derived from the same data set, their uncertainties are correlated, which has to be taken into account when using them together. The total prompt energy spectrum, $ s_{\rm{total}} $, is the sum of isotopic prompt energy spectra times the corresponding effective fission fractions:

The spectrum of Pu combo is used in this analysis because it has a smaller uncertainty than the total uncertainty of the extracted prompt energy spectrum of ²³⁹Pu and the correlation between the fission fractions of ²³⁹Pu and ²⁴¹Pu in most commercial low-enriched uranium (LEU) nuclear power reactors are similar. The ²³⁹Pu spectrum can be extracted by subtracting prediction of the prompt energy spectrum based on the ²⁴¹Pu model.

Based on the covariance matrix among the isotopic prompt energy spectra from the extraction algorithm [32], the total prompt energy spectrum and its correlation with the other two isotopic prompt energy spectra ($ s_{235} $ and $ s_{\rm{combo}} $) are obtained with Eq. (4) through standard error propagation. The total and isotopic prompt energy spectra, as well as the covariance matrix are provided in the supplemental materials. Dominant components of the energy-dependent uncertainties for total and isotopic (²³⁵U and Pu combo) prompt energy spectra are shown in Fig. 2. Detector response uncertainty contains the contribution of the uncertainties from detection efficiency, the energy nonlinearity model, the energy scale difference between ADs, and the IAV effect. Model (²³⁸U, ²⁴¹Pu) uncertainty represents uncertainties of the input ²³⁸U and ²⁴¹Pu model in the analysis to extract the isotopic (²³⁵U and Pu combo) prompt energy spectra. The uncertainties of the isotopic prompt energy spectra are dominated by statistical uncertainty and the uncertainty of the models (for ²³⁸U and ²⁴¹Pu). Other uncertainties (e.g. background uncertainty) have small contribution to the total uncertainty and they are not shown in the plot. The calculated total measured prompt energy spectrum and its uncertainty in this analysis are consistent with the previous results [32]. The uncertainty of the total measurement is dominated by the detector response uncertainty, especially the uncertainty of the detection efficiency (1.19%). The primary not-fully-energy-correlated uncertainties for the total measurement are dominated by the uncertainty of the energy nonlinearity model. There is no contribution of the energy-dependent uncertainties from the input ²³⁸U and ²⁴¹Pu models to the total measurement since it is a direct measurement.

Figure 2.  (color online) Dominant components of the energy-dependent uncertainties for isotopic (²³⁵U and Pu combo) and total prompt energy spectra. Here relative uncertainty represents the square root of the diagonal elements of the covariance matrix. "Detector" uncertainty (detector response uncertainty) contains the contribution of the uncertainties from detection efficiency ("Det. Eff."), the energy nonlinearity model ("Nonlinearity"), energy scale difference between ADs ("Rel. ES"), and the IAV effect ("IAV"). "Model (²³⁸U, ²⁴¹Pu)" uncertainty represents uncertainties of the models for ²³⁸U and ²⁴¹Pu in the analysis to extract the isotopic (²³⁵U and Pu combo) prompt energy spectra. The inset shows the components of the detector response uncertainty for the total prompt energy spectrum, with statistical uncertainty shown for comparison. Other uncertainties (e.g. background uncertainty) have a small contribution to the total uncertainty and are not shown.

The correlation matrix of the isotopic and total prompt energy spectra is shown in Fig. 3. The total uncertainty including both the rate and spectral shape of the total prompt energy spectrum is ~1.3% in 2 to 5 MeV energy region, with large bin-to-bin correlation due to the dominant detection efficiency uncertainty. The relative uncertainty in the spectral shape is smaller (<0.5% in 2-5 MeV energy region) than the overall uncertainty for the total prompt energy spectrum [32]. The correlation between the same bin of the prompt energy spectra for ²³⁵U and Pu combo is mostly less than 0 since the isotopic prompt energy spectra constitute the total measurement and their statistical fluctuations are anticorrelated.

Figure 3.  (color online) Correlation matrix of the isotopic (²³⁵U and Pu combo) and total prompt energy spectra. The first, middle, and last 26 bins of the correlation matrix represent the prompt energy spectra of ²³⁵U, Pu combo, and total measurement in Fig. 2 respectively.

The total prompt energy spectrum and the individually extracted isotopic prompt energy spectra of the IBD reaction can be compared with antineutrino models after taking into account the detector response, and are used to obtain the generic antineutrino energy spectra of the IBD reaction in Sec. IV.

The systematic uncertainty dominates the total uncertainty for the total prompt energy spectrum with the bin width of 0.25 MeV as shown in Fig. 2. Since the statistical uncertainties are within sub-percent level at all energies, the bin width is further reduced to better match the energy resolution of the Daya Bay detectors. This can also provide more spectral information, thus more studies can be carried out. The Appendix contains a discussion of the measurement and applications of the finely-binned prompt energy spectrum.

To provide a model-independent reactor $ \overline{\nu}_{e} $ energy spectrum prediction, unfolding techniques are used to obtain generic antineutrino energy spectra of the IBD reaction (or generic $ \overline{\nu}_{e} $ energy spectrum weighted by the IBD cross section). Unfolding transforms the prompt energy spectrum to the $ \overline{\nu}_{e} $ energy spectrum for direct comparison with spectra from other experiments as well as with theoretical predictions. Daya Bay previously used the SVD regularization method to produce the first generic $ \overline{\nu}_{e} $ energy spectrum [31]. In this Section, unfolding using a new method, the Wiener-SVD method [36], is described and used to produce total and isotopic (²³⁵U and Pu combo) $ \overline{\nu}_{e} $ energy spectra. The spectra are compared to results obtained with the SVD [46] and Bayesian iteration [47] methods.

Because of the finite energy resolution, the response matrix is found to be ill-conditioned with close-to-zero singular values. If the inverse of the detector response matrix is naively used to obtain the $ \overline{\nu}_{e} $ energy spectrum, statistical fluctuations and systematic variations will be significantly amplified. This phenomenon would lead to meaningless results. Thus the unfolding techniques are utilized to suppress the fluctuations and solve this problem. In this work, the unfolded results of Bayesian iteration method and SVD method are based on an unfolding package RooUnfold [48], while the results of Wiener-SVD method are based on the algorithm shared in the GitHub [36]. The Bayesian iteration method needs an initial guess of the $ \overline{\nu}_{e} $ energy spectrum as a starting value and the spectrum is updated iteratively based on the Bayes' theorem, which takes into account the response matrix and the observed prompt energy spectrum. The iteration is stopped when the change of the unfolded spectrum is small and before the unphysical fluctuations appear [47]. The SVD method is based on matrix inversion calculation, reducing or suppressing the small singular values in the detector response matrix with a regularization parameter. This method is improved with $ a\; priori $ information about the solution to overcome the instability of unfolding [46]. The Wiener-SVD approach achieves the unfolding by maximizing the signal to noise ratio in the effective frequency domain (a Wiener filter constructed based on an input $ \overline{\nu}_{e} $ model and the observed spectrum) thus avoiding the need for a regularization parameter used by other methods. It also has smaller mean square error (MSE, which averages the total bias and variance) than traditional SVD regularization method [36].

As described above, a model of $ \overline{\nu}_{e} $ energy spectrum is needed as an input in each of these unfolding procedures. As shown in Refs. [9, 10, 32, 49], both ²³⁵U and ²³⁹Pu are likely to be responsible for the spectral distortion, rather than only the ²³⁵U contribution. In this case, the Huber-Mueller model with spectral correction is used as an input model. The correction of the spectral distortion is the ratio of the total measurement of the prompt energy spectrum over the prediction based on the Huber-Mueller model at Daya Bay. The ratio defined above is based on prompt energy, and the prompt energy ($ E_{\rm{p}} $) can be shifted to $ \overline{\nu}_{e} $ energy ($ E_{\bar{\nu}_e} $) using the following formula: $ E_{\bar{\nu}_e}\approx E_{\rm{p}}+0.78\; {\rm{MeV}} $. The ratio obtained above is applied to all Huber-Mueller model predictions (total, ²³⁵U, and Pu combo $ \overline{\nu}_{e} $ energy spectra) to construct an input model for the unfolding method. However, the input model is not identical to the true $ \overline{\nu}_{e} $ spectrum, and this effect will induce additional bias in the unfolded results. Unfolding methods suppress the fluctuations as well as the fine structures in the original $ \overline{\nu}_{e} $ energy spectra. With different variations added on the input model, changes in the unfolded results are observed. A study was done to consider the effect of fine structures based on the SM2018 model [28] which found a ~2% variation compared with a smooth spectrum (smearing with the energy resolution of Daya Bay detectors). This variation is consistent with the sawtooth distortions found in Ref. [50]. Based on this study, the possible variation due to the difference between the input model and the true spectrum is set to 3% (bin-to-bin uncorrelated) conservatively.

The relative performance of these methods also depends on the specific convergence criteria applied. The Bayesian iteration method needs the number of iterations as an input, while the SVD method needs input for the regularization parameter. To optimize the input parameters for these methods and compare the MSE between these methods, dedicated tests with toy Monte Carlo simulations are done. In the tests, $ \overline{\nu}_{e} $ energy spectra based on the Huber-Mueller model ($ {{S}}^{\bar{\nu}_e} $) are used to generate 10000 samples of total and isotopic prompt energy spectra with the response matrix and fluctuations based on the covariance matrix in Sec. III. The total and isotopic unfolded $ \overline{\nu}_{e} $ energy spectra ($ {{S}}^{\rm{unfold}} $) are obtained with each unfolding method for each sample. Based on a dedicated study on the response matrix, the binning method for the unfolded $ \overline{\nu}_{e} $ energy spectra is optimized with totally 25 energy bins, while originally the prompt energy spectra have 26 energy bins in total. For each toy Monte Carlo test, total and isotopic $ \overline{\nu}_{e} $ energy spectra are unfolded individually, but they are combined to one spectrum in order to evaluate their correlation more conveniently. Then $ {{S}}^{\rm{unfold}} $ are compared to the true $ \overline{\nu}_{e} $ energy spectra ($ {{S}}^{\bar{\nu}_e} $) to evaluate the bias and uncertainty of the unfolded result. Here $ {{S}}^{\rm{unfold}} $ and $ {{S}}^{\bar{\nu}_e} $ are defined to contain the total, ²³⁵U and Pu combo spectra, with 75 (= 25$ \times $3) energy bins chaining all these three spectra together (each with 25 energy bins). In each unfolding procedure, a 3% bin-to-bin uncorrelated fluctuation is added on the true $ \overline{\nu}_{e} $ energy spectra to construct different input $ \overline{\nu}_{e} $ spectra, in order to account for the possible variation from input models. For input $ \overline{\nu}_{e} $ isotopic spectra, an additional 10% systematic fluctuation (bin-to-bin correlated) is added to account for the possible size of the spectral distortion. An element of the covariance matrix is calculated by

where $ t $ is the sample number and runs from 1 to $ N $ = 10000, $ i $ represents the $ i $-th energy bin in the spectra. $ {{Cov}} $ contains the effect of the bias from the input model and the uncertainty from the measurement. The best number of iterations for the Bayesian iteration method is determined to be 1 while the regularization parameter for the SVD method is determined to be 21, with the least MSE for toy Monte Carlo tests. Here the MSE is defined as the summation of the square root of the diagonal terms in $ {{Cov}} $. $ {{Cov}} $ contains the information of the correlation among total and isotopic unfolded $ \overline{\nu}_{e} $ energy spectra.

Figure 4 shows the unfolded $ \overline{\nu}_{e} $ energy spectra as well as their relative uncertainties. The total uncertainties of the unfolded $ \overline{\nu}_{e} $ energy spectra from the Wiener-SVD method in the 3 to 6 MeV energy region are 2% for the total $ \overline{\nu}_{e} $ energy spectrum and 3.5% (5%) for the ²³⁵U (Pu combo) $ \overline{\nu}_{e} $ energy spectrum. The Bayesian iteration and SVD methods yield consistent results. The result from the Wiener-SVD method has the least MSE in the 3 to 6 MeV energy region, especially for the total $ \overline{\nu}_{e} $ energy spectrum. This is presumably due to the optimized signal to noise ratio in the effective frequency domain for the Wiener filter, as found in Ref. [36]. Since the Wiener filter has small suppression in the low effective frequency domain, the result from the Wiener-SVD method has larger fluctuation and bigger uncertainties than traditional Bayesian iteration method and SVD method in the low and high energy regions, where the original statistical fluctuations are large.

Figure 4.  (color online) (Top panels) Unfolded isotopic (²³⁵U and Pu combo) and total $\overline{\nu}_{e}$ energy spectra weighted by IBD cross section. The error bars on the data points represent the square root of diagonal elements of the covariance matrix, which is generated by Monte Carlo methods. (Middle panels) Relative uncertainty of the unfolded isotopic and total $\overline{\nu}_{e}$ energy spectra. (Bottom panels) Enlarged plot of the relative uncertainty of the unfolded isotopic and total $\overline{\nu}_{e}$ energy spectra. The result from the Wiener-SVD method has the smallest MSE in the 3 to 6 MeV energy region, especially for the total $\overline{\nu}_{e}$ energy spectrum.

To gain insight on the contribution from each component of the energy-dependent uncertainties to the unfolded $ \overline{\nu}_{e} $ energy spectra, dedicated tests with toy Monte Carlo simulations were done. Toy Monte Carlo samples of total and isotopic prompt energy spectra were generated with the response matrix and fluctuations based on the error budget in Fig. 2. In each unfolding procedure, the input $ \overline{\nu}_{e} $ energy spectra are exactly the same as the true $ \overline{\nu}_{e} $ energy spectrum for all toy Monte Carlo tests. The total and isotopic unfolded $ \overline{\nu}_{e} $ energy spectra ($ {{S}}^{\rm{unfold}} $) are obtained with each unfolding method and are used to calculate elements of the covariance matrix based on the following formula:

Here $ \bar{{{S}}}^{\rm{unfold}} = \dfrac{1}{N}\displaystyle\sum\nolimits_t^N{{{S}}^{\rm{unfold}}_{t}} $ is the average of the unfolded $ \overline{\nu}_{e} $ energy spectra $ {{S}}^{\rm{unfold}} $. With toy Monte Carlo samples incorporating fluctuations based on the covariance matrixes of different components in Fig. 2, contributions from the uncertainties due to statistics, detector, and model (²³⁸U, ²⁴¹Pu) to the total uncertainties of the unfolded results are obtained. The total covariance matrix without the bias from the input model ($ ({{Cov}}^{\prime})^{\rm{total}} $) is obtained based on toy Monte Carlo samples taking into account the total uncertainty from the measurement, using Eq. (6). Since the covariance matrix ($ {{Cov}}^{\rm{total}} $) from Eq. (5) takes into account the bias from the input model, the additional bias from the input model ($ {{Cov}}^{\rm{bias}} $) is calculated by the following equation:

Based on these tests, the dominant components of the energy-dependent uncertainties for total and isotopic unfolded $ \overline{\nu}_{e} $ energy spectra based on the Wiener-SVD method are shown in Fig. 5. Contributions from statistics and the models for ²³⁸U and ²⁴¹Pu dominate the total uncertainties for isotopic spectra. The input model used by the unfolding method and the detector response uncertainty both induce an uncertainty on the total unfolded $ \overline{\nu}_{e} $ energy spectrum at the 1% level and dominate the total uncertainty in the spectrum.

Figure 5.  (color online) Dominant components of the energy-dependent uncertainties for isotopic (²³⁵U and Pu combo) and total unfolded $\overline{\nu}_{e}$ energy spectra. The inset shows the details of the energy-dependent uncertainties for the total unfolded $\overline{\nu}_{e}$ energy spectrum. "Detector" uncertainty (detector response uncertainty) represents the contribution of the uncertainties from detection efficiency, energy scale difference between ADs, energy nonlinearity model, and IAV effect. Model (²³⁸U, ²⁴¹Pu) uncertainty represents uncertainties of the input model in the analysis to extract the isotopic (²³⁵U and Pu combo) prompt energy spectra. The contribution of the uncertainty from unfolding method takes into account the bias from the fluctuations of the input model.

The unfolded $ \overline{\nu}_{e} $ energy spectra are the $ \overline{\nu}_{e} $ energy distributions weighted by IBD cross section for one fission reaction and the subsequent beta decay reactions. They can thus be directly compared with theoretically predicted $ \overline{\nu}_{e} $ energy spectra of the IBD reaction. The unfolded $ \bar\nu_e $ energy spectra, both with and without IBD cross section weighting, are provided in the supplemental material.

As introduced above, the unfolded $ \overline{\nu}_{e} $ energy spectra of the IBD reaction can provide a model-independent input for other reactor antineutrino experiments. Furthermore, it can serve as an additional resource for validating theoretical models and standard nuclear databases. In this section, the application of the unfolded $ \overline{\nu}_{e} $ energy spectra is shown with detailed examples.

The generic unfolded $ \overline{\nu}_{e} $ energy spectra of the IBD reaction can be directly used by other reactor antineutrino experiments utilizing the IBD reaction for $ \overline{\nu}_{e} $ detection. Experiments which do not utilize the IBD reaction can in principle remove the IBD cross section to obtain original $ \overline{\nu}_{e} $ spectra from beta decay of fission fragments in the reactor. The antineutrino energy spectrum of the IBD reaction for an experiment with fission fractions ($ f_{235}^A $, $ f_{238}^A $, $ f_{239}^A $, $ f_{241}^A $) can be predicted as:

Here $ A $ represents an arbitrary reactor antineutrino experiment; $ \Delta f_i $ is the fission fraction difference between experiments for $ i $-th fissile isotope, $ \Delta f_i = f^{A}_i-f^{DB}_i $, where $ f^{DB} $ is the average effective fission fraction of the Daya Bay experiment (Sec. III); $ {{S}}_{\rm{total}} $, $ {{S}}_{235} $, and $ {{S}}_{\rm{combo}} $ are the total and isotopic unfolded $ \overline{\nu}_{e} $ energy spectra respectively. When the other experiment has different fission fractions than Daya Bay, small corrections based on the theoretical model inputs of $ {{S}}_{238} $ and $ {{S}}_{241} $ are necessary. For experiments detecting the $ \overline{\nu}_{e} $ from commercial LEU reactors, the Pu combo spectrum is recommended since it has smaller uncertainty and the fission fractions of ²³⁹Pu and ²⁴¹Pu are strongly correlated.

To estimate the uncertainty of the prediction using Eq. (8), one should consider the correlation between total and isotopic spectra. Equation (8) is converted to a matrix version in order to calculate the covariance matrix of the prediction more easily. A transformation matrix $ {{R}} $ is constructed as follows:

Here $ {{R}} $ is a $ 25\times 125 $ matrix, transforming 5 spectra into a prediction with different fission fractions; $ {{I}}_{a} $ means the $ a\times a $ identity matrix. The predicted $ \overline{\nu}_{e} $ energy spectrum corrected for fission fractions can be obtained by

In the latter example, $ {{S}}_{238} $ is chosen to be the Mueller ²³⁸U model [15] with a 15% bin-to-bin uncorrelated uncertainty and $ {{S}}_{241} $ is chosen to be the Huber ²⁴¹Pu model [16, 17] with a 10% bin-to-bin uncorrelated uncertainty. They both have little contribution to the total predicted spectrum because the fission fraction of ²³⁸U is relatively stable over time and the ratio between the fission fractions of ²³⁹Pu and ²⁴¹Pu is normally close to 0.183. Other models for ²³⁸U and ²⁴¹Pu with reasonable uncertainties could be chosen as well. The uncertainties of the prediction can be obtained by error propagation. The covariance matrix $ {{Cov}}_{\rm{pred}} $ of $ {{S}}_{\rm{pred}} $ is calculated by the formula:

Here $ {{Cov}}_{\rm{unfold}} $ is constructed based on the covariance matrix of the unfolded $ \overline{\nu}_{e} $ energy spectra from Eq. (5), with the correlation between total and isotopic spectra taken into account. In addition, $ {{Cov}}_{\rm{unfold}} $ also includes the covariance matrix of the models for ²³⁸U and ²⁴¹Pu, without consideration on their correlation with other spectra. In this way, the prediction of the $ \overline{\nu}_{e} $ energy spectrum for the other experiment can be obtained. Then the $ \overline{\nu}_{e} $ energy spectrum can be converted to the measurement of the prompt energy spectrum with detector response of the other experiment taken into consideration.

To demonstrate how to utilize the unfolded $ \overline{\nu}_{e} $ energy spectra and calculate the uncertainty of the prediction, examples are shown for predictions with different fission fractions than Daya Bay. The measurement of the prompt energy spectrum at EH1 AD1 is divided into 20 groups ordered by the ²³⁹Pu effective fission fraction in each week, as was done in Ref. [32]. The data in different groups represent different periods corresponding to different reactor burnup. Predictions of the prompt energy spectrum at EH1 AD1 in the earlier period ($ 5^{\rm{th}} $ group, $ f_{235}:f_{239}:f_{238}:f_{241} = 0.597: 0.278: 0.076: 0.049 $), intermediate period ($ 10^{\rm{th}} $ group, $ f_{235}:f_{239}:f_{238}:f_{241} = 0.567: $$ 0.301: 0.076: 0.056 $) and later period ($ 15^{\rm{th}} $ group, $f_{235}\,:\,f_{239}\,:\,f_{238}\,:\,f_{241} = 0.541\,:\, 0.322\,:\, 0.076\,:\, 0.061$) are performed. The predictions of the measured prompt energy spectrum based on unfolded $ \overline{\nu}_{e} $ energy spectra are obtained by applying the response matrix to transform the $ \overline{\nu}_{e} $ energy to the prompt energy of IBD events. Since the prompt energy spectra for Daya Bay are predicted, we can also do a prediction based on total and extracted isotopic prompt energy spectra $ {{S}}_{\rm{prompt}} $ in Sec. III directly for comparison. The results based on $ {{S}}_{\rm{prompt}} $ can be used to validate the prediction from unfolded $ \overline{\nu}_{e} $ energy spectra. The detector information (detector DAQ live time, detection efficiency, etc.), reactor information (time-varying thermal power, etc.), and survival probability from neutrino oscillation are taken into account in the prediction. The non-equilibrium effect is also considered in the prediction based on the reactor running status as was done in Ref. [32]. The predictions based on $ {{S}}_{\rm{unfold}} $ and $ {{S}}_{\rm{prompt}} $ are shown in Fig. 6. The measurements of the prompt energy spectrum with statistical uncertainty at EH1 AD1 are also shown in Fig. 6 for comparison with ~40000 IBD events excluding backgrounds in each period. The difference between the predictions based on $ {{S}}_{\rm{unfold}} $ and $ {{S}}_{\rm{prompt}} $ is small and well within uncertainties. The predictions are also consistent with the measurement in the 1 to 7 MeV energy region within statistical uncertainties. At lower or higher energy region the predictions are also consistent with the measurement when considering systematic uncertainties. This suggests the feasibility of utilizing unfolded $ \overline{\nu}_{e} $ energy spectra to get prediction for other reactor antineutrino experiments. The application to the experiments exposed to a larger range of fission fractions can be done in a similar way.

Figure 6.  (color online) (top) Prediction of the prompt energy spectra with different fission fractions at Daya Bay based on unfolded $\overline{\nu}_{e}$ energy spectra ${{S}}_{\rm{unfold}}$ or prompt energy spectra ${{S}}_{\rm{prompt}}$. The measurements with corresponding fission fractions from EH1 AD1 are shown for comparison with ~40000 IBD events in each period. Only statistical uncertainty on the measurement is shown in the plot. The legends on the top right represent the effective fission fractions of ²³⁵U, ²³⁹Pu, ²³⁸U, ²⁴¹Pu for each period respectively. (bottom) Comparison between prediction and measurement at Daya Bay. The error bars on the histogram only contain the uncertainties from prediction. Statistical uncertainty for the measurement is shown with the grey band. The difference between the prediction and the data are consistent within statistical uncertainty in the 1 to 7 MeV energy region.

Another application of the unfolded $ \overline{\nu}_{e} $ spectra is the comparison with theoretical reactor antineutrino models directly. Since the Huber-Mueller model provides the uncertainty and has a long history of comparison with measurements, this model is used as an example to compare with the unfolded $ \overline{\nu}_{e} $ energy spectra from Sec. IV. The detailed comparison between the Huber-Mueller model and the latest measurement of the prompt energy spectrum is shown in Ref. [32]. Here the focus is on the comparison based on $ \overline{\nu}_{e} $ energy spectra. The covariance matrix of the Huber-Mueller model ($ {{Cov}}^{\rm{HM}} $) is obtained by doing toy Monte Carlo based on the bin-to-bin correlated uncertainty and uncorrelated uncertainties provided in Refs. [15-17], and the calculation method is the same as in previous publications [31, 32]. The difference between the Huber-Mueller model ($ {{S}}^{\rm{HM}} $) and unfolded $ \overline{\nu}_{e} $ energy spectra ($ {{S}}^{\rm{unfold}} $) is evaluated based on the $ \chi^2 $ value defined as:

where $ {{S}}^{\rm{unfold}} $ represents the ²³⁵U, Pu combo, or total unfolded $ \overline{\nu}_{e} $ energy spectrum (each with 25 bins in this case), and $ {{Cov}}^{\rm{unfold}} $ represents their respective covariance matrix; $ {{S}}^{\rm{HM}} $ represents the prediction of the total or isotopic $ \overline{\nu}_{e} $ energy spectrum based on the Huber-Mueller model (each with 25 bins). $ {{Cov}}^{\rm{all}} $ is the total covariance matrix of $ {{S}}^{\rm{unfold}} $ and $ {{S}}^{\rm{HM}} $, which is the sum of $ {{Cov}}^{\rm{unfold}} $ and $ {{Cov}}^{\rm{HM}} $. The $ \chi^2 $ is equal to 57.2, 12.3, and 98.5 for ²³⁵U, Pu combo, and total $ \overline{\nu}_{e} $ energy spectra respectively. Since both the rate and spectral shape of the $ \overline{\nu}_{e} $ energy spectra are compared, the number of degrees of freedom in the above comparison is exactly the number of bins (25). Thus, the corresponding $ p $-values are $ 2.52\times 10^{-4} $, 0.984, and $ 1.14\times 10^{-10} $. The corresponding significance of deviation for the total $ \overline{\nu}_{e} $ energy spectrum is 6.4$ \sigma $, confirming the flux and spectral difference between measurement and the Huber-Mueller model.

Since the bias from input models in the unfolding method is included when evaluating the covariance matrix of the unfolded $ \overline{\nu}_{e} $ energy spectra, the sensitivity of the unfolded $ \overline{\nu}_{e} $ energy spectra is less powerful for reactor $ \overline{\nu}_{e} $ spectrum model discrimination than the original measurements of prompt energy spectrum. To achieve a more discriminating comparison with the Daya Bay data by avoiding the impact of unfolding bias, another method can be used to do the comparison. In principle, the unfolded $ \overline{\nu}_{e} $ energy spectrum with the Wiener-SVD unfolding method ($ {{S}}^{\rm{Wiener}} $) is obtained based on the following formula:

Here $ {{M}} $ is the response matrix of the detector; $ {{S}}^{\rm{prompt}} $ is the measurement from the detector; $ {{A}}_c $ is a smearing matrix used during the Wiener-SVD unfolding procedure, and it is constructed based on the signal to noise ratios in the effective frequency domain given expectations of signal and noise. The smearing matrix $ {{A}}_{c} $ suppresses the unphysical fluctuations but introduces biases as well during the unfolding procedure. The covariance matrix of the unfolded $ \overline{\nu}_{e} $ energy spectrum ($ {{Cov}}^{\rm{Wiener}} $) is calculated by error propagation from the Wiener-SVD unfolding package. $ {{Cov}}^{\rm{Wiener}} $ does not include the bias from the unfolding method, and it is different from $ {{Cov}}^{\rm{unfold}} $ based on toy Monte Carlo tests in Sec. IV. A more detailed description is found in Ref. [36]. The smearing matrices $ {{A}}_c $ are provided in the supplemental material.

To compare with the unfolded $ \overline{\nu}_{e} $ energy spectra, the theoretical $ \overline{\nu}_{e} $ energy spectra are smeared in the same way as the Wiener-SVD unfolding method: $ {{S}}^{\rm{smear}} = {{A}}_c{{S}}^{\rm{HM}} $. The covariance matrix of $ {{S}}^{\rm{smear}} $ is calculated by $ {{Cov}}^{\rm{smear}} = {{A}}_c{{Cov}}^{\rm{HM}}{{A}}_c^T $. In this way, the $ \chi^2 $ between the Huber-Mueller model and the unfolded $ \overline{\nu}_{e} $ energy spectra is calculated using the following formula:

where $ {{Cov}}^{\rm{sum}} = {{Cov}}^{\rm{Wiener}}+{{Cov}}^{\rm{smear}} $. Based on the above procedure, the $ \chi^2 $ is equal to 64.9, 22.3, and 108.8 for ²³⁵U, Pu combo, and total $ \overline{\nu}_{e} $ energy spectra, respectively. The corresponding $ p $-values are $2.11 \times 10^{-5} $, 0.618, $ 1.99\times 10^{-12} $, indicating greater discrepancy between the measurements and the model than considering the impact of bias in the unfolded $ \overline{\nu}_{e} $ energy spectra. This conclusion is further confirmed by a statistical check based on a dedicated toy Monte Carlo test.

The inclusion of the smearing matrix in the comparison with theoretical models maintains the discriminating sensitivity of the original data by avoiding the impact of bias in the unfolded $ \overline{\nu}_{e} $ energy spectra. This method can also be used to discriminate different hypotheses, with the same sensitivity of discrimination based on the original measurements of prompt energy spectra. To validate the feasibility of this method, two different hypotheses are considered: the "normalized Huber-Mueller model" with the same flux as the original measurements, and the "normalized distorted Huber-Mueller model" with an additional Gaussian (maximum: 10%, mean: 6 MeV, sigma: 0.6 MeV) added on the Huber-Mueller model, with the flux normalized as well. Different hypotheses are compared with the unfolded $ \overline{\nu}_{e} $ energy spectra based on the Eq. (14) (method 1), while the comparison between data and models based on the prompt energy spectra (method 2) is also done to verify the method 1. The $ \chi^2 $ based on ²³⁵U spectrum with method 1 (method 2) for the "normalized Huber-Mueller model" is 45.0 (45.1), while for the "normalized distorted Huber-Mueller model" the $ \chi^2 $ is 19.3 (19.3). To quantify the comparison between these two hypotheses, a frequentist approach is used based on 10000 toy experiments generated for the "normalized distorted Huber-Mueller model" hypothesis. A probability distribution function (PDF) is formed from the $ \Delta \chi^2 $ values for each toy experiment between the two models for both methods. The observed experimental $ \Delta \chi^2 $ value is compared to its respective PDF to determine quantitative preference of the observation for a particular model. The $ \Delta \chi^2 $ is found to be 25.7 (25.8) of the "normalized Huber-Mueller model" with respect to the "normalized distorted Huber-Mueller model" based on method 1 (method 2), correspond to $ p $-value of $ 7.1\times 10^{-3} $ ($ 8.1\times 10^{-3} $), indicating 2.7$ \sigma $ (2.6$ \sigma $) disfavoring of the "normalized Huber-Mueller model" hypothesis for ²³⁵U spectrum. Similarly, the $ p $-values are 0.276 (0.273) and 0.020 (0.023) for the Pu combo spectrum and the total spectrum with method 1 (method 2). The similar behavior in the discrimation of hypotheses between these two methods validates the feasibility of the comparison based on unfolded $ \overline{\nu}_{e} $ energy spectra using features of the Wiener-SVD unfolding approach.

A discussion of the possible usefulness of the unfolded $ \overline{\nu}_{e} $ energy spectra is listed below.

1. Generation of models of the $ \overline{\nu}_{e} $ flux and energy spectrum for different reactor types based on the measurement at Daya Bay:

1a) For experiments utilizing commercial LEU pressurized-water reactors [6, 11, 51], even with fission fractions different from those of Daya Bay, the prediction can be made at 2% precision with very little dependence on other isotopic $ \overline{\nu}_{e} $ flux models.

1b) For experiments utilizing highly-enriched uranium (HEU) research reactors [7, 8], the unfolded ²³⁵U $ \overline{\nu}_{e} $ energy spectrum can be used directly.

1c) For mixed-oxide (MOX) or other fuel/reactor types, which exhibit substantially different fission fractions than LEU reactor cores, some modest level of dependence on $ \overline{\nu}_{e} $ flux models of sub-dominant fission isotopes is expected.

1d) For reactor-based CE$ \nu $NS experiments at both HEU and LEU reactors [34, 35, 52], the CE$ \nu $NS cross section can be added on the generic $ \overline{\nu}_{e} $ energy spectra with the IBD cross section removed.

2. Comparison to other theoretical models in a variety of possible formats:

2a) The theoretical models can be compared with unfolded $ \overline{\nu}_{e} $ energy spectra directly using error matrices in the supplemental materials, at the expense of reduced precision from additional unfolding uncertainties.

2b) The comparison between theoretical models and unfolded $ \overline{\nu}_{e} $ energy spectra can be done with better precision using features of the Wiener-SVD unfolding method.

2c) The comparison can also be done by getting prediction of the prompt energy spectrum with better precision using Daya Bay's detector response matrix. This procedure has a long history, and some knowledge about the detector response of Daya Bay detectors is needed.

3. Antineutrino-based reactor monitoring applications:

3a) The unfolded $ \overline{\nu}_{e} $ energy spectra provide inputs for reactor power monitoring [53] or determination of the reactor fuel types [54]. Experiments with different detector response than Daya Bay can utilize the unfolded $ \overline{\nu}_{e} $ energy spectra to improve their precision.

3b) The unfolded $ \overline{\nu}_{e} $ energy spectra can also be used to reduce the likelihood of an undetected diversion of irradiated nuclear material [55, 56]. In addition, with inputs from isotopic $ \overline{\nu}_{e} $ energy spectra, the fission fractions of different isotopes for other reactors can be extracted.

Daya Bay has produced a significantly improved estimate of the $ \overline{\nu}_{e} $ energy spectrum from commercial nuclear reactors using 3.5 million inverse beta decay reactions. The estimate benefits from a reduction of the uncertainty in the energy response model from 1% to 0.5% [40] and a tripling in the IBD statistics compared to the previous results [31]. In addition, the measurement of prompt energy spectrum for all $ \bar\nu_e $ and for the $ \bar\nu_e $ resulting from the main fissioning isotopes, ²³⁵U, ²³⁹Pu and ²⁴¹Pu, along with a new unfolding method based on the Wiener-SVD method [36], enables estimates of the respective $ \bar\nu_e $ energy spectra and flux. A technique together with the relevant spectra and covariance matrices have been provided to allow a data-based prediction of the $ \bar\nu_e $ flux and energy spectra for arbitrary fission fractions from a commercial reactor. The unfolded $ \bar\nu_e $ energy spectra can also be directly compared with new theoretical predictions. To avoid bias from the unfolding method and to have better precision of the comparison between theoretical models and unfolded $ \overline{\nu}_{e} $ energy spectra, a new method is proposed using features of the Wiener-SVD unfolding method.

The $ \overline{\nu}_{e} $ energy spectrum from a nuclear reactor is a sum of the $ \overline{\nu}_{e} $'s from about 800 fission products, and it contains fine structures because of the different end-point energies (leading to sharp cutoffs) in the individual $ \overline{\nu}_{e} $ energy spectra from fission products with large fission yields [57, 58]. Recently an experiment with unprecedented energy resolution was proposed to search for such fine structures [51]. These fine structures may also be important for experiments performing precision measurements of neutrino oscillation [11]. Thus, we provide a prompt energy spectrum with fine energy bins.

To better match the energy resolution of the Daya Bay detectors, the bin width is reduced to 0.05 MeV from the previous choice of 0.25 MeV from 1 to 8 MeV. From 0.7 to 1 MeV the uncertainty is still dominated by systematic uncertainty, so only one bin is set in this energy range. The standard deviation of the energy resolution at Daya Bay (Eq. (1)) is shown in Fig. A1. The 0.05 MeV bin width is below the standard deviation of the energy resolution at all energies.

Figure A1.  (color online) The standard deviation of the energy resolution at Daya Bay. The 0.05 MeV bin width is shown in red dash line, which is below the standard deviation of the energy resolution at all energies.

To estimate the precision of the measured prompt energy spectrum, the statistical uncertainty and systematic uncertainties are evaluated based on the understanding of the detector discussed in Sec. II. The systematic uncertainties contain the contribution from detector (energy nonlinearity model, energy scale difference between ADs, IAV effect) and backgrounds (accidental background, cosmogenic ⁹Li and ⁸He beta-decays, fast neutrons, Am-C neutron sources, ¹³C($ \alpha $, n)¹⁶O reactions, and SNF background). The backgrounds were subtracted from the total measurement of IBD samples. Here the spectral uncertainty of the measurement is focused on while the flux information (~1.5% flux uncertainty) is not taken into account, thus the uncertainty of detection uncertainty has no impact on the spectral shape because of its energy independence. The total covariance matrix of the measured prompt energy spectrum is constructed as following:

where $ V^{\rm{corr}} $ ($ V^{\rm{uncorr}} $) represents bin-to-bin correlated (uncorrelated) uncertainty. $ V^{\rm{uncorr}} $ contains the contribution of statistical uncertainty, and they can be evaluated by constructing covariance matrix with non-zero diagonal elements only. On the other hand, $ V^{\rm{corr}} $ is evaluated by toy Monte Carlo as:

where $ N^{\rm{toy}} $ is the number of toy Monte Carlo samples, $ N_i^{\rm{ran(nom)}} $ is the random (nominal) predicted number of events at prompt energy bin $ i $. In each toy Monte Carlo test the total number of events in the random predicted spectra ($ \sum\nolimits_i N_i^{\rm{ran}} $) are normalized to the nominal predicted spectrum ($ \sum\nolimits_i N_i^{\rm{nom}} $). Figure A2 shows the fractional size of the diagonal elements of the covariance matrix, $ \sqrt{V_{ii}}/N^{\rm{nom}}_i $, for each component in each prompt energy bin. The elements of the correlation matrix, $ V_{ij}/\sqrt{V_{ii}V_{jj}} $, for the total uncertainty is also shown in Fig. A2. The relative uncertainty ($ \sqrt{V_{ii}}/N^{\rm{nom}}_i $) of the shape measurement is less than 1% from 2 to 5 MeV of prompt energy. The uncertainty is dominated by the statistical uncertainty in the 2 to 8 MeV energy region. The measured finely-binned prompt energy spectrum and its covariance matrix are provided in the supplemental materials.

Figure A2.  (color online) The fractional size of the diagonal elements of the covariance matrix, $\sqrt{V_{ii}}/N^{\rm{pred}}_i$, for each component in each prompt energy bin. Inset: the elements of the correlation matrix: $V_{ij}/\sqrt{V_{ii}V_{jj}}$.

Figure A3 shows the finely-binned prompt energy spectrum and the comparison with models normalized to the measured number of events. Here the Huber-Mueller model [15-17] from the conversion method and the SM2018 model [28] from the summation method are used for comparison. The formula of the IBD cross section from Ref. [59] is used to evaluate the detection probability of $ \overline{\nu}_{e} $ with different energy. The response matrix of Daya Bay detector introduced in Sec. II is used to map $ \overline{\nu}_{e} $ energy to reconstructed prompt energy of IBD events. The Huber-Mueller model is provided in 0.25 MeV bins below 8 MeV of $ \overline{\nu}_{e} $ energy, so the content within each 0.25 MeV bin is calculated using exponential interpolation and the uncertainty of the prediction is not shown. The comparison with the measurement above 7.2 MeV of prompt energy is not shown since the content above 8 MeV of $ \overline{\nu}_{e} $ energy for the Huber-Mueller model is from extrapolation. The SM2018 model is based on the summation method and lacks an uncertainty estimate. Non-equilibrium effects are considered for the Huber-Mueller model based on the correction from Ref. [15], which contributed as much as ~6% additional IBD candidates for specific fissile isotopes and energy bins. The non-equilibrium effect leads to ~0.71% (with relative uncertainty of 30%) more IBD candidates in total energy ranges. For the SM2018 model, the calculation with a 450-day irradiation duration ensures that the contribution from long-lived fission fragments reached equilibrium. Other information like detector DAQ livetime and detection efficiency, reactor power and working time, survival probability due to the $ \theta_{13} $ driven neutrino oscillation, is all taken into account in the calculation. A spectral distortion is observed through the whole energy region for both the Huber-Mueller model and the SM2018 model, revealing that there are still things missing in current antineutrino prediction models.

Figure A3.  (color online) (a) Measured prompt energy spectrum and the predictions from the Huber-Mueller (H-M) model and the SM2018 model (normalized to the number of measured events). The error bars on the data points (~0.6%) represent the square root of diagonal elements of the covariance matrix for the measurement, which contains both statistical and systematic uncertainties. (b) Ratio of the normalized predicted spectra and the measured prompt energy spectrum. Red (blue) histogram represents the comparison between measurement and the Huber-Mueller (SM2018) model. The error bars on the data points represent the uncertainties from measurement. (c) The continuity ($R_i = \dfrac{S_i}{S_{i+1}}$) of the measurement (black) and predictions from the Huber-Mueller model (red) and the SM2018 model (blue). (d) The continuity ($R_i^{\prime} = \dfrac{S_i+S_{i+2}-2S_{i+1}}{S_{i+1}}$) of the measurement (black) and predictions from the Huber-Mueller model (red) and the SM2018 model (blue).

With the finely-binned spectrum more studies can be carried out. In the following, we discuss one example in detail.

To study the fine structures and evaluate the continuity of the measurement, a method following Ref. [57] is used:

where $ S_i $ is the number of events in prompt energy bin $ i $ of the measurement. The $ R_i $ from a smooth spectrum without fine structures (e.g. Huber-Mueller model) was used to compare with the one from measurement based on the above method to evaluate the continuity.

To avoid the comparison with a smooth model when evaluating the continuity, we propose another method based on the following formula:

$ R_i $ and $ R_i^{\prime} $ have similar behavior of evaluating the continuity of the spectra, and $ R_i^{\prime} $ is more straightforward to reveal the unsmooth structures ($ R_i^{\prime}\neq0 $). The continuity of the measurement is shown in bottom two panels of Fig. A3. The evaluated continuities of the predictions from the Huber-Mueller model and the SM2018 model are quite similar, and the $ R_i^{\prime} $ of predictions is close to 0 at most of the energy region. To evaluate the sensitivity to fine structures at Daya Bay, the continuities of the $ \overline{\nu}_{e} $ energy spectra with different energy resolution are evaluated based on the SM2018 model, which are shown in Fig. A4. The SM2018 is based on summation method, with fine structures in the original $ \overline{\nu}_{e} $ energy spectrum. In the absence of resolution effects, $ R_i^{\prime} $ can reach $ 3\times 10^{-3} $ around 4.5 MeV on the $ \overline{\nu}_{e} $ energy spectra of the SM2018 model. After considering energy resolution of Daya Bay, $ R_i^{\prime} $ decreases to $ 5\times 10^{-4} $. For experiments with ~3%/$ \sqrt{E_{\rm{p}} [{\rm{MeV}}]} $ energy resolution [11, 51], the change on $ R_i^{\prime} $ is small when considering energy resolution smearing with respect to the original $ \overline{\nu}_{e} $ energy spectrum. This observation suggests that the Daya Bay measurement is not sensitive to fine structures calculated from current nuclear databases because of the finite energy resolution. To confirm this conclusion, a dedicated test was done by increasing the SM2018 prediction by ~10% below 4.5 MeV to mimic the shape cut-off caused by a decay branch. While $ R_i^{\prime} $ increases to 0.05 with respect to the original $ \overline{\nu}_{e} $ energy spectrum without energy resolution smearing, $ R_i^{\prime} $ increases to $ 2\times 10^{-3} $ only after taking into account the energy resolution smearing. Therefore, we conclude that the measurement at Daya Bay is not sensitive to the fine structures in the original $ \overline{\nu}_{e} $ energy spectrum. Nevertheless, we report the significance of deviation from continuity as $ \sigma^{\prime}_i = R_i^{\prime}/\Delta R_i^{\prime} $, where $ \Delta R_i^{\prime} $ is the uncertainty of $ R_i^{\prime} $. In addition, toy Monte Carlo tests are done by fluctuating the original prompt energy spectrum predicted by the SM2018 and Huber-Mueller models taking into account the measurement uncertainty. The average distribution of $ \sigma^{\prime}_i $ of 10000 prompt energy spectra from toy Monte Carlo tests is consistent with the distribution of the measurement. The $ \sigma^{\prime}_i $ distributions of the measurement and toy Monte Carlo tests are summarized in the Table A1. There are two adjacent points around 7.4 MeV with larger than 2$ \sigma $ deviation from continuity on the measurement. The distribution of $ \sigma^{\prime}_i $ using measured data is consistent with the normal distribution, suggesting that the unsmooth structures on the measurement are consistent with statistical fluctuations. No evidence of fine structures in reactor antineutrino spectrum based on Daya Bay measurement is found because of its finite energy resolution.

Figure A4.  (color online) (top panel) Original $\overline{\nu}_{e}$ energy spectrum weighted by the IBD cross section based on the SM2018 model, $\overline{\nu}_{e}$ energy spectrum with the energy resolution from Daya Bay, and $\overline{\nu}_{e}$ energy spectrum with 3%/$\sqrt{E_{\rm{p}} [{\rm{MeV}}]}$ energy resolution. Here $E_{\rm{p}}$ is the prompt energy, while the energy resolution of $\overline{\nu}_{e}$ energy $E_{\bar{\nu}_e}$ is calculated with $E_{\bar{\nu}_e}\approx E_{\rm{p}}+0.78\; {\rm{MeV}}$. (middle panel) The continuity ($R_i$) of the original $\overline{\nu}_{e}$ energy spectrum, the $\overline{\nu}_{e}$ energy spectrum with the energy resolution from Daya Bay, and the $\overline{\nu}_{e}$ energy spectrum with 3%/$\sqrt{E_{\rm{p}} [{\rm{MeV}}]}$ energy resolution. (bottom panel) The continuity ($R_i^{\prime}$) of the original $\overline{\nu}_{e}$ energy spectrum, the $\overline{\nu}_{e}$ energy spectrum with the energy resolution from Daya Bay, and the $\overline{\nu}_{e}$ energy spectrum with 3%/$\sqrt{E_{\rm{p}} [{\rm{MeV}}]}$ energy resolution.

	$<1\sigma$	$1\sim2\sigma$	$2\sim3\sigma$
Data	119	18	2
Toy Monte Carlo (average)	117.6	20.7	0.7

Table A1.  The significance of deviation from continuity ($\sigma^{\prime}_i$) based on the measurement and toy Monte Carlo tests.

We acknowledge Yellow River Engineering Consult- ing Co., Ltd., and China Railway 15th Bureau Group Co., Ltd., for building the underground laboratory. We are grateful for the ongoing cooperation from the China Guangdong Nuclear Power Group and China Light & Power Company.