Reactor fuel fraction information on the antineutrino anomaly

We analyzed the evolution data of the Daya Bay reactor neutrino experiment in terms of short-baseline active-sterile neutrino oscillations taking into account the theoretical uncertainties of the reactor antineutrino fluxes. We found that oscillations are disfavored at 2.6σ with respect to a suppression of the 235U reactor antineutrino flux and at 2.5σ with respect to variations of the 235U and 239Pu fluxes. On the other hand, the analysis of the rates of the short-baseline reactor neutrino experiments favor active-sterile neutrino oscillations and disfavor the suppression of the 235U flux at 3.1σ and variations of the 235U and 239Pu fluxes at 2.8σ. We also found that both the Daya Bay evolution data and the global rate data are well-fitted with composite hypotheses including variations of the 235U or 239Pu fluxes in addition to active-sterile neutrino oscillations. A combined analysis of the Daya Bay evolution data and the global rate data shows a slight preference for oscillations with respect to variations of the 235U and 239Pu fluxes. However, the best fits of the combined data are given by the composite models, with a preference for the model with an enhancement of the 239Pu flux and relatively large oscillations.


Introduction
Nuclear reactors are copious sources of electron antineutrinos (ν e ), which have been detected by many experiments starting from the historical Reines and Cowan experiment in 1953-56 (see the reviews in refs. [1][2][3]). The electron antineutrinos are produced by β decays of the neutron-rich nuclei generated by the fissions of the four fissionable isotopes 235 U, 238 U, 239 Pu, and 241 Pu. The totalν e flux emitted by each reactor is the sum of the fluxes generated by the four fissionable isotopes weighted by the effective fuel fractions which are monitored in time by the reactor managers.
The time evolution of nuclear reactors is divided in cycles of length that can go from about a month to one or two years. At the beginning of each cycle all or part of the fuel is replaced with fresh fuel, which is typically composed by uranium enriched with the fissile isotope 235 U with respect to the natural abundance, which is about 99.3% of 238 U and 0.7% of 235 U. At the start of each reactor cycle the main contribution to theν e flux comes from the fissions of 235 U, with a small contribution of the fissionable isotope 238 U. The neutron flux produced by the fissions generate mainly 239 Pu and in smaller quantity 241 Pu. Hence, as 235 U is consumed during a cycle, the 235 U contribution to theν e flux decreases and the contributions 239 Pu and 241 Pu increase, with a dominance of the 239 Pu contribution, which can be comparable with the 235 U contribution towards the end of each cycle.
For our discussion it is useful to distinguish two types of nuclear reactors: research reactors and power reactors. In research reactors the fresh fuel is made of almost pure 235 U and the cycles are short and giveν e fluxes which are practically due only to 235 U.
In most commercial power reactors the fresh fuel is made of uranium enriched by some percent of 235 U and the cycles are typically between one and two years. Therefore, the time evolution of the fuel composition must be taken into account in the calculation of thē ν e flux of power reactors.
Theν e fluxes generated by the four fissionable isotopes have been calculated most recently in 2011 [4,5], where they were found to be a few percent larger than in previous JHEP10(2017)143 calculations [6][7][8]. The resulting predictions for the rates ofν e -induced events measured in several neutrino experiments are smaller than the experimental rates. The average deficit of about 5% is the so-called "reactor antineutrino anomaly" [9], which can be due to two causes: a) some mistake in theν e flux calculations and/or b) the disappearance ofν e during their propagation from the reactor to the detector. This disappearance is most likely due to active-sterile neutrino oscillations (see the review in ref. [10]).
The Daya Bay collaboration presented recently [11] the results of the measurement of the correlation between the reactor fuel evolution and the changes in the antineutrino detection rate in the Daya Bay experiment, which detectsν e 's produced by two complexes of power reactors at distances of about 360 m and 500 m. Theν e detection rate is quantified by the cross section per fission σ f , given by where F i and σ f,i are the effective fission fractions and the cross sections per fission of the four fissionable isotopes 235 U, 238 U, 239 Pu, 241 Pu, denoted, respectively, with the label i = 235, 238, 239, 241. The Daya Bay collaboration presented in figure 2 of ref. [11] the values of σ f for eight values of the effective 239 Pu fission fraction F 239 . They fitted these data allowing variations of the two main cross sections per fission σ f,235 and σ f,239 , with the assumption that σ f,238 and σ f,241 have the Saclay+Huber theoretical values [4,5,9] with enlarged 10% uncertainties. They also compared the best-fit of this analysis with the best-fit obtained under the hypothesis of active-sterile neutrino oscillations, which predicts the same suppression for the four cross sections per fission with respect to their theoretical value. They obtained ∆χ 2 /NDF = 7.9/1, corresponding to a p-value of 0.49%, which disfavors the active-sterile oscillations hypothesis by 2.8σ. In this calculation the uncertainties of the theoretical calculation of the four cross sections per fission were not taken into account.
In this paper we present the results of analyses of the Daya Bay evolution data [11] with least-squares functions that take into account explicitly the uncertainties of the theoretical calculation of the four cross sections per fission. Moreover, we consider additional models with independent variations of the 235 U and 239 Pu fluxes with and without active-sterile neutrino oscillations, and we extend the analysis taking into account also the information on the cross sections per fission of all the other reactor antineutrino experiments which have different fuel fractions. We also perform proper statistical comparisons of the non-nested models under consideration through Monte Carlo estimations of the p-values.
Given a set of data labeled with the index a on the cross section per fission for different values of the fuel fractions, we write the theoretical predictions as where i = 235, 238, 239, 241 and σ SH f,i are the Saclay+Huber cross sections per fission. The coefficients r i are introduced in order to take into account the uncertainties of the

JHEP10(2017)143
Saclay+Huber cross sections per fission or to study independent variations of the antineutrino fluxes from the four fissionable isotopes with respect to the Saclay+Huber theoretical values [4,5,9].
We consider the following models: 235 A variation of the cross section per fission of the antineutrino flux from 235 U only.
In this case, we analyze the data with the least-squares statistic where σ exp f,a are the measured cross sections per fission, V exp is the experimental covariance matrix, and V SH is the covariance matrix of the fractional uncertainties of the Saclay-Huber theoretical calculation of the antineutrino fluxes from the four fissionable isotopes (given in table 3 of ref. [12]).
In this analysis there is only one parameter determined by the fit: r 235 . The parameters r 238 , r 239 , and r 241 are nuisance parameters.
235+239 Independent variations of the cross sections per fission of the antineutrino fluxes from 235 U and 239 Pu.
In this case, we analyze the data with the least-squares statistic In this analysis there are two parameters determined by the fit: r 235 and r 239 . The parameters r 238 and r 241 are nuisance parameters.
OSC Active-sterile neutrino oscillations, in which the measured cross sections per fission are suppressed with respect to the theoretical cross sections per fission σ th f,a by the survival probability P ee which is independent of the 239 Pu fraction F 239 .
In this case, we analyze the data with the least-squares statistic In the analysis of the Daya Bay evolution data there is only one parameter determined by the fit: P ee . The parameters r 235 , r 238 , r 239 , and r 241 are nuisance parameters.
In the analysis of the other reactor antineutrino data we take into account that P ee

JHEP10(2017)143
depends on the neutrino mixing parameters ∆m 2 41 and sin 2 2ϑ ee in the simplest 3+1 active-sterile neutrino mixing model (see ref. [10]). Hence, in this case there are two parameters determined by the fit: ∆m 2 41 and sin 2 2ϑ ee .
235+OSC A variation of the cross section per fission of the antineutrino flux from 235 U and active-sterile neutrino oscillations with a survival probability P ee as in the OSC model.
In this case, we analyze the data with the least-squares statistic In the analysis of the Daya Bay evolution data there are two parameters determined by the fit: r 235 and P ee . The parameters r 238 , r 239 and r 241 are nuisance parameters.
In the analysis of the other reactor antineutrino data we take into account that P ee depends on ∆m 2 41 and sin 2 2ϑ ee as in the OSC model. Therefore, in this case there are three parameters determined by the fit: r 235 , ∆m 2 41 , and sin 2 2ϑ ee .
239+OSC This model is similar to the 235+OSC model, with 235 U 239 Pu. The number of parameters determined by the fit is two in the analysis of the Daya Bay evolution data (r 239 and P ee ) and three in the analysis of the other reactor antineutrino data (r 239 , ∆m 2 41 , and sin 2 2ϑ ee ).
Note that: 1. In all the models we constrained the values of σ f,238 and σ f,241 around the corresponding theoretical values with the theoretical uncertainties. In principle it would be interesting to consider also σ f,238 and σ f,241 as free parameters to be determined by the fit of the data. However, this is not possible in practice because the data do not constrain them in a sufficient way. Indeed, also the Daya Bay collaboration [11] constrained σ f,238 and σ f,241 around the theoretical values, albeit with an arbitrary 10% uncertainty. We adopt the theoretical uncertainties because they are well defined and motivated, and it is possible that the corresponding predictions are correct, whereas those of σ f,235 and σ f,239 are not.
2. We did not consider the case 235+239+OSC, because this model is not constrained by the data. This is due to the fact that arbitrarily large values of σ f,235 and σ f,239 can be counterbalanced by an arbitrary small P ee .
The plan of the paper is as follows: in section 2 we analyze the Daya Bay evolution data, in section 3 we analyze the reactor antineutrino data which were available before the release of the Daya Bay fuel evolution data in ref. [11], in section 4 we perform the combined analysis, and in section 5 we draw our conclusions. Table 1. Fits of the Daya Bay evolution data [11].  Figure 1. Fits of the Daya Bay evolution data [11] normalized to the Saclay-Huber theoretical predictions [4,5,9]. The error bars show only the uncorrelated statistical uncertainties.

Daya Bay evolution
The results of the different fits of the Daya Bay evolution data are given in table 1 where we list the values of the minimum χ 2 , the number of degrees of freedom and the goodness-of-fit. In table 1 we also list the best-fit values of the fitted parameters. Figure 1 shows the comparison of the different fits with the Daya Bay evolution data normalized to the Saclay-Huber theoretical cross sections per fission [4,5,9]. Note that the Daya Bay evolution data have the following two important features: F1 A suppression of σ f with respect to σ SH f in agreement with the reactor antineutrino anomaly. This feature can be fitted with at least one of the r i and P ee smaller than one (if the others are equal to one).

JHEP10(2017)143
From table 1 one can see that all the fits have acceptable goodness-of-fit, but the OSC fit corresponding to active-sterile oscillations has a goodness-of-fit which is significantly lower than the others, because it corresponds to a constant σ f /σ SH f and cannot fit feature F2.
The results of our analysis agree with the conclusion of the Daya Bay collaboration [11] that the 235 model fits well the data and little is gained by allowing also the variation of σ f,239 in the 235+239 model. The shift in figure 1 of the line corresponding to the 235 model with respect to an ideal line fitting the data by eye is allowed by the large correlated systematic uncertainties of the Daya Bay bins [11].
The excellent fit in the 235 model is due to the fact that it can fit the two features of the Daya Bay evolution data listed above. It can obviously fit feature F1 with r 235 < 1. It can also fit feature F2, because for r 235 < 1 and r 238 = r 239 = r 241 = 1 the condition ( This condition is satisfied, because numerically the left-hand side is about 1.30 and the right-hand side is between 0.20 and 0.24.
Obviously, the 235+OSC model can provide a fit which is at least as good as the 235 model, with the additional possibility to improve the fit of feature F1 with P ee < 1.
It is maybe more surprising that also the 239+OSC model fits better than the 235 model for r 239 > 1. This can happen because the condition (2.3) for fitting feature F2 is always satisfied for r 239 > 1 and r 235 = r 238 = r 241 = 1. Then, a sufficiently small value of P ee < 1 allows us to fit feature F1 in spite of the increase of σ th f,a due to r 239 > 1. Nested models can be compared in the frequentist approach by calculating the p-value of the χ 2 min difference, which has a χ 2 distribution corresponding to the difference of the number of degrees of freedom of the two models. With this method we can compare only the nested models 235 and 235+OSC, because the χ 2 in eq. (1.3) can be obtained from that in eq. (1.6) with the constraint P ee = 1. In this comparison, we have ∆χ 2 = 0.  Also non-nested models can be compared considering the χ 2 min difference, but one must calculate the p-value with a Monte Carlo. In this case one must consider as the null hypothesis the model which has the higher χ 2 min and generate many sets of synthetic data assuming the null hypothesis. The fits of all the sets of synthetic data with the two models under consideration gives the distribution of the χ 2 min difference from which one can calculate the p-value of the observed χ 2 min difference. We do not bother to consider the comparison of the 235 and 235+239 models, since the small ∆χ 2 min = 0.2 cannot lead to the rejection of the null hypothesis 235. On the other hand, it is interesting to compare the OSC and 235 models which have ∆χ 2 min = 5.7. According to our Monte Carlo simulation, the p-value of the null hypothesis OSC is 0.85%. Hence, the comparison of the OSC and 235 models disfavors the OSC model at the 2.6σ level.
We also compared with a Monte Carlo the OSC and 235+239 models which have ∆χ 2 min = 5.9. We found that the null hypothesis OSC has a p-value of 1.3%, which is larger than in the previous case because the 235+239 model has one parameter more than the 235 model. Thus, in this case, the OSC model is disfavored at the 2.5σ level, which is slightly less stringent than the 2.8σ obtained by the Daya Bay collaboration [11] without considering the theoretical uncertainties.  In the case of the 235+239 model, figure 3 show that the Daya Bay evolution data indicate a larger suppression of σ f,235 than σ f,239 , in agreement with the results of the analysis of the Daya Bay collaboration [11]. We obtained  (10% in the calculation of the Daya Bay collaboration and the Saclay+Huber theoretical values [4,5,9] 8.15% and 2.60% in our calculation).
The vertical lines in figure 4 show the bounds on sin 2 2ϑ ee = 2(1 − P ee ) obtained in the OSC analysis of the Daya Bay evolution data, in which oscillations are averaged because of the large source-detector distance. One can see that sin 2 2ϑ ee = 0.12 ± 0.06, (2.12) and there is no lower bound at 2σ, because oscillations are favored over the no-oscillation case only at the 1.9σ level. Figure 5 shows that the variation of r 235 in the 235+OSC model causes a shift of the allowed region for sin 2 2ϑ ee towards lower values with respect to figure 4 obtained with the OSC. In figure 5 there is no lower bound at 1σ, because oscillations are favored over the no-oscillation case only at 0.4σ. This is due to the preference for values of r 235 smaller than one, as shown by the best-fit value in table 1.
On the other hand, in figure 6 corresponding to the 239+OSC there is a shift of the allowed region for sin 2 2ϑ ee towards larger values with respect to figure 4 obtained with the OSC, because P ee is smaller in order to compensate the increase of σ th f,a due to r 239 > 1.  Figure 5. Allowed regions in the sin 2 2ϑ ee -∆m 2 41 plane obtained from the fit of the Daya Bay evolution data [11] (Daya Bay), from the fit of the reactor rates (Rates), and from the combined fit (Combined) with the 235+OSC model. The best fit points are indicated by crosses, except for the fit of the Daya Bay evolution data for which the best fit is the vertical dash-dotted line. For the Daya Bay and Rates fits the 1σ, 2σ, and 3σ allowed regions are limited, respectively, by solid, dashed, and dotted lines.
The results of the fits with the models described in section 1 are listed in table 2       which gives σ f,235 = 6.28 ± 0.08.
This is a determination of σ f,235 with smaller uncertainty than that obtained in eq. (2.7) from the Daya Bay evolution data. Figure 3 shows that in the case of the 235+239 model the determination of r 235 and r 239 is quite different in the analyses of the Daya Bay evolution data and the reactor rates. In the first analysis r 235 and r 239 are correlated, whereas in the second analysis they are slightly anticorrelated. Moreover, the analysis of the reactor rates prefers a larger value of r 235 and a smaller value of r 239 than the analysis of the Daya Bay evolution data. The results of the analysis of the reactor rates are  Figure 4 show the allowed region in the sin 2 2ϑ ee -∆m 2 41 plane obtained from the fit of the reactor rates in the OSC model. One can see that there in only one region allowed at 1σ around the best-fit point given in table 2, but the 2σ allowed regions do not have an upper bound for ∆m 2 41 . The 3σ allowed region does not have a lower bound for sin 2 2ϑ ee , because oscillations are favored over the no-oscillation case only at the 2.7σ level.
From a comparison of figures 4, 5, and 6 one can see that the variations of r 235 and r 239 in the 235+OSC and 239+OSC models, respectively, have small effects on the allowed region in the sin 2 2ϑ ee -∆m 2 41 plane, in agreement with the best-fit values close to one of r 235 and r 239 in table 2.

Combined analysis
In this section we present the results of the combined fits of the reactor rates in table 1 of ref. [12] (without the 2016 Daya Bay rate) and the 2017 Daya Bay evolution data [11].
The results of the fits with the models described in section 1 are listed in  Table 3. Fits of the reactor rates in table 1 of ref. [12] (without the 2016 Daya Bay rate) and the 2017 Daya Bay evolution data [11].
The comparison of the nested models 235 and 235+OSC give ∆χ 2 = 5.1 with two degrees of freedom. Hence, the p-value of the null hypothesis 235 is 7.8% and it can be rejected in favor of the introduction of active-sterile neutrino oscillations only at 1.8σ. As a check, with a Monte Carlo simulation we obtained a p-value of 5.1%, which corresponds to 1.9σ.
The  which gives σ f,235 = 6.25 ± 0.07. (4.2) This is a determination of σ f,235 with smaller uncertainty than that obtained in eq. (2.7) from the Daya Bay evolution data and that obtained in eq. (3.2) from the reactor rates. Figure 3 shows that in the case of the 235+239 model the determination of r 235 and r 239 from the combined fit improves the uncertainties of the two parameters with respect to those obtained from the separate analyses of the Daya Bay evolution data and the reactor rates The results are

JHEP10(2017)143
Within the uncertainties, these results are compatible with those obtained in ref. [30] with different assumptions on the uncertainties of σ f,238 and σ f,241 . Note that here we performed a full analysis of the Daya Bay evolution data using the complete information available in the Supplemental Material of ref. [11] whereas in ref. [30] the Daya Bay evolution data have been taken into account with a Gaussian approximation of the χ 2 distribution in figure 3 of ref. [11]. Figure 4 show the allowed region in the sin 2 2ϑ ee -∆m 2 41 plane in the OSC model. The allowed regions are smaller than those obtained from the fit of the reactor rates and there is a 3σ lower bound for sin 2 2ϑ ee , because oscillations are favored over the no-oscillation case at the 3.1σ level. However, there is no upper bound for ∆m 2 41 at 2σ, because at that confidence level the data can be fitted with an averaged oscillation probability which does not depend on the source-detector distance.
Comparing figures 4 and 5, one can see that the variation of r 235 in the 235+OSC enlarges the allowed regions towards lower values of sin 2 2ϑ ee and there is no lower bound for sin 2 2ϑ ee at 2σ, because oscillations are favored over the no-oscillation case only at 1.4σ. This is due to the preference for values of r 235 smaller than one, as shown by the best-fit value in table 3. Figure 6 shows that the best-fitting model 239+OSC gives the strongest indication in favor of oscillations, which are favored over the no-oscillation case at 3.0σ. This is due to the preference for values of r 239 larger than one, as shown by the best-fit value in table 3.

Conclusions
In this paper we analyzed the Daya Bay evolution data [11] in the 235, 235+239, OSC, 235+OSC, and 239+OSC models described in section 1, which allow to compare the fits of the data under the hypotheses of variations of the 235 U and 239 Pu reactor antineutrino fluxes with respect to the Saclay+Huber theoretical value [4,5,9] and short-baseline activesterile neutrino oscillations, taking into account the theoretical uncertainties of the reactor antineutrino fluxes. We found that the best explanation of the Daya Bay evolution data is the 235 model with a variation of the 235 U flux with respect to the Saclay+Huber theoretical value [4,5,9]. Comparing the OSC model of active-sterile neutrino oscillations with the 235 model, we found that it is disfavored at 2.6σ.
We also compared the OSC model with the 235+239 model which allows independent variations of the 235 U and 239 Pu fluxes with respect to Saclay+Huber theoretical values [4,5,9]. We found that the OSC model is disfavored at 2.5σ. This result is slightly less stringent than the 2.8σ obtained by the Daya Bay collaboration [11] without considering the theoretical uncertainties.
The Daya Bay evolution data can also be fitted well with the 235+OSC model, with a suppression of the 235 U flux and neutrino oscillations, or with the 239+OSC model, with an enhancement of the 239 Pu flux and relatively large neutrino oscillations.
We also performed a similar analysis of the reactor antineutrino data which were available before the release of the Daya Bay fuel evolution data in ref. [11]. In this case, we found that the best explanation of the data is the OSC model with active-sterile neutrino JHEP10(2017)143 oscillations, which depend on the source-detector distance and fit the rates measured by reactor experiments with a source-detector distance between about 10 and 100 m better than the distance-independent suppression of the reactor antineutrino flux given by suppressions of the 235 U and 239 Pu fluxes. In this case, the 235 model with a suppression of the 235 U flux only is disfavored at 3.1σ and the 235+239 model with independent suppressions of the 235 U and 239 Pu fluxes is disfavored at 2.8σ. As with the fit of the Daya Bay evolution data, composite models including both variations of the 235 U or 239 Pu fluxes and active-sterile oscillations provide good fits to the global reactor rate data.
Finally, we performed combined fits of the Daya Bay evolution data and the other reactor rates and we found that all the considered models fit well the data. The OSC model has a better goodness-of-fit than the 235 and 235+239 models, which are almost equivalent. We obtained better fits of the data with the composite 235+OSC and 239+OSC models. In particular, the best-fit model is 239+OSC, with an increase of the 239 Pu flux with respect to the Saclay+Huber theoretical value [4,5,9] and relatively large activesterile neutrino oscillations.
In conclusion, although the recent Daya Bay evolution data [11] disfavor short-baseline active-sterile neutrino oscillations over a suppression of the 235 U reactor antineutrino flux or independent suppressions of the 235 U and 239 Pu fluxes, the result is reversed in the analysis of the other available reactor antineutrino data. Both sets of data are individually well-fitted by composite models with variations of the 235 U or 239 Pu fluxes and activesterile neutrino oscillations. The combined data set indicates a preference for the composite models and, in particular, the best fit is obtained with the 239+OSC model, through an enhancement of the 239 Pu flux and relatively large oscillations. However, while these combined fits suggest a preference for models including sterile neutrinos, the significant uncertainties in the reactor rate measurements and the high goodness-of-fits observed for models both with and without sterile neutrinos make it clear that the search for the explanation of the reactor antineutrino anomaly [9] still remains open. We hope that it will be solved soon by the new short-baseline reactor neutrino experiments which will measure the reactor antineutrino flux from reactors with different fuel compositions: highly enriched 235 U research reactors for PROSPECT [31], SoLid [32], and STEREO [33], and commercial reactors with mixed fuel compositions for DANSS [34] and Neutrino-4 [35].