Diagnosing the Reactor Antineutrino Anomaly with Global Antineutrino Flux Data

TL;DR

This paper analyzes recent reactor antineutrino measurements to understand the anomaly, finding that data favor incorrect isotope flux models over sterile neutrino explanations, with no strong preference for oscillations.

Contribution

It provides a comprehensive global fit of new antineutrino flux data, evaluating hypotheses about the anomaly and the role of sterile neutrinos and isotope flux modeling.

Findings

Global fits favor isotope flux mis-modeling over sterile neutrino oscillations.

No significant preference for sterile neutrino oscillations when considering new data.

Hybrid models with both oscillations and flux inaccuracies are only mildly favored.

Abstract

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.

Tables

Table 1. Table 1: Results of fitting Daya Bay and RENO flux evolution datasets with three oscillation-excluding hypotheses in Section  III regarding the origin of the reactor anomaly. Each hypothesis name denotes the unconstrained isotopic IBD yields in the fit. The best-fit values and 1 σ 𝜎 \sigma ranges are given for unconstrained parameters, while the parenthetical values denote the best-fit values of the constrained fit parameters.

	235			239			235+239
	Daya Bay	RENO	DB+RENO	Daya Bay	RENO	DB+RENO	Daya Bay	RENO	DB+RENO
${\chi }_{\text{min}}^2$	$3.8$	$5.2$	$9.0$	$10.2$	$10.5$	$16.1$	$3.6$	$4.8$	$8.7$
NDF	$7$	$7$	$15$	$7$	$7$	$15$	$6$	$6$	$14$
GoF	$80\%$	$64\%$	$88\%$	$18\%$	$16\%$	$38\%$	$73\%$	$58\%$	$85\%$
$r_5$	$0.927\pm 0.016$	$0.917\pm 0.026$	$0.925\pm 0.015$	$\left(0.950\right)$	$\left(0.961\right)$	$\left(0.945\right)$	$0.926\pm 0.016$	$0.913\pm 0.027$	$0.923\pm 0.015$
$r_8$	$\left(0.988\right)$	$\left(0.985\right)$	$\left(0.985\right)$	$\left(0.935\right)$	$\left(0.949\right)$	$\left(0.929\right)$	$\left(1.001\right)$	$\left(1.000\right)$	$\left(1.003\right)$
$r_9$	$\left(0.993\right)$	$\left(0.989\right)$	$\left(0.990\right)$	$1.015\pm 0.034$	$0.992\pm 0.052$	$1.014\pm 0.030$	$0.981\pm 0.036$	$0.957\pm 0.054$	$0.975\pm 0.032$
$r_1$	$\left(0.994\right)$	$\left(0.990\right)$	$\left(0.991\right)$	$\left(0.947\right)$	$\left(0.959\right)$	$\left(0.942\right)$	$\left(1.001\right)$	$\left(1.000\right)$	$\left(1.000\right)$

Table 2. Table 2: As Tab.  1 , for the three oscillation-including hypotheses in Section  III .

	OSC			235+OSC			239+OSC
	Daya Bay	RENO	DB+RENO	Daya Bay	RENO	DB+RENO	Daya Bay	RENO	DB+RENO
${\chi }_{\text{min}}^2$	$9.5$	$6.5$	$16.3$	$3.7$	$4.7$	$8.8$	$3.8$	$4.7$	$8.8$
NDF	$7$	$7$	$15$	$6$	$6$	$14$	$6$	$6$	$14$
GoF	$22\%$	$48\%$	$37\%$	$72\%$	$58\%$	$85\%$	$70\%$	$59\%$	$84\%$
$P_{ee}$	$0.944\pm 0.025$	$0.935\pm 0.030$	$0.939\pm 0.024$	$0.988\pm 0.021$	$0.972\pm 0.034$	$0.986\pm 0.022$	$0.932\pm 0.025$	$0.921\pm 0.031$	$0.929\pm 0.024$
$r_5$	$\left(1.003\right)$	$\left(1.001\right)$	$\left(1.004\right)$	$0.937\pm 0.030$	$0.943\pm 0.047$	$0.938\pm 0.029$	$\left(1.000\right)$	$\left(1.000\right)$	$\left(1.000\right)$
$r_8$	$\left(1.037\right)$	$\left(1.009\right)$	$\left(1.045\right)$	$\left(1.004\right)$	$\left(1.000\right)$	$\left(1.001\right)$	$\left(1.003\right)$	$\left(1.000\right)$	$\left(1.002\right)$
$r_9$	$\left(1.006\right)$	$\left(1.002\right)$	$\left(1.007\right)$	$\left(1.000\right)$	$\left(1.000\right)$	$\left(0.999\right)$	$1.094\pm 0.048$	$1.088\pm 0.070$	$1.094\pm 0.043$
$r_1$	$\left(1.003\right)$	$\left(1.002\right)$	$\left(1.004\right)$	$\left(1.000\right)$	$\left(1.000\right)$	$\left(1.000\right)$	$\left(1.000\right)$	$\left(1.000\right)$	$\left(1.000\right)$

Table 3. Table 3: Results of fitting all time-integral and evolution flux measurements with hypotheses regarding the origin of the reactor anomaly. Each hypothesis name denotes the unconstrained parameters in the fit. The best-fit values and 1 σ 𝜎 \sigma ranges are given for unconstrained parameters, while the parenthetical values denote the best-fit values of the constrained fit parameters.

	235	239	235+239	235+238+239	235=239=241+238	OSC	235+OSC	239+OSC
${\chi }_{\text{min}}^2$	$34.6$	$41.6$	$34.1$	$29.9$	$38.6$	$33.1$	$29.5$	$26.9$
NDF	$39$	$39$	$38$	$37$	$38$	$38$	$37$	$37$
GoF	$67\%$	$36\%$	$65\%$	$79\%$	$44\%$	$69\%$	$80\%$	$89\%$
$r_5$	$0.933\pm 0.010$	$\left(0.941\right)$	$0.932\pm 0.009$	$0.952\pm 0.014$	$0.941\pm 0.013$	$\left(1.014\right)$	$0.984\pm 0.025$	$\left(1.014\right)$
$r_8$	$\left(0.890\right)$	$\left(0.868\right)$	$\left(0.914\right)$	$0.672\pm 0.135$	$0.926\pm 0.096$	$\left(1.021\right)$	$\left(0.969\right)$	$\left(0.956\right)$
$r_9$	$\left(0.987\right)$	$0.997\pm 0.029$	$0.969\pm 0.030$	$1.042\pm 0.046$	$0.944\pm 0.015$	$\left(1.019\right)$	$\left(1.026\right)$	$1.099\pm 0.040$
$r_1$	$\left(0.989\right)$	$\left(0.938\right)$	$\left(1.003\right)$	$\left(1.001\right)$	$0.942\pm 0.013$	$\left(1.015\right)$	$\left(1.024\right)$	$\left(1.015\right)$
$\Delta m_{41}^2$						${0.49}_{-0.03}^{+0.02}$	${0.48}_{-0.03}^{+0.05}$	$0.49\pm 0.02$
${\sin }^{2}2{\vartheta }_{ee}$						$0.15\pm 0.04$	${0.10}_{-0.04}^{+0.05}$	$0.16\pm 0.04$

Equations

$${\sigma }_f(t)=\sum _i{F}_i(t){\sigma }_i,$$

$${\sigma }_f(F_{239})={\bar{\sigma }}_f+\frac{d{\sigma }_f}{d{F}_{239}}(F_{239}-{\overline{F}}_{239}),$$

$$\begin{gathered}\chi^{2}=\sum_{a,b}\bigg{(}\sigma_{f,a}-P_{ee}^{a}\sum_{i}r_{i}F_{i,a}\sigma_{i}\bigg{)}\textrm{V}^{-1}_{ab}\\ \times\bigg{(}\sigma_{f,b}-P_{ee}^{b}\sum_{i}r_{i}F_{i,b}\sigma_{i}\bigg{)}\\ +\sum_{i,j}(\sigma^{th}_{i}-\sigma_{i})\textrm{V}^{-1}_{\textrm{HM},ij}(\sigma^{th}_{j}-\sigma_{j}),\end{gathered}$$

Full text

††preprint: APS/123-QED

I Introduction

Within the particle physics community, there remains enduring interest in the observed deficit of detected reactor antineutrino ($\overline{\nu}_{e}$) fluxes relative to the commonly-used conversion predictions Huber (2011); Mueller et al. (2011). This deficit, called the reactor antineutrino anomaly Mention et al. (2011); Abazajian et al. (2012), 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 $\overline{\nu}_{e}$ into unobservable sterile types via one or more new mass-squared splittings of the order of 1 eV2 (see the review in Ref. Gariazzo et al. (2016)). 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 de Gouvêa et al. (2015); Klop and Palazzo (2015); Gandhi et al. (2015) and neutrinoless double beta decay experiments Barry et al. (2011); Li and Liu (2012); Girardi et al. (2013); Giunti and Zavanin (2015).

The reactor antineutrino anomaly could also be caused by inaccuracies in the beta-converted $\overline{\nu}_{e}$ flux models of the fission isotopes 235U, 239Pu, and 241Pu and the ab initio model of 238U  Giunti (2017a); A. Hayes, J. Friar, G. Garvey, G. Jungman, and G. Jonkmans (2014); Hayes et al. (2015); Hayes and Vogel (2016); Dwyer and Langford (2015); Hayes et al. (2017); Gebre et al. (2018); Littlejohn et al. (2018); Wang and Hayes (2017); Wang et al. (2016). When converting measured fission beta spectra from the BILL spectrometer into attendant antineutrino spectra for 235U, 239Pu, and 241Pu Mampe et al. (1978); Schreckenbach et al. (1985); Von Feilitzsch et al. (1982); Hahn et al. (1989), some inaccuracies could produce errors common to flux predictions of all isotopes: for example, non-consideration of important beta spectrum shape corrections A. Hayes, J. Friar, G. Garvey, G. Jungman, and G. Jonkmans (2014). 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 Onillon (2018). 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 Giunti et al. (2017); Gariazzo et al. (2018); Dentler et al. (2017). 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 235U-enriched (HEU) reactors provide indications that 235U predictions could be incorrect Giunti (2017a, b), 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 235U predictions over sterile neutrino oscillations as the sole cause of the anomaly An et al. (2017). 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 Giunti et al. (2017). Recent short-baseline measurements of the ratios of IBD energy spectra at different distances in the NEOS Ko et al. (2017) and DANSS Alekseev et al. (2018) experiments appear to also exhibit some preference for sterile neutrino oscillations Gariazzo et al. (2018), while other recent short-baseline measurements, PROSPECT Ashenfelter et al. (2018) and STEREO Lhuillier (2015), 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 Bak et al. (2018), and Daya Bay and Double Chooz have provided improved IBD yield measurements An et al. (2018); Bezerra (2018). 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 oscillation-excluding ones.

II Experimental Inputs

Reactor antineutrino fluxes, sometimes reported experimentally as IBD yield, or the average flux times the IBD cross-section per fission, vary over time in a manner dependent on the fuel content of nearby reactor cores:

[TABLE]

where $\sigma_{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 235U, 238U, 239Pu, and 241Pu, respectively). A number of experiments have provided single measurements of time-integral IBD yields; in this case the measured IBD yield $\bar{\sigma}_{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 $\sigma_{f}$ values and corresponding fission fractions, one can attempt to determine IBD yields for the individual fission isotopes, $\sigma_{i}$.

For the global IBD yield fits presented here, we use as input the existing body of measurements from Ref. Gariazzo et al. (2018), with a few exceptions. This includes time-integral measurements from ILL Kwon et al. (1981); Hoummada et al. (1995), Savannah River Greenwood et al. (1996), Krasnoyarsk Vidyakin et al. (1987, 1990, 1994), and Nucifer Boireau et al. (2016) at 235U-burning HEU reactor cores, time-integral measurements from conventional low-enriched cores from Gosgen Zacek et al. (1986), Rovno Afonin et al. (1988); Kuvshinnikov et al. (1991), Bugey-3 Achkar et al. (1995), Bugey-4 Declais et al. (1994), Palo Verde Boehm et al. (2001), and Chooz Apollonio et al. (2003), and the flux evolution measurement of Daya Bay An et al. (2017).

In addition to these, we examine the inclusion of the new flux evolution measurement reported by the RENO collaboration Bak et al. (2018), and the improved reactor flux measurements provided by Daya Bay An et al. (2018) and Double Chooz Bezerra (2018). RENO’s new evolution result provides highly-correlated flux measurements at eight different fission fraction values, while Daya Bay and Double Chooz flux measurements have been improved in precision to the 1.5% and 1.0% level. To account for a 0.3% shift in the new Daya Bay time-integral flux with respect to the normalization of its older flux evolution result, all Daya Bay flux evolution data points are coherently shifted by this amount in our analysis.

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 similar-sized correlated uncertainty bands and uncorrelated statistical uncertainties. To further illustrate this comparison, we fit both experiments’ data to linear functions as given in Ref. An et al. (2017):

[TABLE]

where $\bar{\sigma}_{f}$ is the time-integral IBD yield defined above, and $\frac{d\sigma_{f}}{dF_{239}}$ is the change in IBD yield per unit change in 239Pu fission fraction $F_{239}$. Measured time integrated yields $\bar{\sigma}_{f}$ are (5.93$\pm$0.09) $\times 10^{-43}$ cm2/fission and (5.84$\pm$0.12) $\times 10^{-43}$ cm2/fission for Daya Bay and RENO, respectively, while measured slopes $\frac{d\sigma_{f}}{dF_{239}}$ are (-1.87$\pm$0.18) $\times 10^{-43}$ cm2/fission and (-1.93$\pm$0.29) $\times 10^{-43}$ cm2/fission. Yields and slopes are consistent within 1$\sigma$ between the two experiments. The modestly larger Daya Bay fission fraction range and smaller Daya Bay correlated (1.5% versus 2.0%) and uncorrelated ($\sim$0.1% versus $\sim$0.2% per data point) uncertainties produce smaller uncertainties in its measurement of the time-integral yield and the slope.

Yield and slope values provided by the 235U, 239Pu, and 241Pu predictions of Ref. Huber (2011) and the 238U prediction of Ref Mueller et al. (2011) are also pictured in Fig. 1. The level of disagreement of Daya Bay and RENO with predicted time-integral yields, 2.0$\sigma$ and 2.2$\sigma$, respectively, are relatively similar. For relative slopes, $\frac{1}{\bar{\sigma}_{f}}\frac{d\sigma_{f}}{dF_{239}}$, Daya Bay and RENO show differing levels of consistency with predictions (3.1$\sigma$ versus 1.8$\sigma$, respectively), despite similar central values.

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 $\chi^{2}$ definition:

[TABLE]

where the experimental inputs $F_{i}$ and $\sigma_{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 $\sigma_{f}$ values is based on the uncertainties provided in Refs. Gariazzo et al. (2017); An et al. (2017), with alterations that take into account the new Daya Bay systematic uncertainty An et al. (2018) and the uncertainties of the new RENO evolution data Bak et al. (2018). 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 $\sigma^{th}$ given in Refs. Huber (2011) for $\sigma^{th}_{5,9,1}$ and Ref. Mueller et al. (2011) for $\sigma^{th}_{8}$, with a theoretical uncertainty matrix $V_{\textrm{HM}}$ given in Table 3 of Ref. Gariazzo et al. (2017). The primary fit parameters are the ratios $r_{i}$ between the best-fit and predicted yields, and the neutrino mixing parameters $\sin^{2}2\vartheta_{ee}\equiv 4|U_{e4}|^{2}\left(1-|U_{e4}|^{2}\right)$ and $\Delta m^{2}_{41}\equiv m_{4}^{2}-m_{1}^{2}$, where $U$ is the neutrino mixing matrix and $m_{k}$ is the mass of the massive neutrino $\nu_{k}$, which determine the averaged $\overline{\nu}_{e}$ oscillation survival probability $P_{ee}^{a}$ 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}$.

• •

239: Constrain all $r_{i}$ except $r_{9}$.

• •

235+239: Constrain only $r_{8}$ and $r_{1}$.

• •

235+238+239: Constrain only $r_{1}$.

• •

235=239=241+238: Require common scaling of $r_{5,9,1}$; allow variations between $r_{5,9,1}$ within the Huber uncorrelated uncertainties.

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}$.

• •

239+OSC: Same, but removing constraint on $r_{5}$.

• •

235+OSC: Same, but removing constraint on $r_{9}$.

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 235U and 239Pu 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. 1 and 2. Measured 235U IBD yields with respect to predictions are 0.926$\pm$0.016 and 0.913$\pm$0.027 for Daya Bay and RENO, respectively, while values for 239Pu are 0.981$\pm$0.036 and 0.957$\pm$0.054. The improvement in Daya Bay’s detection efficiency has improved its IBD yield measurements An et al. (2017): errors on 235U and 239Pu yields have been reduced by 0.9% and 1.0%, respectively, with respect to Ref. Giunti et al. (2017), which uses identical $\chi^{2}$ definitions and theoretical IBD yield uncertainties.

RENO’s best-fit 235U and 239Pu IBD yields are quite similar to those obtained from the new Daya Bay dataset. Both experiments observe a substantial difference in 235U yield relative to Huber-Mueller, but not in 239Pu 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. 1 and 2111 In the analysis of Daya Bay and RENO evolution data alone we consider the averaged survival probability $P_{ee}=1-\sin^{2}2\vartheta_{ee}/2$, because the source-detector distance is much larger than the oscillation length for $\Delta m^{2}_{41}\gtrsim 0.1\,\text{eV}^{2}$. . While RENO data prefers hypotheses involving incorrect fluxes to the pure OSC hypothesis, all hypotheses in Tabs. 1 and 2 save 239 exhibit a $\chi^{2}_{\text{min}}$ within 1.8 of the overall minimum. Meanwhile, Daya Bay data shows substantial preference for incorrect flux modelling: for example, a $\Delta\chi^{2}_{\text{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 $\Delta\chi^{2}_{\text{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 Giunti et al. (2017); Gebre et al. (2018), this corresponds to a change in the preference of the 235 model against the OSC model from $2.6\sigma$ for Daya Bay to $2.9\sigma$ 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. 3.

In the absence of oscillations, the allowed regions for the IBD yields of 235U, 238U, and 239Pu (235+239+238 hypothesis) are pictured in Fig, 3, along with the previously-determined best-fit values Gebre et al. (2018). The best-fit $r$ values $r_{5}=0.952\pm 0.014$, $r_{8}=0.672\pm 0.135$, and $r_{9}=1.042\pm 0.046$ are obtained for 235U, 238U, and 239Pu, respectively. The addition of improved Daya Bay, RENO and Double Chooz datasets has modestly improved the combined IBD yield constraints for 238U and 239Pu, that before were given by $r_{8}=0.695\pm 0.163$ and $r_{9}=1.034\pm 0.064$ Gebre et al. (2018), whereas the uncertainty of the 235U IBD yield is practically unchanged. As in previous fits neglecting oscillations, measured IBD yields for 235U and 238U disagree substantially with predicted central values, now at the $3.6\sigma$ and $2.4\sigma$ level, while the yield for 239Pu remains consistent with its predicted value.

Comparing the different oscillation-excluding fits in Tab. 3, we find substantially higher $\chi^{2}$ values provided by the 235=239=241+238 and 239 hypotheses. The former hypothesis corresponds to a common inaccuracy being present in all beta-conversion antineutrino flux predictions, while the latter corresponds to 239Pu 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 235U-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 $\Delta\chi^{2}$ between models corresponds to 3.0$\sigma$ 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 Onillon (2018). In particular, the IBD yield data points to incorrect calibration of results for 235U as well as incorrect ab initio prediction of the 238U flux.

Considering now all possible oscillation-including or excluding hypotheses, we see that lowest $\chi^{2}_{\text{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\sigma$ and $1.8\sigma$, 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\sigma$ and $2.3\sigma$, respectively, or with the oscillation-excluding hypothesis 235, which is disliked by the data at $1.4\sigma$ and $2.0\sigma$, respectively. We also note that the 239+OSC hypothesis is preferred to the 235+OSC hypothesis by only $1.5\sigma$.

To allow examination of oscillation-including hypotheses in more detail, the allowed regions in the plane of the oscillation parameters $\Delta m^{2}_{41}$ and $\sin^{2}2\vartheta_{ee}$ for the three oscillation-including 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\vartheta_{ee}$ space: 235+OSC and 239+OSC exhibit the lowest and highest best-fit $\sin^{2}2\vartheta_{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\sigma$ and $2.5\sigma$, respectively.

VI Summary

We have performed global fits of the complete set of reactor $\overline{\nu}_{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$\sigma$ preference for incorrect 235U 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 $\overline{\nu}_{e}$ flux measurements produces improved IBD yield constraints over those reported in previous publications. We find that all $\overline{\nu}_{e}$ flux data, in the absence of oscillations, now disfavor a common inaccuracy among all beta conversion predictions at 3.0$\sigma$ confidence level. We also find that flux data, alone, currently does not provide a sizable preference for oscillation-including or oscillation-excluding hypotheses.

Acknowledgements.

We would like to thank Yonas Gebre, Soo-Bong Kim, Marco Laveder, and Xianyi Zhang for useful discussions. The work of Y.F. Li was supported in part by the National Natural Science Foundation of China under Grant No. 11835013, by the Strategic Priority Research Program of the Chinese Academy of Sciences under Grant No. XDA10010100, and by the CAS Center for Excellence in Particle Physics (CCEPP). The work of B.R. Littlejohn and P. T. Surukuchi was supported by the DOE Office of Science, under award No. DE-SC0008347.