Diagnosing the Reactor Antineutrino Anomaly with Global Antineutrino Flux Data

We have examined the impact of new Daya Bay, Double Chooz, and RENO measurements on global fits of reactor antineutrino flux data to a variety of hypotheses regarding the origin of the reactor antineutrino anomaly. In comparing RENO and Daya Bay measurements of inverse beta decay (IBD) yield versus $^{239}$Pu fission fraction, we find differing levels of precision in measurements of time-integrated yield and yield slope, but similar central values, leading to modestly enhanced isotopic IBD yield measurements in a joint fit of the two datasets. In the absence of sterile neutrino oscillations, global fits to all measurements now provide 3{\sigma} preference for incorrect modeling of specific fission isotopes over common mis-modeling of all beta-converted isotopes. If sterile neutrino oscillations are considered, global IBD yield fits provide no substantial preference between oscillation-including and oscillation-excluding hypotheses: hybrid models containing both sterile neutrino oscillations and incorrect $^{235}$U or $^{239}$Pu flux predictions are favored at only 1-2{\sigma} with respect to models where $^{235}$U, $^{238}$U, and $^{239}$Pu are assumed to be incorrectly predicted.


I. INTRODUCTION
Within the particle physics community, there remains enduring interest in the observed deficit of detected reactor antineutrino (ν e ) fluxes relative to the commonly-used conversion predictions [1,2]. This deficit, called the reactor antineutrino anomaly [3,4], has been observed over a wide range of baselines and reactor fission fractions.
It has been hypothesized that the observed deficit could be the result of oscillation of reactor ν e into unobservable sterile types via one or more new mass-squared splittings of the order of 1 eV 2 (see the review in Ref. [5]). Active-sterile oscillations should produce deficits in detected inverse beta decay (IBD) rates that are dependent on the baseline of the experiment and independent of the fission fractions in the observed reactors. An oscillation-based origin for the reactor antineutrino anomaly would have far-reaching experimental implications in neutrino physics, impacting the interpretation of prominent future long-baseline [6][7][8] and neutrinoless double beta decay experiments [9][10][11][12].
The reactor antineutrino anomaly could also be caused by inaccuracies in the beta-converted ν e flux models of the fission isotopes 235 U, 239 Pu, and 241 Pu and the ab initio model of 238 U [13][14][15][16][17][18][19][20][21][22]. When converting measured fission beta spectra from the BILL spectrometer into attendant antineutrino spectra for 235 U, 239 Pu, and 241 Pu [23][24][25][26], some inaccuracies could produce errors common to flux predictions of all isotopes: for example, non-consideration of important beta spectrum shape corrections [14]. On the other hand, other issues could produce errors specific to individual fission isotopes: for example, inconsistent calibration of neutron fluxes between different BILL beta spectrum measurements [27]. A model-based origin to the reactor anomaly should be reflected in a deficit in IBD detection rates that is not dependent on baseline and may or may not depend on the fission fractions of the experiment's reactor core.
Hybrids of these two origins have also been highlighted in the literature [28][29][30]. Such a scenario would produce dependencies of the measured IBD rate deficit on both fission fraction and baseline.
Recent studies have analyzed the global IBD yield dataset to provide measurements of individual isotopic IBD yields and to assess the consistency of these datasets with hypotheses regarding the sources of the reactor flux anomaly. Analyses including IBD yield measurements from highly 235 Uenriched (HEU) reactors provide indications that 235 U predictions could be incorrect [13,31], assuming the absence of active-sterile oscillations. Yield measurements from periods of differing observed fission fractions from Daya Bay, termed its 'flux evolution measurement', provide a distinct preference for incorrect 235 U predictions over sterile neutrino oscillations as the sole cause of the anomaly [32]. Meanwhile, combined analyses of both Daya Bay evolution and global IBD yield measurements investigating a wider variety of hypotheses has shown that best fits to these data are produced by a hybrid of both incorrect flux predictions and active-sterile neutrino oscillations [28]. Recent short-baseline measurements of the ratios of IBD energy spectra at different distances in the NEOS [33] and DANSS [34] experiments appear to also exhibit some preference for sterile neutrino oscillations [29], while other recent short-baseline measurements, PROSPECT [35] and STEREO [36], do not.
Recently, the community has seen the release of new results that are relevant to the investigation of these reactor anomaly hypotheses. In particular, the RENO collaboration has provided its first flux evolution measurement [37], and Daya Bay and Double Chooz have provided improved IBD yield measurements [38,39]. The goal of this paper is to provide a comparison between Daya Bay and RENO flux evolution results, and to determine the impact of recent flux results on the global preference for the hybrid model of active-sterile oscillations and incorrect flux predictions. We find that RENO and Daya Bay results provide a generally consistent picture of reactor flux evolution, but differ in their precision and their ability to differentiate between sterile-and model-related deficit hypotheses. We also show that the addition of RENO and the new absolute flux results enables some improvement in the precision of isotopic IBD yield measurements. Finally, the global flux fits are found to exhibit only marginal preference for oscillation-including hypotheses over oscillationexcluding ones.

II. EXPERIMENTAL INPUTS
Reactor antineutrino fluxes, sometimes reported experimentally as IBD yield, or the average flux times the IBD crosssection per fission, vary over time in a manner dependent on the fuel content of nearby reactor cores: where σ i is the IBD yield per fission for each parent fission isotope and F i (t) is the fission fraction of fission isotope i in the measured reactor core (i = 5, 8, 9, 1 for 235 U, 238 U, 239 Pu, and 241 Pu, respectively). A number of experiments have provided single measurements of time-integral IBD yields; in this case the measured IBD yieldσ f is dictated by the average fission fractions of nearby reactor cores over the measurement time period. Some experiments have also provided multiple IBD yield measurements from different time periods of varying fission fraction; given the high degree of detector stability exhibited in these experiments, these measurements are highly systematically correlated. Using the measured σ f values and corresponding fission fractions, one can attempt to determine IBD yields for the individual fission isotopes, σ i . For the global IBD yield fits presented here, we use as input the existing body of measurements from Ref. [29], with a few exceptions. This includes time-integral measurements from ILL [40,41], Savannah River [42], Krasnoyarsk [43][44][45], and Nucifer [46] at 235 U-burning HEU reactor cores, timeintegral measurements from conventional low-enriched cores from Gosgen [47], Rovno [48,49], Bugey-3 [50], Bugey-4 [51], Palo Verde [52], and Chooz [53], and the flux evolution measurement of Daya Bay [32].
In addition to these, we examine the inclusion of the new flux evolution measurement reported by the RENO collaboration [37], and the improved reactor flux measurements provided by Daya Bay [38] and Double Chooz [39]. RENO  To compare characteristics of the flux evolution data provided by Daya Bay and RENO, these results are overlaid in Fig. 1. It can be seen that the two set of measurements span roughly equivalent fission fraction ranges and show similarsized correlated uncertainty bands and uncorrelated statistical uncertainties. To further illustrate this comparison, we fit both experiments' data to linear functions as given in Ref. [32]: whereσ f is the time-integral IBD yield defined above, and Yield and slope values provided by the 235 U, 239 Pu, and 241 Pu predictions of Ref. [1] and the 238 U prediction of Ref [2] are also pictured in Fig. 1

III. IBD YIELD FITTING PROCEDURE
To determine the best-fit isotopic IBD yields from the experimentally-provided IBD yields and fission fractions, we use the following χ 2 definition: where the experimental inputs F i and σ f are those described above, with the indices a and b denoting the different experiments. The covariance matrix V ab describing the uncertainties of the measured σ f values is based on the uncertainties provided in Refs. [32,54], with alterations that take into account the new Daya Bay systematic uncertainty [38] and the uncertainties of the new RENO evolution data [37]. Reduced fully-correlated systematic uncertainties in the new Daya Bay flux measurement are propagated to the flux evolution dataset via a proportional subtraction from all on-and off-diagonal elements of Daya Bay's uncertainty covariance matrix. The covariance matrix for the RENO evolution measurement is formed from the quoted statistics along the diagonal and a flat contribution to all elements from their quoted 2.1% correlated systematic uncertainty. The isotopic IBD yields to be freely fitted are removed from the second expression, which constrains the isotopic IBD yields to the predicted yields σ th given in Refs. [1] for σ th 5,9,1 and Ref. [2] for σ th 8 , with a theoretical uncertainty matrix V HM given in Table 3 of Ref. [54]. The primary fit parameters are the ratios r i between the bestfit and predicted yields, and the neutrino mixing parameters sin 2 2ϑ ee ≡ 4|U e4 | 2 1 − |U e4 | 2 and ∆m 2 41 ≡ m 2 4 − m 2 1 , where U is the neutrino mixing matrix and m k is the mass of the massive neutrino ν k , which determine the averaged ν e oscillation survival probability P a ee for each experiment a in the 3+1 neutrino mixing scheme.
In order to test a range of scenarios regarding the origin of the reactor flux anomaly, we apply a variety of constraints on the fit parameters r i and P ee . The first set of hypotheses assumes flux predictions are the sole origin of the flux anomaly; this is achieved by adding the constraint P ee = 1, as well as the following additional constraints on various isotopes' yields: • 235: Constrain all r i except r 5 .
Additional hypotheses including free fits of sterile neutrino oscillation parameters are also considered. These scenarios correspond to the reactor flux anomaly being caused by sterile neutrino oscillations alone, or a hybrid combination of oscillations and incorrect flux modelling.
• OSC: Constrain all r i and freely fit P ee .
Rather than present best fits for all hypotheses and data combinations, we will highlight noteworthy results for each considered data combination.

IV. COMPARISON OF RENO AND DAYA BAY RESULTS
The allowed regions for the IBD yields of 235 U and 239 Pu in the absence of oscillations (235+239 hypothesis) are shown in Fig. 2 for the new Daya Bay and RENO datasets, with best-fit values also overviewed in Tabs. I and II. Measured 235 U IBD yields with respect to predictions are 0.926±0.016 and 0.913±0.027 for Daya Bay and RENO, respectively, while values for 239 Pu are 0.981±0.036 and 0.957±0.054. The improvement in Daya Bay's detection efficiency has improved its IBD yield measurements [32]: errors on 235 U and 239 Pu yields have been reduced by 0.9% and 1.0%, respectively, with respect to Ref. [28], which uses identical χ 2 definitions and theoretical IBD yield uncertainties.
RENO's best-fit 235 U and 239 Pu IBD yields are quite similar to those obtained from the new Daya Bay dataset. Both experiments observe a substantial difference in 235 U yield relative to Huber-Mueller, but not in 239 Pu yield. The result of a combined fit of the Daya Bay and RENO datasets to the 235+239 hypothesis is also pictured in Fig. 2. The combined fit produces minor improvements in precision for r 5 (from 1.6% to 1.5%) and r 9 (from 3.6% to 3.2%) over those obtained by Daya Bay alone. We note that we have not considered the 235+238+239 hypothesis here, as the combined RENO and Daya Bay data are not sufficient to constrain it.
Results from fits of other hypotheses to the RENO and Daya Bay datasets are also overviewed in Tabs. I and II 1 . While RENO data prefers hypotheses involving incorrect fluxes to the pure OSC hypothesis, all hypotheses in Tabs. I 1 In the analysis of Daya Bay and RENO evolution data alone we consider the averaged survival probability Pee = 1 − sin 2 2ϑee/2, because the source-detector distance is much larger than the oscillation length for  and II save 239 exhibit a χ 2 min within 1.8 of the overall minimum. Meanwhile, Daya Bay data shows substantial preference for incorrect flux modelling: for example, a ∆χ 2 min of 5.7 is seen between OSC and 235 models. As discussed above, this difference in model discrimination power is due to the differences in experimental uncertainties between experiments, as opposed to substantial differences in the central values of best-fit parameters.
Given the similarity of the best-fit parameters from the two datasets, a combined fit yields enhancements in the preferences against the OSC model. The ∆χ 2 min between 235 and OSC increased to from 5.7 for Daya Bay to 7.3 in the combined fit. Using a frequentist Monte Carlo statistical analysis [19,28], this corresponds to a change in the preference of the 235 model against the OSC model from 2.6σ for Daya Bay to 2.9σ in the combined fit of RENO and Daya Bay.

V. GLOBAL FLUX FITS
We now turn to global fits of all time-integral and evolution IBD yield measurements. A comparison of global flux fit results for various oscillation-including or -excluding hypotheses introduced above are summarized in Tab. III.
In the absence of oscillations, the allowed regions for the IBD yields of 235 U, 238 U, and 239 Pu (235+239+238 hypothesis) are pictured in Fig, 3, along with the previouslydetermined best-fit values [19]. The best-fit r values r 5 = 0.952±0.014, r 8 = 0.672±0.135, and r 9 = 1.042±0.046 are obtained for 235 U, 238 U, and 239 Pu, respectively. The addition of improved Daya Bay, RENO and Double Chooz datasets has modestly improved the combined IBD yield constraints for 238 U and 239 Pu, that before were given by r 8 = 0.695±0.163 and r 9 = 1.034 ± 0.064 [19], whereas the uncertainty of the 235 U IBD yield is practically unchanged. As in previous fits neglecting oscillations, measured IBD yields for 235 U and 238 U disagree substantially with predicted central values, now at the 3.6σ and 2.4σ level, while the yield for 239 Pu remains consistent with its predicted value.
Comparing the different oscillation-excluding fits in Tab. III, we find substantially higher χ 2 values provided by the 235=239=241+238 and 239 hypotheses. The former hy- pothesis corresponds to a common inaccuracy being present in all beta-conversion antineutrino flux predictions, while the latter corresponds to 239 Pu being the sole cause of the reactor flux anomaly. The overall worst fit is provided by the 239 hypothesis, which is in tension with existing flux constraints from 235 U-burning HEU reactors, as well as with the Daya Bay and RENO evolution datasets. If the 'common inaccuracy' 235=239=241+238 hypothesis is quantitatively compared with the hypothesis of uncorrelated inaccuracies between isotopes (235+238+239 hypothesis) using a frequentist Monte Carlo statistical approach, the 8.7 ∆χ 2 between models corresponds to 3.0σ preference for the latter hypothesis. Thus, if sterile neutrinos do not contribute to the reactor flux anomaly, the global IBD yield data favors model inaccuracies that are particular to specific fission isotopes. This conclusion is supported by recent work suggesting inconsistent calibration of neutron fluxes between fission beta spectrum measurements made by the BILL spectrometer [27]. In particular, the IBD yield data points to incorrect calibration of results for 235 U as well as incorrect ab initio prediction of the 238 U flux.
Considering now all possible oscillation-including or excluding hypotheses, we see that lowest χ 2 min values are delivered by the hybrid 235+OSC and 239+OSC models, as well by as the 235+238+239 model discussed above. Using the frequentist statistical approach, we find that the two hybrid models 235+OSC and 239+OSC are preferred at 1.1σ and 1.8σ, respectively, to the most-preferred oscillation-excluding hypothesis 235+238+239. The global flux fit thus does not provide definitive preference for or against the existence of sterile neutrino oscillations. There is only a small improvement in comparing the two hybrid models 235+OSC and 239+OSC with the oscillation-excluding hypothesis 235+239, which is disfavored at 1.6σ and 2.3σ, respectively, or with the oscillation-excluding hypothesis 235, which is disliked by the data at 1.4σ and 2.0σ, respectively. We also note that the 239+OSC hypothesis is preferred to the 235+OSC hypothesis by only 1.5σ.
To allow examination of oscillation-including hypotheses in more detail, the allowed regions in the plane of the oscillation parameters ∆m 2 41 and sin 2 2ϑ ee for the three oscillationincluding hypotheses are pictured in Fig. 4. These hypotheses involving either zero or one incorrectly predicted fluxes provide similar best fit regions, except that these regions are shifted relatively between one another in sin 2 2ϑ ee space: 235+OSC and 239+OSC exhibit the lowest and highest bestfit sin 2 2ϑ ee values, respectively. This indicates that the null-oscillation IBD energy spectrum ratio results reported by DANSS, NEOS, PROSPECT, and STEREO are likely to have the most substantial impact on the oscillation parameter space suggested by the 239+OSC hypothesis. We also note that although all three oscillation-including hypotheses fit the data well, the hybrid 235+OSC and 239+OSC hypotheses are preferred to the pure oscillation hypothesis (OSC) by 1.9σ and 2.5σ, respectively.

VI. SUMMARY
We have performed global fits of the complete set of reactor ν e flux data, including new flux evolution data from RENO and new time-integrated flux measurements from Daya Bay and Double Chooz. We find that the RENO and Daya Bay flux evolution datasets are similar in their preferred central values of absolute IBD yield and yield slopes, but differ in their level of precision, due to lower statistical and systematic uncertainties and wider fission fraction ranges provided by Daya Bay. A joint fit of the two datasets leads to an increased 2.9σ preference for incorrect 235 U predictions over sterile neutrino oscillations as the sole source of the reactor antineutrino flux anomaly. A global fit of these two evolution results and all time-integrated ν e flux measurements produces improved IBD yield constraints over those reported in previous publications. We find that all ν e flux data, in the absence of oscillations, now disfavor a common inaccuracy among all beta conversion predictions at 3.0σ confidence level. We also find that flux data, alone, currently does not provide a sizable preference for oscillation-including or oscillation-excluding hypotheses.