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\sigma$ with respect to a suppression of the $^{235}\text{U}$ reactor antineutrino flux and at $2.5\sigma$ with respect to variations of the $^{235}\text{U}$ and $^{239}\text{Pu}$ 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 $^{235}\text{U}$ flux at $3.1\sigma$ and variations of the $^{235}\text{U}$ and $^{239}\text{Pu}$ fluxes at $2.8\sigma$. 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 $^{235}\text{U}$ or $^{239}\text{Pu}$ 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 $^{235}\text{U}$ and $^{239}\text{Pu}$ 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 $^{239}\text{Pu}$ flux and relatively large oscillations.


I. INTRODUCTION
The Daya Bay collaboration presented recently [1] the results of the measurement of the correlation between the reactor fuel evolution and the changes in the antineutrino detection rate which is quantified by the cross section per fission σ f , given by where F a 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 Fig. 2 of Ref. [1] 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 [2][3][4] 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 [1] 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 arXiv:1708.01133v1 [hep-ph] 3 Aug 2017 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 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 [2][3][4].
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 leastsquares 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. [5]).
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 leastsquares 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 leastsquares 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 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. [6]). 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 leastsquares 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 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 ).
In Section II we analyze the Daya Bay evolution data, in Section III we analyze the reactor antineutrino data which were available before the release of the Daya Bay fuel evolution data in Ref. [1], and in Section IV we perform the combined analysis.

II. DAYA BAY EVOLUTION
The results of the different fits of the Daya Bay evolution data are given in Tab of the minimum χ 2 , the number of degrees of freedom and the goodness-of-fit. In Tab. I 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 [2][3][4]. 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).

F2:
where The inequality (7) is satisfied for with From Tab. I one can see that all the fits have acceptable goodness-of-fit, but the OSC fit corresponding to activesterile 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 [1] 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 Fig. 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 [1].
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 (9) becomes This condition is satisfied, because numerically the lefthand 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 (9) 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. (3) can be obtained from that in Eq. (6) with the constraint P ee = 1. In this comparison, we have ∆χ 2 = 0.2 with one degree of freedom. Hence, the null hypothesis 235 cannot be rejected in favor of the alternative more complex hypothesis 235+OSC.
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 [1] without considering the theoretical uncertainties.
which determines the 235 U cross section per fission to be σ f,235 = 6.20 ± 0.15.
In the case of the 235+239 model, Fig. 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 [1]. We obtained which imply σ f,235 = 6.17 ± 0.16,  Fig. 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,

The vertical lines in
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 Fig. 4 obtained with the OSC. In Fig. 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 Tab. I.
On the other hand, in Fig. 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 Fig. 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.
The results of the fits with the models described in Section I are listed in Tab. II and the fit of the data is    From Tab. II one can see that all the model have an excellent goodness-of-fit, but the models OSC, 235+OSC, and 239+OSC with active-sterile neutrino oscillations have a significantly lower value of χ 2 min . This is due to the different source-detector distances in the experiments. As one can see from Fig. 7, where the reactor data are ordered by increasing values of the source-detector distance L. One can see that active-sterile oscillations can fit better the data of the short-baseline experiments which have a source detector distance between about 10 and 100 m. On the other hand, the poor fit of the data with the 235 model is explained by the lack of a trend Fig. 8 Figure 2 shows the marginal ∆χ 2 = χ 2 − χ 2 min for the factor r 235 obtained from the fit of the reactor rates in the 235 model. The result is which gives σ f,235 = 6.28 ± 0.08.
This is a determination of σ f,235 with smaller uncertainty than that obtained in Eq. (13) 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 which imply σ f,235 = 6.36 ± 0.09, (23) σ f,239 = 3.84 ± 0.28.
(24) 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 Tab. II, 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 Figs. 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 Tab. II.

IV. COMBINED ANALYSIS
In this section we present the results of the combined fits of the reactor rates in Table 1 of Ref. [5] (without the 2016 Daya Bay rate) and the 2017 Daya Bay evolution data [1].
The results of the fits with the models described in Section I are listed in Tab. III.
From Tab. III one can see that all the models have an excellent goodness-of-fit. The OSC model has a better goodness-of-fit than the 235 model. There is little improvement of the goodness-of-fit from the 235 model to the 235+239 model, whereas the goodness-of-fit improves significantly in the 235+OSC model and especially in the 239+OSC model.
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 activesterile 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 235 and 235+239 models have ∆χ 2 min = 2.3 and 1.8 with respect to the OSC model and our Monte Carlo comparison disfavors them at 1.7σ and 2.2σ, respectively.
From Fig. 2 one can see that in the 235 model the combined fit indicates a value of r 235 intermediate between those obtained from the analyzes of the Daya Bay evolution data and the reactor rates. The result is which gives σ f,235 = 6.25 ± 0.07.
This is a determination of σ f,235 with smaller uncertainty than that obtained in Eq. (13) from the Daya Bay evolution data and that obtained in Eq. (20) 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 r 235 = 0.934 ± 0.009, (27) which give  Table 1 of Ref. [5] (without the 2016 Daya Bay rate) and the 2017 Daya Bay evolution data [1].
Within the uncertainties, these results are compatible with those obtained in Ref. [24] 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. [1] whereas in Ref. [24] the Daya Bay evolution data have been taken into account with a Gaussian approximation of the χ 2 distribution in Fig. 3 of Ref. [1]. 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 Figs. 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 Tab. III. 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 Tab. III.

V. CONCLUSIONS
In this paper we analyzed the Daya Bay evolution data [1] in the 235, 235+239, OSC, 235+OSC, and 239+OSC models described in Section I, 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 [2][3][4] and short-baseline active-sterile 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 [2][3][4]. 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 [2][3][4]. 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 [1] 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. [1]. In this case, we found that the best explanation of the data is the OSC model with active-sterile neutrino oscillations, which depend on the source-detector distance and fit the rates measured by reactor experiments with a sourcedetector 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 [2][3][4] and relatively large active-sterile neutrino oscillations.
In conclusion, although the recent Daya Bay evolution data [1] 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 [3] still remains open. We hope that it will be solved soon by the new shortbaseline reactor neutrino experiments which will measure the reactor antineutrino flux from reactors with different fuel compositions: highly enriched 235 U research reactors for PROSPECT [25], SoLid [26], and STEREO [27], and commercial reactors with mixed fuel compositions for DANSS [28] and Neutrino-4 [29].