Analysis of the reactor antineutrino spectrum anomaly with fuel burnup

Recently, three successful antineutrino experiments (Daya Bay, Double Chooz, and RENO) measured the neutrino mixing angle θ13; however, significant discrepancies were found, both in the absolute flux and spectral shape. Much effort has been expended investigating the possible reasons for the discrepancies. In this paper, the change of neutrino energy spectrum with burnup is analyzed from the point of view of the change of neutrino energy spectrum with burnup. An accurate method for calculating neutrino energy spectrum is proposed. The non-equilibrium correction is studied by using this method. It is found that the non-equilibrium correction contributes not only to the energy region less than 4.0 MeV, but also to the energy region greater than 4.0 MeV, with a maximum correction of about 3%.


Introduction
The antineutrino experiment, such as Daya Bay [1], Double Chooz [2] and RENO [3], have measured the neutrino mixing angle θ 13 . However, significant discrepancies were found both in the absolute flux [4] and spectral shape [5][6][7]. For the antineutrino spectral shape, a 2.9σ deviation was found in the measured inverse beta decay position energy spectrum compared to predictions. In particular, an excess of events at energies of 4 -6 MeV was found in the measured spectrum [5][6][7], with a local significance of 4.4σ. These results have brought home the notion that neutrino fluxes are not as well understood as had been thought. At present, it is not clear what physical processes give rise to the neutrino spectra bump. In general, the time dependence of the antineutrino flux and spectrum ofν e from a reactor can be estimated using where W th (t) is the reactor thermal power, f i (t) is the fission fraction associated with each isotope, e i is the thermal energy release per fission for each isotope, S i (E) is a function of theν e energy E signifying theν e yield per fission for each isotope, c ne i (E, t) is the non-equilibrium correction of the long-live fission fragment isotopes, and S S NF (E, t) is the yield from the spent nuclear fuel (SNF) [8]. In this flux and spectrum prediction model, the spectrum of each isotopes are assumed as constant in each energy bin, and only the fission fraction of each isotopes vary with burnup. If the half life time of each fission product of isotopes are very small and its neutron absorbtion cross section are relative large, the constant spectrum assumption is reasonable, on the the other hands, it is no reasonable. The antineutrino spectrum of each isotopes are estimated * e-mail: maxb@ncepu.edu.cn from ILL electron spectra which were measured in the reactor, and the measured isotopes are only irradiated 12h for U and 36h for Pu. For neutrino reactor experiments, the irradiation time scale would rather be reactor cycle duration, typically 1 year or 1.5 years. Among the fission products, about 10% of the have a β-decay lifetime long enough to keep accumulating after several days. Moreover, the neutron energy spectrum of a standard Pressurized Water Reactor (PWR) have more important epithermal and fast neutron energy components than in the ILL measurements. Because of these reasons, the antineutrino spectrum of each isotopes has to be corrected, this is called as 'Nonequilibrium corrections'. Nonequilibrium corrections were also studied in the reference, and its corrections with fuel burnp are not discussed in detail. The purpose of this study is to verify whether the constant spectrum prediction model is accurate and what is the relation between this model and the antineutrino spectrum anomaly. In this study, a more precise spectrum calculation method are used based on the abinital method. The nonequilibrium corrections vary with fuel burnup are given and its possible of having contribution to the Bump are also discussed.
Antineutrino spectrum calculation method was introduced in Section 2, and a precise antineutrino spectrum calculation method was proposed. Reactor simulation was discussed in Section 3, and total antineutrino spectra of the reactor at different burnup with the precise method was given. The spectra comparison at different burnup was discussed in the Section 4. The last section is the conclusion.

Antineutrino spectrum calculation
There are four important isotopes in the antineutrino spectrum prediction, such as 235 U, 238 U, 239 Pu, 241 Pu. Total fission fraction of these four isotopes are more than 99%. Therefor, when we predict the antineutrino spectrum of a reactor, only these four isotopes are taking into account. The antineutrino spectrum of a reactor usually calculated by equation 2.
the four β-spectra of the main fuel isotopes are measured directly without the knowledge of the underlying branch-level processes.
• method-2 : The nuclear database and Fermi theory are used to calculate the antineutrino spectra of each fission product, then summate each product spectra using the cumulative fission yield as weight. This method is always called ab initial method.
• method-3 : Equation 3 and 4 are used to calculate the total antineutrino spectrum.
In the summation method, the decay rate of each fission product are needed. For the short-lived isotopes, and its neutron absorption cross section, the decay rate can be approximately equal total fission rate of parent isotope times the cumulative yield of the isotopes i. In this study, we proposed the method-3, which can avoid this approximation because the decay rate has bee calculated by equation 4 accurately.
The advantage of method-3 is that is taken into account the nonequilibrium corrections, and the difference between method-3 and method-1 can be explain the nonequilibrium correction. This will be discussed in the later.

Reactor simulation and spectra
To calculate the antineutrino spectra by applying equation 3, the atomic density must be obtained. Reactor simulation code is necessary because of a strong energy dependence in the neutron cross section, a severe local angular dependence, and a temperature dependence in the neutron cross section. The Daya Bay reactor operates with 157 fuel assemblies producing a total thermal power of 2895 MW. The assembly is a 17 × 17 desigin, for 289 rods, which contains 264 fuel rods, 24 control rods and one control rod guide tube. The enrichment of the new fuel of the Daya Bay reactor core is 4.45% for the 18month reload core design [9]. A typical PWR pin cell with a fuel enrichment 4.45% was made and Reactor Monte Carlo code (RMC) [10] was used to solving the transport equation. A database containing about 1946 nuclide (including 1119 fission product) is used in the depletion calculation, and the Evaluated Nuclear Data File ENDF/B-VII.0 [11]are used in the transport calculation. Fission fraction of the main four isotopes as a function of burnup is shown in figure 1. As shown in the figure 1, the fission fraction of 235 U and 238 U decrease with fuel burunp increasing because of its atomic density in the core decrease with fuel burnup increasing. For the 239 Pu and 241 Pu, the fission fraction increase with fuel burnup because of its atomic density increasing.
According to the reaction simulation, the atomic density of each fission product is obtained from the simulation results. Beside the atomic density, the decay constant can be found in the nuclear database [16].
The beta-decay spectrum S n (E ν ) for a single transition in the nucleus (Z,A) with end-point energy E 0 is 3 , E e (p e ) is the electron total energy (momentum), F(E e , Z, A) the Fermi function needed to account for the Coulomb interaction of the outgoing electron with the charge of the daughter nucleus, and C(E e ) a shape factor for the forbidden transitions arising from additional electron momentum terms. The term δ(E e , Z, A) represents a fractional correction to the spectrum. The primary corrections to the beta-decay analysis are radiative δ rad , finite size δ FS , and weak magnetism δ W M . For the allowed transitions C(E e ) = 1, the radiative, finite size, and weak magnetism corrections are taken from [12]. For the forbidden transitions, the shape factor and fractional weak magnetism corrections are taken from [13]. In the present work, the cumulative fission yields are taken from the ENDF/B-VII.1 library[16]. The decay data are taken from a version of JENDL-3.1 [14].

Results and analysis
After the atomic density, decay constant and the antineutrino spectrum of each fission product, the total antineutrino spectra at different burnup can be calculated by applying equation 3 and equation 4 as shown in figure 2. As shown in figure 2, the spectra difference of different fuel   burnup is not much. In order to investigate the fuel burunp to the antineutrino spectra, we define the spectra difference in equation 6 The spectra difference is the spectra which the spectra of fuel burnup t minus that of beginning of cycle. The spectra difference is shown in figure 3. As shown in figure 3, the spectra difference increases with increasing fuel burnup.     In ILL (Grenoble, France) experiment, the β spectrum of 235 U, 239 Pu, and 241 Pu were measured, and then the antineutrino spectrum were inferred from a beta spectrum with the assumption of 25 virtual beta spectrum [17][18][19].
In order to dealing with the reactor anomaly, the reactor antineutrino flux were examined based on the ILL online measurements of the integral beta spectrum of the fission products. The improvements on the earlier analysis of ILL integral measurements caused to an increased energy of antineutrino flux [20,21], which was verified in Huber's analysis [12]. The spectra difference using equation 2 with Huber-Mueller model ( 235 U, 239 Pu, 241 Pu from Huber model, 238 U from Mueller model) is shown in figure  4. In the method-3 calculation, as we discussed before, the nonequilibrium corrections has been included in the spectrum calculation. However, it has not been included in the Huber-Mueller model because of the irradiated time for U and Pu is relative short comparison with reactor core life. If we define ∆S (E ν , t) should be represent the nonequilibrium corrections, and it is shown in figure 5. As shown in figure 5, the nonequilibrium corrections appear not only in the low energy (below 4.0 MeV ), but also in the high energy region (above 4.0MeV ). This is not consistent with previous study [8,21,22] because of different assumption. In previous study, we always regards as only long-lived fission products which having the property of beta end energy point larger than 1.8 MeV and T 1/2 ≥ 10 h are chosen. The end energy point of beta decay of these isotopes, such as 90 Y, 106 RH, 144 Pr, and so on, are less than 4.0 MeV, therefor, the nonequilibrium corrections appeared in the below 4.0 MeV energy region. If we only take these isotopes contribution, the nonequilibrium corrections are shown in figure 6. As shown in figure 5 and figure 6, the corrections of below 4.0 MeV has similar trend because of the same isotopes contribution. The corrections above 4.0 MeV may caused by the isotopes of E d ≥ 1.8 MeV and T 1/2 < 10 h. In previous study, we always assume that those isotopes have been reach the equilibrium and its contribution to the antineutrino spectrum have been included in the ILL spectrum. If this is true, method-2 will be the same as method-3. It can be seen from our results that there are still some isotopes of E d ≥ 1.8 MeV and T 1/2 < 10 h which have no been equilibrium, and its contribution to the nonequilibrium corrections have about maximum −3%. Antineutrino energy (MeV) The difference between two method(%)  Comparing with the measured antineutrino spectrum in Daya Bay, it is found that the neutrino energy spectrum in the middle of cycle agrees better with the measured energy spectrum in Daya Bay. The energy spectrum at the beginning of cycle is lower than that at 4.0 MeV. On the contrary, the neutrino energy spectrum at the end of cycle is higher than that at 4.0 MeV. This is mainly due to the different fission fraction of 239 Pu and 235 U. Compared with Huber-Mueller model, the neutrino energy spectrum in the lifetime agrees better with the experimental data.

Conclusion
In order to consider the relationship between neutrino energy spectrum and burnup, a new method based on ab initio calculation is proposed, which can consider the contributions of all equilibrium and non-equilibrium nuclides. By comparing the results of this method with those of Huber-Mueller model, the influence of non-equilibrium state correction is analyzed. It is found that the nonequilibrium correction contributes not only to the energy region less than 4.0 MeV, but also to the energy region greater than 4.0 MeV, with a maximum correction of about 3%. At the same time, the neutrinos calculated by this method are compared with the Daya Bay experiment results.