Ab initio calculations of reactor antineutrino fluxes with exact lepton wave functions

New \textit{ab initio} calculations of the isotopic reactor antineutrino fluxes are provided with exact numerical calculations of the lepton wave functions, assuming all the decay branches are allowed GT transitions. We illustrate that the analytical Fermi function and finite size effect each could have the largest spectral deviation of $\mathcal{O}(10\%)$, whereas the effect of their combination could result in spectral deviations at the level of 5\%---10\%. Meanwhile, we also find that several forms of the extended charge distributions have negligible effects on the spectral variation. Using the state-of-the-art nuclear databases, compared to usual \textit{ab initio} calculations using the analytical single beta decay spectrum, our new calculation shows sizable but opposite spectral deviations at the level of 2\%---4\% for the cumulative antineutrino and electron energy spectra which may partially contribute to the observed spectral excess in the high energy antineutrino range. Finally we observe that the effect of analytical beta decay spectrum approximation is rather universal for all the four fissionable isotopes.


Introduction
Electron antineutrinos from nuclear reactors have been widely used to study the fundamental properties of massive neutrinos [1]. Reactor antineutrinos are the beta-decaying products of fission fragments associated with four main fissionable isotopes 235 U, 238 U, 239 Pu, and 241 Pu. Predicting the reactor antineutrino flux and spectrum is always an important prerequisite for the reactor antineutrino experiments . There are two different methods to obtain the theoretical calculation of the reactor antineutrino flux and spectrum. The first one employs the ab initio method [3][4][5][6][7][8][9][10][11][12][13][14] by a direct summation of all the beta decay branches using the available information from the latest nuclear databases. The other method uses an effective conversion procedure of virtual branches [15][16][17][18][19][20][21][22][23] based on the measurements of the integral electron energy spectra of the fissionable isotopes.
During the reactor fuels burning, more than six thousand beta decay branches contribute to the antineutrinos produced in each fissionable isotope. In the first method of the ab initio calculation, the fission yield, the endpoints and branching ratios of each fission fragment can be accessed from the nuclear databases. The reactor antineutrino flux can be obtained by a direct summation of all the beta decay branches using an analytical description of the single beta decay spectrum [3][4][5][6][7][8][9][10][11][12][13][14]. Therefore, the calculating accuracy depends on the uncertainties of the fission yields and the beta decay information. The second effective conversion method has been developed using the corresponding electron spectrum associated with the thermal neutron induced fission of 235 U, 239 Pu, and 241 Pu at ILL, Grenoble, France in 1980s [15][16][17], and the fast neutron induced fission of 238 U at FRMII in Garching, Germany in 2011 [18]. The antineutrino flux is calculated by assuming dozens of virtual branches and fitting the electron spectrum of the fission isotope to be consistent with the measurement.
In the ab initio method, the fission yields have been evaluated by different nuclear databases, large uncertainties and even incompleteness for many important fragments are still a problem for the latest nuclear databases. In addition, the beta decay information including the branching ratios, endpoints and the quantum numbers of the initial and final states are not always known for the existing fragments. There is a well known problem of the Pandemonium effect [34] for some nuclear beta decay data that has been proved to play important roles in the reactor antineutrino flux calculations [8,14]. Moreover, although theory of the analytical description of single beta decay spectrum, based on the series expansions of the nuclear charge and radius, as well as the electron and antineutrino energies, has been developed since 1960s [35], the accuracy and validity of the analytical description on reactor antineutrino predictions, which have beta decay branches with large nuclear charge and/or high endpoint energies, have not been carefully studied until recently [36].
In this work, by assuming all the decay branches are allowed GT transitions we employ a new ab initio calculations of the isotopic reactor antineutrino fluxes using the exact calculations of the numerical lepton wave functions. We make a systematic test on the reliability and accuracy of the analytical description of the single beta decay spectrum, and observe that depending on the endpoint energies and nuclear charges, the analytical Fermi function and finite size effect could have the largest spectral deviation of O(10%), whereas the effect of their combination may result in spectral deviations at the level of 5%-10%. Using the stateof-the-art nuclear databases, compared to the usual ab initio calculations using the analytical single beta decay spectrum, our new calculation of reactor antineutrinos shows sizable but opposite spectral deviations at the level of 2%-4% and the effect of analytical beta decay spectrum approximation is rather universal for all the four fissionable isotopes.
This work is organized as follows. In Sec. 2 we review the description of beta decays theory and discuss how to make the approximation towards the usual analytical calculations. Then we investigate the dependence of these approximation on the beta decay endpoint energies and nuclear charges in Sec. 3. In Sec. 4 the new exact calculations of the single beta decay spectrum are applied to the four reactor fission isotopes, 235 U, 238 U, 239 Pu, and 241 Pu. Finally the concluding remarks are presented in Sec. 5.

Theoretical description of the beta decay
In this section, we want to review the theoretical description of the beta decay and discuss how to make the approximation towards the standard analytical calculations.
In the Standard Model, the nuclear beta decay is mediated by the charged SU (2) L gauge boson W whose mass is generated by the spontaneous symmetry breaking. For the nuclear decay with energy at the range of several MeV, the weak interaction can be effectively described by the contact interactions of hadron currents and lepton currents following conventions in [37]: where G β = G F cos θ C , with G F being the Fermi constant and θ C being the Cabbibo angle. ψ n , ψ p ψ e , and ψ ν are the wave functions of the neutron, proton, electron and antineutrino respectively. The bounded nucleus system usually have the spherical symmetry, thus in order to describe such system, the spherical coordinate is the most appropriate one. Therefore, the contact interaction can be decomposed to summations of various angular transfer interactions (We choose the definitions of spherical harmonics as in Appendix VIII of [38]): where Y M L is the spherical harmonic function,r i = r i /r i with r i = | r i | (i = 1, 2). Thus the Hamiltonian with the definite angular momentum and parity can be written as: where Ω N and Ω L are the solid angle coordinates of the nucleus and lepton systems respectively. T M KLs are the irreducible multi-pole tensor operators of rank K and in general can be expressed as the products of spherical harmonics and Dirac matrices, whose expressions can be found in Refs. [38,39]. K, L and s are the quantum numbers of the total, orbital and spin angular momenta of the mediating W boson, respectively, and M is the magnetic quantum number of the total angular momentum.
For a definite weak transition with known spin and parity of the initial and final nuclei, only those transition operators with certain selection rules are relevant. A naive estimation shows that the electron and antineutrino wave functions are proportional to (p e R) le and (p ν R) lν respectively when p e R and p ν R are small, where R is the nuclear radius, p e (p ν ) and l e (l ν ) are the momentum and orbital angular momentum of the outgoing electron (antineutrino) respectively. Therefore the larger angular momentum transfer would imply the smaller beta decay strength. This suggests that the largest weak decay happens for transitions with l e = 0 and l ν = 0, which, corresponding to the selection rules of ∆J π = 0 + and 1 + , are usually defined as the allowed Fermi and allowed Gamow-Teller (GT) transitions, respectively. Transitions with the larger angular momentum transfer are called the forbidden decays.
In this work we assume all transitions associated with fission fragments are allowed transitions and temporarily neglect the effects of forbidden decays § . Since reactor antineutrinos are produced from beta decays of the fission fragments with relatively large nuclear charges, and due to isospin symmetry and coulomb interactions between protons, their 0 + states are usually high-lying and outside the Q-value window. Thus here we only consider the allowed GT decays for the reactor antineutrino flux calculations.
For these allowed GT decays, we have the leading contribution of K = 1, L = 0 and s = 1, and the differential decay width can be calculated as [38,39]: where are the mass (spin, parity) of the initial and final nucleus respectively, p e and E e are the momentum and energy of the outgoing electron. J π f f ||σ||J π i i is the nuclear matrix element of the GT transition, where the double bras and kets refer to integrations over the solid angle only. is the lepton matrix element, where φ κe (Z) and φ κν are radial wave functions of the electron and antineutrino respectively, where Z is the nuclear charge, κ e and κ ν are defined as where i = (e, ν), j e and j ν are the total angular momenta of the electron and antineutrino respectively.
To proceed, one needs to calculate each of the nuclear and lepton wave functions, however, the nuclear wave functions are rather complicated and rely on the nuclear structure models. In order to factorize the nuclear and lepton matrix elements, a normal assumption is to treat the lepton radial wave functions as constants at the surface of the nuclear distribution, namely, the surface approximation (SA) taken as where R is the nuclear radius. This will simplify the differential decay width as where the GT transition strength is defined as: Notice that F κe,κν can be calculated using the antineutrino radial wave function j lν of the spherical Bessel function, and the electron radial wave functions f κe , g κe , which can be derived by using the following evolution equations [38,39]: where V (r) is the Coulomb potential of centrally distributed nuclear charges. The next step is to evaluate the lepton wave functions under the Coulomb potential, which can be calculated numerically, and in some special cases can be analytically expressed [37]. In this work, we employ the Radial package [36,41] and use Eq. (8) to calculate the exact numerical solutions for the radial electron wave function for any given nuclear charge distributions. In contrast, to obtain the well known analytical Fermi function and finite size effect, several additional approximations have to be implemented: • The neutrino long-wave approximation (νLA). Antineutrino produced from the nuclear beta decay can be described by the plane wave under the multi-pole expansions. The neutrino long-wave approximation is defined by taking the s-wave neutrino component to be unity and neglecting all other components: where j lν (p ν r) is the spherical Bessel function with the angular momentum l ν . Then the differential decay width will be reduced to where f 1 (p e R) and g −1 (p e R) are radial electron wave functions at the radius R with κ e = 1 and κ e = −1 respectively. • where γ = 1 − (αZ) 2 and y = αZE e /p e with α being the fine structure constant.
• The finite size effect. For the extended charge distributions in the nucleus, such as the uniform, the Gaussian or the Fermi charge distribution, besides the Fermi function of the point charge potential, there will additional terms to account for the finite size effect. There are many different analytical calculations of finite size corrections based on different methods and assumptions [42,43]. In this work, in order to compare with the exact numerical calculations, we choose the analytical finite size effect obtained by Wilkinson [43] using a uniform charge distribution within the nuclear radius R: with being the mass number of the nucleus. The definition of a n (n = −1, 0, · · · , 5) can be found in Ref. [43].
After all the above simplification, one can finally arrive at the famous analytical expression for the allowed GT transition: The antineutrino energy spectrum can be obtained by the direct replacement E e → Q − E ν . Notice that there exist some additional corrections [20] in Eq. (17), such as those of the screening effect, weak magnetism and radiative corrections, but to make our study easier, we temporarily neglect their contributions to the beta decay spectrum in this work.

Single beta spectrum investigation
Before calculating reactor antineutrino spectra with the ab initio method, we want to evaluate the accuracy of the analytical expressions derived in the previous section by comparing with the exact numerical solutions of the electron radial wave functions. Three different approximation effects, including the νLA, the analytical Fermi function and the finite size effect, will be tested by selecting four representative beta decay branches with different endpoint energies and nuclear charges, which are summarized in Table 1. Let us first consider the point charge distribution of the Coulomb potential and test the validity for the approximations of the νLA and the analytical Fermi function. The numerical calculations according to the Eq. (8) will be used as our benchmark to make the comparison, in which the GT transition strength is normalized to reproduce the experimental decay lifetime. Results of the νLA and the Fermi function are calculated using Eq. (13) and Eq. (14) respectively. Comparisons between the numerical calculations and analytical approximations of four representative decay branches are illustrated in Fig. 1, where the spectral ratios are defined as numerical calculations of the point charge distribution (PD) with (red) and without (blue) the νLA to the analytical Fermi function approximation (FF) for each beta decay branch. The solid and dashed lines are for the electron and antineutrino spectra respectively. Several comments are provided as follows: • The νLA is intended to approximate the radial neutrino wave function using the s-wave component of the orbital angular momentum. According to Fig. 1, one can observe that the deviation of the νLA from the exact numerical solution may distort the spectrum, which will reach around 1% for large Q branches, and at the level of 0.1% for small Q branches. This behaviour can be explained by the property of the spherical Bessel function. The s-wave component at the nuclear surface j 0 (p ν R) is always smaller than unity and it approaches to unity as p ν becomes smaller. Therefore νLA always overestimates the spectrum when the s-wave component is the dominate one. When p ν is large enough and approaching to the endpoint energy, the p-wave component may be sizable and contribute to the additional spectrum distortion near the endpoint energies.
On the other hand, there is no obvious effect for different choices of the nuclear charge Z, which can be simply understood by the neutrality of the neutrino, and its indirect effect is from the nuclear radius when applying the SA. and O(m e R). From Fig. 1, we can observe that the Fermi function overestimates the spectral strength in comparison to the exact numerical one, in particular for the high electron or low antineutrino energy regions. The deviation becomes larger as the electron energy increases (or equivalently the antineutrino energy decreases). This is reasonable because the Fermi function only includes the leading terms of p e R and the approximation would become worse when the electron momentum/energy becomes larger. This property is also applicable when we compare the left and right panels of Fig. 1, where larger endpoint energies Q will have larger spectral deviation. By comparing the upper and lower panels of Fig. 1, another parameter that may have significant contributions to the spectral deviation is the nuclear charge Z, where the deviation would be more significant for the branches with the larger nuclear charge. Therefore, it is clear that the spectral deviation of the Fermi function mainly comes from the endpoint energy Q and the nuclear charge Z, which can reach the magnitude of 20% for the branch with Q 8 MeV and Z 55.
Second, we want to discuss the results for the uniform charge distribution of the Coulomb potential and test the validity of the approximation for the finite size effect. The numerical   (17). Comparisons between the numerical calculations and analytical approximations of four representative decay branches are illustrated in Fig. (2), where the analytical finite size effects (blue) are defined as ratios of the analytical Fermi function (FF) with and without finite size corrections (FS), and the analytical finite size effects (red) are defined as ratios of the numerical spectra for the uniform charge distribution (UD) and point charge distribution (PD) without taking the νLA approximation. The solid and dashed lines are for the electron and antineutrino spectra respectively. The absolute sizes of both analytical and numerical finite size corrections show reduction of the spectral strength in most of the energy range, but the size and slope of analytical finite size corrections are much larger than those of the numerical ones. For the large Q and large Z branch the largest finite size correction will be around 10% for the analytical calculation, whereas it is only 3% for the numerical calculation. For the relative slope, the analytical and numerical finite size effects have similar behaviour that can be enhance the antineutrino spectra in the high energy part but reduce the spectra in the low energy part. However they are more significant for the analytical calculations than those of the numerical ones. From the general consideration of numerical calculations, the analytical Fermi function and finite size correction each has the spectral deviation of O(10%), thus it may be more severe when one considers both effects. However, the situation is the opposite. Since the total decay magnitude is always normalized to the experimental decay rate, it is thus more transparent to consider the relative spectral deviations of both analytical and numerical calculations. The Fermi function shifts more strength to the low antineutrino energies relative to the exact numerical one, whereas the analytical finite size correction shifts more strength to the high antineutrino energies relative to the numerical finite size correction, namely, compared to the exact numerical calculations, the spectral deviations induced by the Fermi function and the analytical finite size effect are in the opposite directions. In Fig. 3, the spectral ratios between the numerical spectra of the uniform charge distribution (UD) without the νLA and the analytical Fermi function (FF) with the finite size effect (FS) are illustrated for the representative branches. In the upper left panel, the largest spectral deviation would approach 15% and 6% for the low and high energy parts of the antineutrino spectra. For the relative deviation, it corresponds to a total spectral variation of 9% for the branch with Q 8 MeV  Finally, before finishing this section we would like to discuss the different effects of different extended charge distributions, where the Gaussian charge distribution (GD) and Fermi charge distributions (FD) are used to compare to the spectra of the uniform charge distribution (UD). The spectral ratios of different numerical calculations are illustrated in Fig. 4, where one can observe that the spectral deviations of GD and FD are at most at the levels of 2‰ and 0.2‰, respectively. Therefore, the effects of different extended charge distributions are pretty small and negligible, and it would be enough to use the uniform charge distribution in the future calculations of reactor antineutrino spectrum. namely, where Y k f (Z, A, t) is the activity of the f -th fission fragment with A and Z at time t, which converges to the cumulative fission yield and is independent of t after sufficient burning time, b f,i stands for the branching ratios of the transition associated with the endpoint Q i . E ν is the energy of the emitted antineutrino. Then the total antineutrino spectrum emitted by a reactor is determined by the summation of the contributions of all four fission fuels: in which α k is the fission fraction of the fission fuel k.
The spectrum associated with each of the fissionable isotopes 235 U, 238 U, 239 Pu and 241 Pu consists of more than 6000 beta decay transitions of the fission fragments. Thus in order to obtain the aggregate spectrum, we require not only the normalized single beta decay spectra, but also the corresponding branching ratios and the cumulative fission yields of corresponding fission fragments. In this work, we employ the latest nuclear database for these relevant nuclear data, where the cumulative fission yield data is taken from the Evaluated Nuclear Data File (ENDF) B-VIII.0 and the beta decay information is from the database of the Evaluated Nuclear Structure Data Files (ENSDF).
In previous ab initio calculations, analytical Fermi function and additional corrections are used to describe single beta decay spectrum. However, as we have demonstrated in the previous section, these analytical approximations are not always accurate, and the largest spectral deviation could be at the level of O(10%). Therefore it is necessary to study how they affect the isotopic antineutrino spectrum calculated from Eq. (17). In this work, we employ the numerical calculations of the radial electron wave functions in Eq. (8) with the uniform charge distribution and the SA approximation, and compare with those using the analytical calculation in Eq. (17). The numerical calculations of the radial electron wave functions are obtained from the Radial package [36,41].
In Fig. 5, ratios of the aggregate spectra between the numerical results obtained from exact calculations of the uniform charge distribution to the analytical ones obtained from approximate calculations of the Fermi function and finite size correction are illustrated for the fissionable isotopes 235 U, 238 U, 239 Pu and 241 Pu. Compared with the aggregate spectra using the numerical radial electron wave function, the analytical calculation would introduce a reduction of the low energy antineutrino spectra and an enhancement in the high antineutrino energy spectra. In the energy range between 2 MeV and 8 MeV, the antineutrino spectral deviations can reach around 2% to 3% at E ν = 8 MeV compared to the numerical calculations. Similar but opposite behavior can be observed for the electron energy spectra, where the spectral deviations of the analytical calculations are around 3% to 4% at E e = 8 MeV compared to the numerical calculations. Note that these spectral deviations are much smaller than those of individual beta decay spectra, which are affected by the small cumulative fission yields, small branching ratios and average effect between thousands of beta decay branches.  Moreover, one can observe that the behavior of spectral variations is rather similar for all the fissionable isotopes and their differences are within the 0.5% level, which indicates that if the single beta decay spectrum is really contribute to some levels of the experimental observed spectral excess, it should be universal for all the four fissionable isotopes.
Finally we would like to comment that the current numerical calculations are still taking the lepton wave functions at the nuclear surface in order to factorize the nuclear matrix element. In order to go beyond the SA approximation, one needs to know the radial dependence of the nuclear matrix element, which will rely on the nuclear many body calculations. Some preliminary studies [36] show that the allowed GT decays are dominated by the transitions between two single particle orbitals near the nuclear surface, but dedicated study on the nuclear structure information should be done for the precise beta spectral calculation. Moreover, even larger spectral deviation from the nuclear structure information has been observed for the first forbidden decays [44], which may have large contributions to the high energy range of the reactor antineutrino spectrum [22]. Therefore, further study on these directions would be interesting and constitutes the main goal of our future work.

Concluding Remarks
The appearance of the reactor flux and spectrum anomalies requires new accurate theoretical predictions for the reactor antineutrinos associated with the fission isotopes 235 U, 238 U, 239 Pu and 241 Pu. Many possible issues in the reactor antineutrino flux calculations have been investigated, but none of them can fully account for all these observed anomalies. In this work, we have proposed a new ab initio calculation of the isotopic reactor antineutrino fluxes with the exact calculation of the radial lepton wave functions, assuming all the decay branches are allowed GT transitions.
With the assumption of the SA approximation, the single beta decay spectrum is numerically calculated using the Radial package, and it has been demonstrated that the approximation of the analytical Fermi function and finite size effect each could have the largest spectral deviation of O(10%), whereas the effect of their combination could result in spectral deviations at the level of 5%-10%. Meanwhile, we also find that several forms of the extended charge distributions have negligible effects on the spectral variation. Using the state-of-theart nuclear databases, our new ab initio calculations have shown sizable but opposite spectral deviations at the level of 2%-4% in the antineutrino and electron spectra which may partially contribute to the observed spectral excess in the high energy range. We also find that the effect of analytical beta decay spectrum approximation is rather universal for all the four fissionable isotopes. Finally we want to emphasize that our new ab initio calculation based on the numerical wave functions can go beyond the SA and even be applied to the first forbidden transitions if nuclear many body calculations can be implemented, which will be the main goal of our future study in this respect.