Lorenzo Perissé

Anthony Onillon

Xavier Mougeot

Matthieu Vivier
email: matthieu.vivier@cea.fr

Thierry Lasserre

Alain Letourneau

David Lhuillier

Guillaume Mention

Lorenzo Périssé

reactor ν_e antineutrino fluxes and spectra

I. INTRODUCTION

The determination of the νe flux and spectrum emitted by a nuclear reactor is a long standing problem, dating back from the first detection of the neutrino in 1956 [1,2]. Reactor νe mostly originate from the β decay of neutronrich products following the fission of uranium and plutonium present in the nuclear fuel. They are also emitted in smaller amounts through the β decay of activation products following neutron irradiation of the fuel and/or structural materials present in the core. Past and present efforts focused on studying reactor νe both at commercial and research reactors. Commercial reactors, designed for the large scale production of electricity, mostly release energy through the fission of 235 U, 239 Pu, 241 Pu and 238 U isotopes while smaller research reactors run almost entirely through the fission of 235 U. The prediction of the νe flux and spectrum following the fission of each of these isotopes is a key ingredient to the study of the weak interaction and of the neutrino fundamental properties at such facilities, and still motivates huge experimental and theoretical efforts. The first attempts to predict a reactor νe spectrum come from the conversion of measured aggregate β spectra resulting from the thermal fission of 235 U [3,4]. This pioneering work led to the so-called conversion method, which was since then refined both with new experimental data [5][6][7][8][9] and an updated conversion procedure [10]. Other calculations were also attempted by adding up the theoretical prediction of all known fission fragment β/ν e spectra using their corresponding yields and decay schemes [11][12][13][14][15][16]. As pointed out from the very beginning, this summation method however mostly suffers from the incompleteness of the evaluated nuclear databases, especially regarding fission fragments having a large total decay energy Q β , hence giving a lower limit on the predicted νe flux. Until recently, the conversion method was therefore considered to provide the most accurate and robust predictions of νe fluxes and spectra at nuclear reactors.

In 2011, two independent revisions of the conversion procedure concluded to an increase in the expected reactor νe fluxes [17,18], leading to a (5.7 ± 2.4)% deficit in the detected νe rates with respect to predictions at short and middle baseline past reactor experiments [19]. This reactor antineutrino anomaly (RAA) was confirmed by the last generation of Inverse Beta Decay (IBD) experiments (Daya Bay [20], Double Chooz [START_REF] De Kerret | Double Chooz Collaboration)[END_REF], and RENO [22]) conducted at commercial reactors and aiming at measuring the θ 13 mixing angle. These experiments also highlighted an unexpected spectrum distortion in the 4.5 to 7.5 MeV energy regime with respect to predictions based on the Huber-Mueller (HM) revised calculations of the actinide fission νe spectrum of 235 U, 238 U, 239 Pu and 241 Pu [START_REF] Abe | Double Chooz Collaboration)[END_REF][24][25]. These anomalies motivated various experimental efforts at short baselines from research highly-enriched uranium (HEU) and commercial lowly-enriched uranium (LEU) reactors both to precisely reassess the νe absolute flux and spectrum as well as to search for a possible ∼1 eV mass sterile neutrino able to explain the RAA (see e.g. [26] and references therein). Although the sterile neutrino interpretation of the RAA is now strongly disfavored (see e.g. [27][28][29][30][31]), the discrepancies between predicted and observed νe fluxes and spectra have since been further confirmed, casting doubts on the reliability of the state-of-the-art conversion predictions. In particular, the absolute normalization of the ILL aggregate β spectra they rely on has recently been seriously questioned [START_REF] Onillon | Proceedings for the 14th installment of Applied Antineutrino Physics (AAP) workshop series[END_REF]. This suspicion was later confirmed by a measurement of the ratio of the 235 U over the 239 Pu aggregate β spectra, indicating for a (5.4 ± 0.2)% excess with respect to that same ratio constructed with the ILL reference data [START_REF] Kopeikin | [END_REF]. At the same time, the RAA and the observed spectral anomalies have motivated many efforts in the nuclear physics community to improve the summation calculations, as this method is a powerful tool for dissecting and assessing systematic effects in reactor νe predictions [34][35][36]. Particular focus have been paid on evaluating the role of the non-unique forbidden transitions [37][38][39][40] and of the quality and incompleteness of the evaluated nuclear databases [41][42][43][44]. One notable finding is that the so-called Pandemonium effect, leading to a biased estimate of many of the fission fragment β decay schemes in modern nuclear evaluated data libraries, still contributes to discrepancies between summation predictions and state-of-the-art measurements of the νe flux at nuclear reactors [45].

Although significant progress has been made, especially with many improvements in the content and the quality of the nuclear evaluated data libraries, the summation method still suffers from approximations and lacks a comprehensive uncertainty budget to truly be a robust prediction tool. Together with the forthcoming generation of experiments aiming at studying Coherent Elastic Neutrino-Nucleus Scattering (see e.g. [START_REF] Abdullah | Contribution to Snowmass 2021[END_REF]), the current and next generations of IBD experiments will reach increasing precision in the measurements of νe fluxes and spectra at reactors, making it necessary to extend and refine the summation calculations. The following article therefore reports about a careful revision of the summation method, improving over current stateof-the-art predictions by using an advanced modeling of β decay and recent evaluated nuclear data. This revision work is based on the software Beta Energy Spectrum Tool for an Improved Optimal List of Elements (BESTIOLE), which was developed for the reevaluation of reactor antineutrino spectra in [17]. Further to these improvements, the largest systematic effects known to impact the calculations of reactor νe fluxes and spectra are here quantitatively studied to propose for the first time a detailed uncertainty model. As a cautionary word, it should already be stressed here that the construction of such an uncertainty model remains a difficult task which has to address both the limitations (quality and incompleteness) of the evaluated nuclear data libraries and the many approximations used in the modeling of the β branches. Especially when the necessary input information are incomplete or missing, the uncertainty associated with each systematic effects is estimated following simple approaches using known and reliable data as a proxy, with the prime care of remaining conservative and as realistic as possible. The robustness and the limitations of this first uncertainty model for the summation method are in this regard discussed and criticized all along the article.

The article is structured as followed: section II lays down the general principles of the summation method along with the uncertainty propagation formalism used throughout this work. Section III details the improvements the present work makes over past modelings of β decay. In particular, the implementation of various electro-magnetic corrections to the Fermi theory as well as realistic calculations of the main forbidden non-unique transitions contributing to the detected flux at a reactor are presented. Section IV describes the construction of a new base of nuclear input data for the modeling of each β-branch contributing to a reactor νe spectrum, using recent Pandemonium-free evaluations available in the literature and in online databases. Following these major revisions, section V then presents the new summation calculations of the four major actinides νe fission spectra ( 235 U, 239 Pu, 241 Pu, 238 U) along with a detailed breakdown of their uncertainties. In light of the RAA and the observed spectral anomalies, these new results are challenged against state-of-the-art predictions and measurements. Finally, section VI summarizes the main results of this work and opens up possible improvements to the summation method, especially in view of addressing the RAA and providing a complete and robust prediction tool for the next generation of reactor νe experiments.

II. THE SUMMATION METHOD

The summation calculation of a reactor β spectrum [47] consists in adding up the β spectrum S p of any fragment p populating the fission process of any of the four major actinides k = 235 U, 239 Pu, 241 Pu, 238 U present in the core:

S tot (E, t) = k,p A k p (t) S p (E), (1)
where A k p (t) is the activity of the p th fission fragment at irradiation time t and E the β particle kinetic energy. The fission fragment activities can be estimated using their cumulative yield Y k p (t) after an irradiation time t, so that eq. 1 can also be written as:

S tot (E, t)= k f k (t) p Y k p (t) S p (E) (2) = k f k (t) S k (E, t), (3)
where f k (t) is the k th actinide fission rate. The term S k (E, t) is here the total β spectrum associated to one fission of the k th actinide, and is usually denoted as actinide fission spectrum. The β spectrum S p of a fission fragment can generally be broken down to the superimposition of N p transitions connecting either the ground state (GS) or an isomeric state (IS) to different excited levels b of its daughter nucleus:

S p (E) = N p b=1 B p b S p b (E, E p 0,b ). ( 4
)
The probability B p b to decay through a given transition is denoted here branching ratio. The transition branching ratios add up to unity if the fission fragment is a pure β-emitter and to less than one otherwise. The endpoint energy E p 0,b of a transition is usually expressed as:

E p 0,b = Q p β + E p IS -E lvl b , (5)
where Q p β is the total β decay energy (corresponding to a ground state to ground state transition), E p IS is the energy of the parent nucleus isomeric state (E p IS = 0 if the parent nucleus β decays from the ground state), and E lvl b is the energy of the daughter nucleus b th excited state. Combining eq. 1, 4 and 5, a reactor β spectrum can hence be considered as the sum of many individual βbranch spectra, each requiring both an accurate β-decay formalism and robust evaluated nuclear data to be modeled. The correspondence between a reactor β and νe spectrum is achieved at the β-branch level using energy conservation and neglecting the daughter nucleus recoil energy by substituting in the above formula the β particle kinetic energy E with the νe energy E ν = E p 0,b -E. As opposed to the conversion method, the summation method then ensures a true and unique correspondence between reactor β and νe spectra while computing them with equal precision.

In the following, many results will be expressed not only through raw νe flux and spectrum calculations, but also computing IBD yields in order to ease the comparison with state-of-the-art νe flux predictions and measurements (see e.g. section V). The IBD yield is here computed as:

σ IBD = Emax E th σ IBD (E ν ) S(E ν ) dE ν , (6)
where σ IBD (E ν ) is the IBD cross-section taken from [START_REF] Vogel | [END_REF] and S(E ν ) is a νe spectrum. Throughout this article, IBD yields could respectively be computed either at the level of β branch, at the level of a fission fragment or for a full actinide fission spectrum. The normalisation of the IBD cross-section was estimated using the Particle Data Group 2022 evaluation of the required fundamental constants [49], also including radiative corrections as derived in [50]. The lower integration bound in eq. 6 was set to the IBD energy threshold E th = 1.806 MeV. The upper integration bound is especially relevant for the computation of the IBD yields associated to the fission of the 235 U, 238 U, 239 Pu and 241 Pu actinides, denoted hereafter as isotopic IBD yields. It was chosen as E max = 10 MeV accordingly to [51]. Extending this bound to energies higher than 10 MeV has been checked to give negligible (<0.1%) changes.

Figure 1 illustrates the summation calculation of a reactor νe spectrum originating from the thermal fission of 235 U. A fission spectrum typically piles the individual β/ν e spectra of ∼800 fragments, and totals more than 10,000 β branches. The exponential decrease of the spectrum as a function of energy results from the underlying Q β distribution of the fission fragments, which makes fewer and fewer β branches to contribute toward high energies. Interestingly, only a handful of fission fragments dominates the associated IBD yield. For instance, about thirty fission fragments exhibit a >1% contribution to the 235 U IBD yield. They are indicated by the red solid lines on Figure 1.

FIG. 1. Summation calculation of the νe spectrum resulting from the thermal fission of 235 U. The total spectrum (black solid line) is broken down to the contributions of all the fission fragments listed in the present day nuclear databases (grey solid lines). The dark grey area, especially visible at high energies, represents the associated total uncertainty. The red solid lines highlight the νe spectrum of 32 fission fragments, each contributing more than 1% of the expected IBD yield (see eq. 6). Summed all together, these fission fragments total 60% of the 235 U isotopic IBD yield.

The total uncertainty associated to the summation calculation of a reactor β/ν e spectrum combines uncertainties from many sources, which are propagated using the following covariance matrix formalism [START_REF] Mathai | Jacobians of matrix transformations and functions of matrix argument[END_REF]. The covariance matrix V s p associated to the binned spectrum S p (E) of a fission fragment p for an uncertainty source s (see eq. 4) is given by:

V s p = i V s i + i,j i =j V s ij , (7)
where the first term sums the covariance matrix of the i th branch spectrum for the parameter s, and the second term sums the cross-term covariance matrices between the spectra of the i th and j th branches. Such cross-term covariance matrices encode correlations between branch spectra sourcing from correlations between parameters associated to one type of uncertainty source. The uncertainty sources considered in the present work are listed in Table VIII. They were defined and sometimes grouped such that they are independent, meaning that the total covariance matrix V p of a fission fragment spectrum can be expressed as:

V p = s V s p , (8)
where the parameter s runs over the different uncertainty sources. The total uncertainty associated to each bin of a fission fragment spectrum is then given by the squared root of the V p diagonal elements. Applying the Jacobian matrix formalism to eq. 2 and eq. 3 for the propagation of uncertainties, the total covariance matrix V k associated to an actinide fission spectrum S k (E) can be approximated at first order as:

V k p,q Y k p Y k q V pq + p,q V Y k p Y k q S p (E) S q (E). (9)
where V pq is the covariance of the binned spectra of the p th and q th fission fragments. As detailed later in section III, correlations between two fission fragment spectra S p and S q originate from common sources of uncertainty in their respective modeling. Finally, the term V Y k p Y k q represents the covariance of the p th and q th fission fragment cumulative yields. The possible sources of correlation among cumulative fission yields are especially discussed in section IV A. In eq. 2 and eq. 3, fission rates f k are usually estimated using reactor simulations. Modeling and nuclear data uncertainties can induce correlations between the estimated f k . Furthermore, interactinide correlations also exist because a fission fragment can contribute to different actinide fission spectra. As such, the total covariance matrix of a reactor spectrum S tot (E) is computed at first order as:

V tot k,l V f k f l S k (E) S l (E) + k,l =k p,q f k f l Y k p Y l q V pq + k f 2 k V k , (10)
where V f k f l is the covariance of the k th and l th actinide fission rates f k and f l . As in eq. 9, the covariance V pq comes from a common source of uncertainty both present in the p th and q th fission fragment spectrum modeling.

In the following work, the construction of the uncertainty budget distinguishes two main classes of uncertainty. The first one gathers uncertainty sources from the modeling of the β branches whereas the second one is related to evaluated nuclear data.

III. BETA DECAY FORMALISM

A. Beta branch modeling

In the present work, the electron spectrum of a β branch is modeled according to the (V-A) theory of weak interaction using an advanced formalism developed by Behrens and Bühring [START_REF] Behrens | Electron radial wave functions and nuclear beta-decay[END_REF]. The most general form of the β-decay Hamiltonian is usually expressed as the product of a lepton current encompassing all information related to the electron and neutrino wave functions and of a nuclear current encoding all the nuclear structure information of the parent and daughter nuclei. In this formalism, a multipole expansion of both these currents is performed to compute the transition matrix elements associated to a nuclear β decay, assuming a spherical symmetry of the system.

The transition matrix element is decomposed such that the lepton kinematic dependency is clearly separated from the pure nuclear structure term. This expansion introduces form factors that can be either of vector (V) or axial-vector (A) type. Their contribution in the transition depends on selection rules related to total angular momentum conservation. Each order of this multipole expansion is associated to a change l in angular momentum between the initial and final states. Usually, only terms with the lowest powers, i.e. l and (l+1), are kept in this expansion, the higher order terms being expected to be orders of magnitude smaller. The spin ∆J = |J f -J i | and parity π f π i changes between the initial and final nuclear states then define which multipole expansion terms contribute to the β-decay probability, leading to the wellknown β-decay selection rules [START_REF] Behrens | Numerical Tables for Beta-Decay and Electron Capture (Landolt-Börnstein New Series[END_REF]. In particular, transitions leaving out a single nuclear form factor in their multipole expansion are called unique transitions, whereas those having several nuclear form factors are called nonunique transitions.

The electron and neutrino wave functions in the lepton current are calculated taking into account the electromagnetic interaction of the outgoing β particle with the static Coulomb potential of the daughter nucleus. They are also expanded in terms of spherical harmonics. The combination of the previously described multipole decomposition of the transition matrix element and the available phase space then leads to the following usual expression for the electron spectrum:

S b (W e , W 0 ) = K p e W e (W 0,b -W e ) 2 × F(Z, W e ) × C(Z, W e ) × (1 + δ e - R + δ WM ), (11)
where W e = E e /m e c 2 + 1 and p e = W 2 e -1 are respectively the electron total energy in units of its rest mass and associated momentum, Z is the atomic number of the daughter nucleus and K a factor normalizing the transition probability to unity. The term F(Z, W e ) denotes the Fermi function. It encodes the distortion of the electron wave function in the static Coulomb potential of the daughter nucleus. In the Behrens and Bühring formalism, the Fermi function is defined as:

F(Z, W e ) = α 2 -1 + α 2 +1 2p 2 e , (12)
where the α k quantities are called Coulomb amplitudes and are related to the normalization of the electron radial wave functions (ERWFs), the integer k being related to the angular momenta as defined in [START_REF] Behrens | Electron radial wave functions and nuclear beta-decay[END_REF]. The calculation of the Coulomb amplitudes is detailed in section III B.

The term C(Z, W e ) in eq. 11 is called the shape factor and includes the nuclear form factors originating from the multipole expansion of the nuclear current. At first order, the shape factor C U L (Z, W e ) of a unique-forbidden transition of degree L can be expressed in a rather simple way since the only contributing nuclear form factor can be factorized out into the normalisation constant K:

C U L (Z, W e ) = L k=1 λ k p 2(k-1) e p 2(L-k) ν (2k -1)! (2(L -k) + 1)! . (13)
The λ k parameters are called Coulomb functions and are related to the Coulomb amplitudes α k through the following relationship:

λ k = α 2 -k + α 2 +k α 2 -1 + α 2 +1 , (14)
while p ν stands for the neutrino momentum. The computation of the shape factor for non-unique forbidden transitions is in the other hand more complicated. It usually requires information about the structure of the initial and final nuclear states. The treatment of the non-unique forbidden transitions is therefore discussed separately in section III C. Finally, the terms δ e - R and δ WM respectively correspond to radiative and weak magnetism corrections. They are detailed in section III D and in section III E.

It is worth mentioning at this stage that the present β-decay formalism includes many refinements with respect to past summation calculations of reactor νe spectra [17,18,41]. First, the calculation of unique forbidden transitions, as described in eq. 13, avoids the usual λ k = 1 approximation. This approximation originates from the fact that the Coulomb functions λ k are almost constant at high energy. However, it was shown to be systematically incorrect compared to unique forbidden shape factor calculations including the full energy dependence of the λ k parameters [START_REF] Mougeot | [END_REF]. Second, the present formalism automatically includes the so-called nucleus finitesize and atomic screening effects through the calculation of the Coulomb amplitudes α k (see section III B). As opposed to past calculations, these electromagnetic corrections are here not only included in the calculation of the Fermi function, but also in the calculation of the shape factor. The impact of these two refinements is discussed in section III B for the calculation of unique forbidden shape factors. Last, the present summation calculations operate a new treatment of the non-unique forbidden transitions. As described in section III C, detailed nuclear structure calculations are used to model the 23 most important non-unique branches contributing to the IBD yield detected at a nuclear reactor (see Table II), hence avoiding the systematic use of the socalled ξ-approximation to compute their corresponding shape factor.

B. Relativistic electron wave functions

The Coulomb amplitudes α k are obtained by solving the Dirac equation for the electron in the static Coulomb potential of the daughter nucleus. The algorithm for the numerical solving of the Dirac equation follows the work of Behrens and Bühring [START_REF] Behrens | Electron radial wave functions and nuclear beta-decay[END_REF], and uses local powerseries expansions of the ERWFs. The algorithm iteratively solves the Dirac equation on a spatial grid, starting from r = 0 up to a reconnection point r = R 2 where the free electron wave functions associated to r → +∞ can be safely recovered. The implementation of this algorithm has been validated against the published tables of Behrens and Jänecke [START_REF] Behrens | Numerical Tables for Beta-Decay and Electron Capture (Landolt-Börnstein New Series[END_REF] up to the last possible decimal for an electrostatic potential generated by a uniformly charged sphere. Any reasonable change in the numerical algorithm parameters such as the choice of the cut-off in the power-series expansion of the ERWFs, the choice of the reconnection point R 2 or the choice of the spatial grid size have been checked to give negligible changes in the value of the Coulomb amplitudes α k [START_REF] Perisse | Modeling of reactor antineutrino spectra[END_REF]. Numerical calculation errors coming from the solving of the Dirac equation are then safely neglected in the computation of a β-branch spectrum.

In this work, the so-called nucleus finite-size and screening corrections to the calculation of a β-branch spectrum are taken into account by considering a nuclear potential V(r) generated by Z charges uniformly distributed within a spherical nucleus of radius R and screened by a cloud of (Z -1) atomic electrons. The potential V(r) takes different forms in the regions respectively delimited by the nuclear radius R and a point R 1 beyond which the screened potential is smoothly transited up to an asymptotic point-like form until reaching the reconnection point R 2 . The complete expression of the screened potential used in this work can be found in [START_REF] Mougeot | [END_REF]. In the region R ≤ r ≤ R 1 , the screened potential is expressed as:

V(r) = - αZ r N i=1 a i e -βir , (15)
where the parameters a i and β i are tabulated from fits of atomic screening functions calculated using the rel-ativistic Dirac-Hartree-Fock-Slater formalism for atoms with Z=1-92 [58]. These screening functions improve over past evaluations, which often relied on the nonrelativistic Thomas-Fermi statistical model of the atom. Although no uncertainties are reported in the literature, they have been shown to provide a very accurate modeling of screened potentials as compared to e.g. electronatom scattering experimental data [59]. As such, no uncertainty in the screened nucleus potential for the calculation of the Coulomb amplitudes α k is considered in this work. The nuclear radius R, also entering the definition of the nuclear Coulomb potential V(r), is either taken from a set of experimental root mean square charge radii evaluated in [60] or estimated using the prescription from [START_REF] Angeli | Table of nuclear root mean square charge radii[END_REF] for nuclei off the stability line. This prescription slightly improves over the Elton formula commonly used in previous summation calculations of reactor νe spectra [17,41] as it better reproduces experimental data for isotopes off the stability line. This new modeling of the nuclear radius has been found to give a negligible difference with respect to summation calculations using the Elton formula [START_REF] Perisse | Modeling of reactor antineutrino spectra[END_REF]. As such, the contribution of the evaluated nuclear radius uncertainties reported in [60], which are smaller than this difference, are safely neglected.

For illustration purposes, a set of transitions with different forbidenness degrees is compared in Figure 2 (a). Changes to unique forbidden shape factors when considering an appropriate calculation of the λ k functions both including nuclear finite-size and atomic screening effects are depicted in Figure 2 (b). Compared to past calculations defaultly applying the λ k = 1 approximation, the shape factor of a unique forbidden transition with a 10 MeV endpoint energy is typically corrected by a O(5%) factor, with even larger corrections for smaller endpoint energy transitions. For instance, unique transitions having an endpoint energy close to the IBD threshold can be corrected up to 60% [START_REF] Perisse | Modeling of reactor antineutrino spectra[END_REF]. Although these corrections are significant at the β-branch level, they are expected to bring small changes to the computation of an actinide fission spectrum since unique forbidden transitions typically contribute to ∼10% of both the corresponding νe and IBD yield (see Table I andFigure 3).

C. Treatment of non-unique forbidden transitions

Non-unique forbidden transitions play an important role in the calculation of reactor νe spectra. For instance, Table I details their contribution to the 235 U, 238 U, 239 Pu and 241 Pu νe fluxes and associated IBD yields, showing that they make a significant 30 to 40% contribution to the latter. As shown by Figure 3, these transitions also dominate the 4-8 MeV portion of a typical actinide fission spectrum. A precise calculation of their contribution is then relevant in light of the recently measured spectral distortions in that same energy range. As mentioned earlier in section III A, non-unique forbidden transitions leave out several form factors in their associated β spec-trum. Their modeling requires the use of a mathematical formalism based on nuclear structure calculations, which are complicated and numerically time-consuming to perform. This section describes how the non-unique forbidden transitions are treated in the present work. Because nuclear structure calculations are computationally heavy, they were used to model the 23 most important nonunique branches participating to the typical IBD yield expected at a nuclear reactor (see Table II). This strategy ensured to compute a significant ∼70% fraction of the total non-unique forbidden transition contribution to the IBD yield, hence providing a good compromise with respect to the overall accuracy of a reactor νe flux prediction. For perspective, the computation of about 70 (resp. 400) additional branches would be necessary to describe 90% (resp. 99%) of this contribution. In section III C 2, the nuclear structure calculations of these 23 transitions are then used to construct an uncertainty model associated to the remaining non-unique transitions, which in the present work still miss realistic nuclear structure calculations and are computed using the ξ-approximation.

Calculations with nuclear structure

In the Behrens and Bühring formalism [START_REF] Behrens | Electron radial wave functions and nuclear beta-decay[END_REF], the shape factor C(Z, W e ) in eq. 11 results from the previously mentioned multipole expansion of both the nuclear and lepton currents, assuming the so-called impulse approximation in which all nucleons are independent particles in a mean-field potential. The convolution of the nuclear structure and the lepton dynamics is embedded in the ad hoc terms M K (k e , k ν ) and m K (k e , k ν ) and is decomposed as a sum of the different multipole orders, giving the following general expression for the shape factor:

C(Z, W e ) = K,ke,kν λ ke M 2 K (k e , k ν ) + m 2 K (k e , k ν ) - 2µ ke γ ke k e W e M K (k e , k ν ) m K (k e , k ν ) , (16)
where any quantities labeled by the lepton quantum numbers k e and k ν depend on the electron and neutrino relativistic wave functions, respectively. The main multipole order K comes from the expansion of the nuclear current and is limited by the change in total angular momentum ∆J between the initial and final nuclear states. It ranges from K min = ∆J to K max = J i + J f . In the present work, the prescription from [START_REF] Behrens | [END_REF] has been followed, considering only the dominant terms with respectively K = K min , K min + 1 and k e + k ν = K + 1, K+2 in eq. 16. The next order terms for a set of high-energy transitions of interest were computed so as to check the validity of this prescription. Their contribution was found to be several orders of magnitude smaller, thus confirming they could be neglected. [START_REF] Plompen | [END_REF] (see section IV A), and with nuclear structure and decay data as described in section IV B. The 1 st non-unique forbidden line also displays in parenthesis the contribution of the 1 st non-unique forbidden transitions computed with nuclear structure calculations (see section III C 1 and Table II). The contribution of the nuclides having no decay information (see section IV D) is also displayed here for completeness, such that all contributions add up to 100%.

In addition to the multipole expansion, the lepton wave functions are also expanded in the Behrens and Bühring formalism in powers of (m e R), (W e R) and (αZ), with R the nuclear radius and α the fine structure constant. This procedure avoids any overlap calculations between the nuclear and lepton wave functions. The nuclear ma-trix elements, also called form factor coefficients, become then independent of the lepton momenta, making the computation of the shape factor extremely fast, about a few milliseconds on a modern computer. Using this approach, the dominant term M K (k e , k ν ) in a non-unique forbidden transition reads: FIG. 3. Contributions of the different types of β -transition to the 235 U and 239 Pu fission νe spectra. Computations are done using cumulative fission yields from JEFF-3.3 [START_REF] Plompen | [END_REF] and using nuclear decay data as described in section IV B. The solid line contributions add up to 100%. The contribution from the 23 1 st non-unique forbidden transitions modeled with nuclear structure calculations (dashed blue line) belongs to the 1 st non-unique forbidden transition contribution (solid blue line). Contributions of other unique forbidden transitions (solid red lines) lie below 0.1% and are not visible.

M K (k e , k ν ) = 1 2 (2K)!! (2K + 1)!! 1 (2k e -1)!(2k ν -1)! (p e R) ke-1 (p ν R) kν -1 × - 2K + 1 K V F (0) K,K-1,1 - αZ 2k e + 1 V F (0) K,K,0 (k e , 1, 1, 1) - W e R 2k e + 1 + p ν R 2k ν + 1 V F (0) K,K,0 - αZ 2k e + 1 K + 1 K A F (0) K,K,1 (k e , 1, 1, 1) - W e R 2k e + 1 - p ν R 2k ν + 1 K + 1 K A F (0) K,K,1 . (17)
The form factor coefficients V/A F (N) K,L,s (k e , m, n, ρ) are either of vector or axial-vector type, as indicated by the left upper script. The L and s orders come with K in the development of the nuclear current and their contributions are summed in M K (k e , k ν ). The other dependencies (N) and (k e , m, n, ρ) label the radial expansion of the lepton wave functions. When ρ = 0, the dependency in (k e , m, n, ρ) disappears and is omitted in the notation, as shown in eq. 17. More details can be found in [START_REF] Behrens | [END_REF]. Such an expansion of the lepton wave functions is only possible for a simple Coulomb potential, e.g. originating from a nucleus modeled as a uniformly charged sphere. In the present work, a different strategy is followed. The formulation of the M K (k e , k ν ) and m K (k e , k ν ) quantities is revisited in order to directly use the numerical ER-WFs described in section III B. This approach ensures more accurate calculations of the shape factor because (i) both the finite-size nucleus and atomic screening effects are inherently taken into account in the calculation of the lepton wave functions and (ii) the lepton wave functions are not expanded anymore. The multipole expansion therefore only remains, the precision of which has been checked to be under control. The cost of such a full numerical treatment for the lepton current is the computational burden, each branch spectrum requiring several tens of minutes to be calculated.

The form factor coefficients

V/A F (N)
K,L,s can be calculated in a very simple way by reducing the β transition to a single nucleon-nucleon transition. However, it is clear that such a description of the nuclear structure is neither realistic nor accurate. In the configuration mixing approach, the many-particle wave function of a nuclear state is described as a linear combination of the singleparticle wave functions. The β decay transition ampli- II. Non-unique forbidden transitions computed with realistic nuclear structure calculations (see section III C 1) and listed by decreasing order of importance according to their relative contribution to the total IBD yield expected at a commercial nuclear reactor. The following fission fractions were assumed: 0.559 for 235 U, 0.088 for 238 U, 0.291 for 239 Pu and 0.062 for 241 Pu. Decay data are extracted from the ENSDF database [START_REF] Nndc | From ENSDF database as of June 2020[END_REF] and corrected with up-to-date TAGS data when indicated by † (see section IV C 1). The last three columns respectively indicate the contribution of these transitions to the total νe flux, the contribution of these transitions to the total IBD yield, and the difference between the transition IBD yield computed with nuclear structure calculations and under the ξ-approximation, expressed relatively to the former. These transitions add up to 69.1% of the total non-unique forbidden transition contribution. They total 26.9% of the expected IBD yield.

Nuclide Q β EIS E lvl BR J π i → J π f φν e φIBD ∆φIBD [MeV] [MeV] [MeV] [%] [10 -1 %] [%] [%]
tude is obtained by evaluating the corresponding onebody spherical tensor operator T λ bracketed between the initial and final nuclear states, with λ the tensor rank. This total amplitude can be expressed as a weighted sum of all the single-particle transition amplitudes that play a role in the β transition [START_REF] Suhonen | From Nucleons to Nucleus: Concepts of Microscopic Nuclear Theory[END_REF]:

ξ f J f |T λ |ξ i J i = 1 √ 2λ + 1 × a,b a|T λ |b ξ f J f | c † a cb λ |ξ i J i . ( 18
)
The single-particle matrix element a|T λ |b describes a nucleon-nucleon transition, weighted by its one-body transition density (OBTD)

ξ f J f | c † a cb λ |ξ i J i .
Eventually, each nucleon-nucleon transition gives a specific form factor coefficient and all the contributions have to be summed up to determine M K (k e , k ν ), as e.g. in eq. 17.

In this work, the OBTD for each nucleon-nucleon transition is computed using the shell model numerical code NuShellX@MSU [START_REF] Brown | [END_REF]. This code can use different interaction Hamiltonians, each fitted in different mass regions to reproduce some experimental data. Nuclear structure wise, nucleons are distributed among a set of low-lying energy levels defining an inert core and a set of high-lying energy levels defining a valence space. Only nucleons present in the valence space can experience a transition in such a model. The proper definition of the inert core, the valence space and the interaction Hamiltonian depend on the nucleus mass number A. For nuclei with A <100, the glepn valence space above the doubly-magic 56 Ni core as well as the recommended interaction with identical name is used [67]. In this scheme, the valence nucleons fill the shells starting from 2p 3/2 up to 2d 3/2 . The valence space is constrained such that the computation of the very large number of configurations remains tractable. It was thus restricted around the magic number 50: proton shells were free up to 1g 9/2 and higher shells were forced to be empty. Neutron shells were forced to be full up to 1g 9/2 and higher shells were free. For nuclei with A >100, the jj56pn valence space above the doubly-magic 132 Sn core along with the recommended interaction khhe adapted from [68] were selected. For nuclei with A ≤ 140, no restriction was necessary on the valence space. For nuclei with A=141 and A=142, the valence space was constrained as follows, with (min,max) number of particles in a given shell: for protons, (2,4) in 1g 7/2 and (0,2) in the other shells; for neutrons, (1,4) in 1h 9/2 , (0,4) in 1i 13/2 , and (0,2) in the other shells. For nuclei with A=144, the proton shell were free except 1h 11/2 that was forced to be empty, and a maximum of three neutrons was allowed in all the orbitals except 1i 13/2 that was also forced to be empty. Once the nuclear levels determined, the OBTD have been calculated for the dominant K values in each β transition.

The single-particle matrix elements entering eq. 18 are computed using nucleon wave functions extracted from harmonic oscillator potentials. For a realistic estimate of their corresponding frequencies, the method depicted in [69] using the proton configurations provided by NuShellX was followed. When available, experimental root mean square charge radii have been taken from [60]. Otherwise, the fitted radius formula off the stability line from [START_REF] Angeli | Table of nuclear root mean square charge radii[END_REF] was used. At this stage, it should be reminded that the present β-decay formalism is totally relativistic. The lepton and nucleon wave functions must thus exhibit large and small components resulting in non-relativistic transition matrix elements combining only large components, and in relativistic transition matrix elements involving small components. The expansion of M K (k e , k ν ) depicted in eq. 17 shows that the first form factor coefficient to arise is the relativistic vector matrix element V F K,K-1,1 , for which an accurate value is therefore of importance. Most of the nuclear structure models, including NuShellX, are however non-relativistic. To circumvent this issue, a simple solution consists in assuming that the large component of the nucleon wave function corresponds to the non-relativistic harmonic oscillator. An estimate of the small component is then deduced from the large one using the non-relativistic limit of the Dirac equation. Nevertheless, the inaccuracy of such an approach has been pointed out for decades (see e.g. [70]). The best current approach consists in using the conserved vector current hypothesis (CVC), a property emerging from the gauge invariance of the weak interaction. The relativistic form factor coefficient V F K,K-1,1 can then be related to the non-relativistic form factor coefficient V F K,K,0 by [START_REF] Behrens | Electron radial wave functions and nuclear beta-decay[END_REF]:

V F K,K-1,1 - R K(2K + 1) V F K,K,0 × W 0 -(m n -m p ) + ∆E C , (19)
where m n and m p are respectively the neutron and proton rest masses, W 0 is the transition endpoint energy and ∆E C is the Coulomb displacement energy. The Coulomb displacement energy can be estimated by modeling the nucleus as a uniformly charged sphere:

∆E C = 6 5 αZ R . ( 20
)
This simple modeling is only an approximation and ∆E C was demonstrated to be sensitive to the mismatch between the initial and final nucleon wave functions [71].

The present calculations rather use the prescription from Behrens and Bühring, which assumes that the singleparticle potential difference is determined by the average of the Coulomb potential [START_REF] Behrens | Electron radial wave functions and nuclear beta-decay[END_REF]. The average of the Coulomb potential involves here the full numerical treatment of the lepton current and the nucleon wave functions as described previously. This Coulomb energy evaluation, which depends on the electron kinetic energy and is specific to each nucleon-nucleon transition, is expected to be more accurate than the usual formula given in eq. 20.

Finally, the axial-vector coupling constant g A , which appears in the definition of axial-vector form factor coefficients, may need to be adjusted to reproduce some experimental observables. Nucleon-nucleon transitions occur in nuclear matter, and a quenched value of g A can partially correct for a mismodeling of the nuclear structure such as for instance the approximate treatment of the many-nucleon correlations. Varying g A has been shown to modify the spectrum shape of some forbidden nonunique transitions, mostly for low-energy transitions [72][73][74]. Without experimental spectrum to compare with, an effective g A value can be estimated from the quenching factor in infinite nuclear matter [75]. However, the accuracy of this value is not guaranteed. Because the forbidden non-unique transitions of interest reported in Table II have large endpoint energies, a free-nucleon value of g A = 1.2763 (15) resulting from the mean of two recent precise measurements [76,77] has been chosen [78].

Calculation of the non-unique forbidden transitions as listed in Table II and using the previously described nuclear structure formalism is computationally heavy. The uncertainty associated to these transitions is therefore modeled in a rather simplistic but conservative way, by using the spectrum difference between the accurate modeling version and an allowed version of the transition obtained under the ξ-approximation (see section III C 2). For each energy bin, the accurate non-unique forbidden and the ξ-approximated spectra are used to define the limits of a uniform distribution. The standard deviation of this distribution is then used to construct an associated covariance matrix. This covariance matrix is also constructed such that the total rate of the branch spectrum is conserved. Energy bins showing equal (resp. opposite) sign in the difference between the accurate modeling and the ξ-approximated version of the spectrum are treated as fully correlated (resp. anti-correlated). With this uncertainty modeling, the IBD yield associated to a non-unique forbidden transition calculated with nuclear structure typically shows a O(10%) uncertainty. The uncertainties associated to these 23 non-unique forbidden transitions computed with nuclear structure induce an uncertainty of ∼0.2% on the isotopic IBD yields when they are assumed to be uncorrelated (see Table VIII). The contribution of this source of uncertainty at the actinide fission spectrum level is also depicted on Figure 6 for 235 U. The impact of fully correlating these 23 nonunique transitions has been checked, and resulted in a ∼0.4% uncertainty on the isotopic IBD yields. Although theoretical correlations may exist, no evidence allows to suggest that these branches are fully correlated. As such, the present work keeps these branches uncorrelated in the propagation of their corresponding uncertainties. The computation of the remaining non-unique forbidden transitions is done using the ξ-approximation and is discussed in the next section.

ξ-approximation

The accurate treatment of the 23 most important nonunique forbidden transitions as described in the previous section cannot be applied on a case-by-case basis to all the remaining hundreds of non-unique transitions present in a reactor νe spectrum. These remaining non-unique forbidden transitions were instead estimated using the ξ-approximation. The ξ-approximation consists in treating a non-unique forbidden transition as its equivalent in term of total angular momentum variation but disregarding any parity change. This means that the shape factor of a L th non-unique forbidden transition is approximated as the shape factor of a (L -1) th unique forbidden transition. Because the vast majority of non-unique forbidden transitions present in a reactor νe spectrum are of the first kind, those were approximated as allowed transitions. The origin and the validity of the ξ-approximation remain poorly documented in the literature. They are discussed in the following. Going back to eq. 16 and keeping only the dominant terms in the multipole expansion, the shape factor of a first non-unique forbidden transition can be expressed as [START_REF] Behrens | Electron radial wave functions and nuclear beta-decay[END_REF]:

C(Z, W e ) = k n (1 + aW e + µ 1 γ 1 b/W e + cW 2 e ), ( 21
)
where the parameters k n , a, b and c are linear combinations of form factor coefficients which are independent of the electron and neutrino momenta. These parameters are defined in Table III. Historically, the notation ξ = αZ/2R was first introduced in previous β-decay formalisms [79][80][81]. The ξ-approximation assumes that the Coulomb energy of the β particle at the nuclear surface is large compared to the total decay energy, i.e. 2ξ W 0 . In this case, terms proportional to (W 0 R) and (W e R) in the shape factor parameters of eq. 21 can be neglected, and terms proportional to (αZ) dominate. The quantities entering the calculation of k n , a, b and c (see Table III) then compare as:

|A 0 | or |C 0 | |RB 0 | |W 0 RC 1 | |RD 0 | |W 0 RE 0 | |W 0 RF 0 | |W 0 RG 0 |, (22)
in which case the shape factor (eq. 21) simplifies to C(Z, W e ) k n (1 + µ 1 γ 1 b/W e ). This shape factor still does not resemble an allowed shape factor at this stage. By comparing this simplified shape factor to the general shape factor as expressed in eq. 16, the second term in the former can be identified to the third term in the latter. This term is proportional to the (M K m K ) product, and also appears in allowed transitions with the same value b = -2R(A 0 B 0 + C 0 D 0 )/(A 2 0 + C 2 0 ). As the M K quantities are usually greater than m K by an order of magnitude, terms connected with m K can be neglected and the shape factor hence becomes independent of the electron energy W e . The ξ-approximation is most likely to be correctly fulfilled for heavy nuclei. For instance, 2ξ is already ∼10 MeV for 70 Br and ∼18 MeV for 241 Pu. Still, many transitions for which the 2ξ W 0 criteria is met can exhibit a clear and significant deviation of their corresponding shape factor with respect to the allowed shape. For instance, a systematic comparison of experimentally measured shape factors with a theoretical prediction using the ξ-approximation showed that the latter often failed to reproduce the former [START_REF] Mougeot | [END_REF]. In such cases, the nuclear structure of the initial and final states are most likely to be responsible for this failure. This can be due to the so-called cancellation effect, where the form factor coefficients in |A 0 | and |C 0 | can compensate each other, making these quantities to be of the same order of magnitude than the others in eq. 22. Another possible situation is a modification of the selection rules in strongly deformed nuclei, the β transition being then denoted as K-hindered [82], with K the projection of the total angular momentum on the symmetry axis.

The same reasoning than the one explained here for first non-unique forbidden transitions can be applied to any higher order forbidden non-unique transitions. If the ξ-approximation criterion is fulfilled, the contribution of the form factor coefficients becomes independent of the electron and neutrino energies and can be factored out of the shape factor. A L th non-unique forbidden transition can then be treated as the corresponding (L -1) th unique forbidden one. Following these theoretical considerations, a conservative criterion to ensure the ξ-approximation to succeed for the calculation of nonunique forbidden transitions could be 2ξ/W 0 > 100, as prescribed in [START_REF] Mougeot | [END_REF]. Checking the decay data of the relevant transitions for the calculation of reactor νe spectra, this criterion is unfortunately met for transitions with endpoint energies smaller than 0.4 MeV only. Furthermore, relaxing this criterion to 2ξ/W 0 > 10 increases the maximum endpoint energy of these transitions to 3 MeV, showing that the ξ-approximation will very likely fail for the vast majority of non-unique forbidden transitions. In order to account for these possible modeling errors, an uncertainty is applied to each of the ξ-approximated transitions. In the same fashion that what is done in section III C 1, this uncertainty is constructed using the difference between a non-unique forbidden spectrum and its ξ-approximated version. Because each of the non-unique forbidden transitions the ξ-approximation is applied to are by definition not known, a simple strategy is to use a reference non-unique spectrum to compute this difference. This reference spectrum is modeled using the 23 non-unique forbidden transitions computed using nuclear

Shape factor C(Z, We) = kn(1 + aWe + µ1γ1b/We + cW 2 e ) with kn = A 2 0 + C 2 0 -2µ1γ1R 2 C1D0 + 1 9 (W0R) 2 (E 2 0 + G 2 0 ) a = R 2C0C1 -2 9 (W0R)(E 2 0 + G 2 0 ) /kn b = -2R(A0B0 + C0D0)/kn c = R 2 C 2 1 + 1 9 (E 2 0 + G 2 0 + λ2F 2 0 + λ2G 2 0 ) /kn where A0 = A F (0) 000 -1 3 αZ A F (0) 011 (1, 1, 1, 1) -1 3 W0R A F (0) 011 B0 = -1 3 A F (0) 011 C0 = -V F (0) 101 -1 3 αZ 1 3 V F (0) 110 (1, 1, 1, 1) -1 3 W0R 1 3 V F (0) 110 -1 3 αZ 2 3 A F (0) 111 (1, 1, 1, 1) + 1 3 W0R 2 3 A F (0) 111 C1 = -2 3 2 3 A F (0) 111 D0 = -1 3 1 3 V F (0) 110 + 2 3 A F (0) 111 E0 = 2 3 V F (0) 110 + 1 3 A F (0) 111 F0 = 2 3 V F (0) 110 -1 3 A F (0) 111 G0 = -A F (0) 211 TABLE III.
Theoretical expression of the shape factor for a first forbidden non-unique transition, as re-written in [START_REF] Behrens | Electron radial wave functions and nuclear beta-decay[END_REF]. Only dominant terms in the multipole expansion have been considered.

structure calculations (see Table II). Figure 4 shows the deviation of their shape factors with respect to their corresponding ξ-approximated version in the normalized νe energy representation E ν /E 0 , where E 0 is the endpoint energy. Because all of these transitions are first nonunique forbidden, these deviations are also deviations from an allowed shape. No general and systematic trend can be identified, making it difficult to describe them all with a unique and simple parametrisation. As such, all of these deviations are crudely approximated by a same linear function. This linear function is chosen such that it reasonably approximates the 93 Rb shape factor, which exhibits the largest deviation among all of these transitions. As illustrated in Figure 4, its amplitude reaches 150% at 0 MeV and vanishes at the endpoint energy. Constructed this way, this linear function bands the deviation of most of the calculated non-unique shape factors with their respective ξ-approximated version. It should hence conservatively describe the modeling errors caused by the ξ-approximation to most of the non-unique forbidden transitions still lacking nuclear structure calculations in the present computation of a reactor νe spectrum. The reference spectrum for each non-unique forbidden branch modeled with the ξ-approximation is then calculated multiplying this linear function by an allowed spectrum, and is used in the following way to construct the covariance matrix associated to a ξ-approximated transition. For each energy bin, the previously described refer-ence spectrum and the ξ-approximated transition spectrum are used to define the limits of a uniform distribution. This uniform distribution is conservatively symmetrized to take into account that non-unique forbidden shape factors could either exhibit positively or negatively slopped deviations. The standard deviation of the resulting uniform distribution is then used to generate a covariance matrix. Furthermore, this covariance matrix is constructed such that energy bins showing equal (resp. opposite) sign in the difference between the reference and ξ-approximated spectra are considered as fully correlated (resp. anti-correlated).

At the level of a branch, the ξ-approximation leads to a typical 30% uncertainty on the corresponding IBD yield. The contribution of ξ-approximated branches to isotopic IBD yields is at the level of 10-15% (see Table I) and their corresponding uncertainty shrinks to ∼0.4% (see Table VIII). The associated fractional uncertainty obtained for the 235 U fission spectrum is displayed in Figure 6, and shows that the uncertainty related to the ξ-approximation mostly dominates the total modeling uncertainty budget at low energies. No correlation among the ξ-approximated transitions has been considered in the previously quoted uncertainties. The impact of fully correlating these transitions has nevertheless been checked. It was found to be significant, increasing the uncertainty at the isotopic IBD yield level to ∼3%. Yet, and for the same reason than the one discussed in the previ-FIG. 4. Ratio of νe spectra of the 23 non-unique forbidden transitions listed in Table II over their associated allowed shape obtained under the ξ-approximation. The x-axis is expressed in normalized νe kinetic energy Eν /E0, where E0 is the transition endpoint energy (see eq. 5). The solid red curve highlights the ratio obtained for the 5/2 -→ 5/2 + transition of 93 Rb. The linear shape factor used to model the uncertainty resulting from the application of the ξ-approximation to this transition is represented by the dashed red curve (see text for more details). The application of the ξ-approximation to this transition results in a 1σ uncertainty band corresponding to the red area centered on unity.

ous section, the present work assumes no transition-totransition correlations. Finally, the robustness of this uncertainty modeling has been checked by computing the 235 U fission spectrum with (i) all non-unique forbidden transitions as treated by the ξ-approximation and (ii) including the 23 most contributing non-unique forbidden transitions as computed in section III C 1, the rest of the non-unique forbidden transitions remaining ξ-approximated. Figure 5 depicts the ratios of these two spectra both in the β and the νe energy representations. The application of the ξ-approximation to these 23 nonunique forbidden branches results in a 1σ uncertainty band (represented in red) which mostly covers the observed differences in both cases. Moreover, the 1.4% decrease in the 235 U IBD yield when going from case (i) to case (ii) is well covered by the 1.8% uncertainty resulting from the application of the ξ-approximation to all non-unique transitions. Similar decreases have also been obtained for the 238 U, 239 Pu and 241 Pu isotopic IBD yields, which equally amount to 1.1%. They are also well covered by the previously described uncertainty resulting from the ξ-approximation.

D. Radiative corrections

Radiative corrections come from non-static Coulomb processes arising at higher perturbation levels such as virtual photon exchange and inner brehmsstrahlung, in which one or multiple photons can be emitted in the final state. These effects can be important in the calculation of a β spectrum, especially for transitions with endpoint energy much larger than the electron mass. Radiative corrections are usually separated between an inner and outer part. The latter being the only one to include energy-dependent terms, no inner radiative correction is then considered in this work. Only outer radiative corrections of order α are taken into account. The higher order terms in α m Z n (m > n) are much smaller and they can be neglected to a first approximation [83]. The O(α) outer radiative correction for the electron can be found in [84]. Although the neutrino is insensitive to the nucleus Coulomb field, virtual photon exchange and energy conservation in inner brehmsstrahlung processes can still indirectly change its kinetic energy. The total O(α) radiative correction for the neutrino has been explicitly calculated in [85], and somewhat differs from the electron case because of a different inner brehmsstrahlung contribution. In the calculation of a νe spectrum at the βbranch level, the electron correction must therefore be replaced by the neutrino correction. Radiative corrections only depend on the transition endpoint energy. Electron radiative corrections typically range from less than 1% for transitions with endpoint energies close to the IBD threshold up to ∼10% for transitions with endpoint energies reaching 10 MeV. On the other hand, neutrino radiative corrections barely reach 1% for any transition with endpoint energy below 10 MeV. Once propagated to the calculation of a full reactor νe spectrum (see eq. 1), they total a 0.5% shape correction and a 0.2% correction on the corresponding IBD yield since most of the β-decay transitions exhibit low endpoint energies. Additionally, transitions with endpoint energies above the IBD threshold display a large range of endpoint energies up to 10 MeV and thus transition spectra average out in a reactor νe spectrum.

The radiative corrections to the calculation of an electron and νe spectrum at the β-branch level have been derived for allowed transitions only [84,85]. The procedure followed in this work and also in other previous summation models [17,18,41] is to apply these same corrections to any other type of β-decay transition. The validity of this procedure is not well-known. As such, an uncertainty associated to the outer radiative correction modeling is conservatively built similarly to what is done for the accurate modeling of non-unique forbidden transitions in section III C. The associated covariance matrix is here constructed using the difference between an actinide fission spectrum (see eq. 3) calculated with and without the application of outer radiative corrections. When propagated to the computation of a full reactor νe spectrum, this source of uncertainty is treated as fully correlated be- tween fission spectra since the corrections always behave in a similar fashion. Using this approach, the radiative correction uncertainty is ∼0.1% when propagated to the IBD yield of each actinide (see Table VIII). The associated fractional uncertainty obtained in the case of the 235 U fission spectrum is displayed in Figure 6, showing that this uncertainty source is negligible with respect to the other β-decay modeling uncertainties.

E. Weak magnetism correction

Weak magnetism (WM) refers to the dominant contribution of a class of additional induced-nuclear currents appearing in the vector part of the β-decay Hamiltonian when taking into account the finite-size and the internal structure of the nucleons. In the Behrens and Bühring formalism [START_REF] Behrens | Electron radial wave functions and nuclear beta-decay[END_REF], a weak magnetism correction to the calculation of a β/ν e spectrum at the branch level in principle only applies to Gamow-Teller allowed and nonunique forbidden transitions, as forbidden transitions of the unique type solely depend on axial nuclear form factors. Equivalent WM corrections for allowed Gamow-Teller transitions were consecutively derived in [50,86], and were extensively applied to any type of transitions in the past summation calculations of reactor νe spectra [17,18,41]. The present work considers the weak magnetism correction derived in [37] for allowed Gamow-Teller transitions. Because no clear prescription about weak magnetism in non-unique transitions yet exists in the literature, the correction for allowed branches is indistinctively applied to the corresponding spectra. The WM correction exhibits a similar dependency to the transition endpoint energy than the radiative corrections (see section III D). The magnitude of the correction is O(1%) (resp. ∼4%) for transitions with E 0 close to the IBD threshold (resp. 10 MeV). Once propagated to the calculation of an actinide fission νe spectrum, weak magnetism typically introduces up to a +0.1% correction below 2 MeV and a linearly decreasing correction reaching -1.5% at 8 MeV [START_REF] Perisse | Modeling of reactor antineutrino spectra[END_REF]. In the same fashion than for radiative corrections, an uncertainty associated to the WM correction is built up by comparing actinide fission νe spectra calculated with or without this correction, and constructing the associated covariance matrix such that the total rate is conserved. This source of uncertainty is fully correlated between each actinide fission spectra since their respective WM correction are similar. This prescription results in a relative uncertainty of about ∼0.3% for each isotopic IBD yield (see Table VIII). The associated fractional uncertainty obtained in the case of the 235 U fission spectrum is displayed in Figure 6. Similarly to the radiative corrections, the weak magnetism correction plays a minor role in the combined uncertainty of a reactor νe spectrum.

IV. EVALUATED NUCLEAR DATA

The summation calculation of a reactor νe spectrum, as depicted in eq. 1, eq. 4 and eq. 5 requires a large set of various nuclear data. In the present work, a new database parsing recent evaluated nuclear databases available online has been built. The construction of this database starts from the fission yield data (section IV A), which defines the list and probability of occurence of all possible products for a fissioning system irradiated by neutrons, such as one of the four major actinides 235 U, 239 Pu, 241 Pu and 238 U present in a reactor core. After the identification of all the β emitters, a database of evaluated decay information is constructed using the Evaluated Nuclear Structure Data Files (ENSDF) [START_REF] Nndc | From ENSDF database as of June 2020[END_REF] as described in section IV B. All modern evaluated nuclear decay databases are known to suffer from the so-called Pandemonium effect, which can significantly impact the calculation of reactor νe fluxes and spectra. Section IV C describes the correction of several important fission products decay information using recent Total Absorption Gamma-ray Spectrometry (TAGS) data available first, and then integral β spectrum measurements. The incompleteness of the evaluated nuclear decay databases leaves many fission products with no decay information. A new modeling of their contribution to the calculation of an actinide fission spectrum is proposed in section IV D. All along these sections, a particular attention is paid to precisely describe the usage and the uncertainty treatment of all the evaluated nuclear data so as to detail and understand how they contribute to the combined uncertainty budget of a fission actinide νe spectrum and flux.

A. Fission yield information

In eq. 1, the activity of each fission product is necessary to properly weight their corresponding contribution into the calculation of a reactor β/ν e spectrum. The total activity of each fission product can be computed over a reactor cycle using independent fission yields together with a reactor evolution code able to estimate a core inventory at irradiation time t. An independent fission yield I k p refers to the probability that a particular nuclide p will be produced directly from the fission of the k th actinide. The activity A k p (t) of such a fission product can be estimated using a cumulative fission yield Y k p (t) such that:

A k p (t) = f k (t) × Y k p (t), ( 23
)
where f k (t) is the fission rate of the k th actinide at irradiation time t. The cumulative yield Y k p (t) can here be understood as the total probability that the nuclide p is present after a time t, meaning it either is due to direct production from a fission event or comes from the decay of a parent fission product.

The fission yield (FY) information required to the computation of a reactor νe spectrum can be extracted from several nuclear data libraries, the most popular being the Joint Evaluated Fission and Fusion Files (JEFF, EU) [87], the Evaluated Nuclear Data Files (ENDF, USA) [88] and the Japanese Evaluated Nuclear Data Library (JENDL, Japan) [89]. For a given fissionning system, both independent and cumulative FYs each estimated for a set of three neutron energies (thermal, epithermal and fast) are provided. The cumulative FYs are those estimated at infinite irradiation time. For a fissioning system k, they are recursively computed as:

Y k,∞ p = lim t→∞ Y k p (t) = I k p + m b mp Y k m , (24)
where b mp is the probability that the parent nuclide m decays to the daughter nuclide p. The FY evaluation process in these libraries pretty much follows the same methodology, where basically a set of selected experimental data are combined with semi-empirical fission models to coherently assess all the independent and cumulative products yields for a given fissioning system. Moreover, the FY evaluated data can nearly be identical from one library to another. A good example is JENDL, which until very recently used to include the FY data from the ENDF evaluation [90].

In the present work, and unless otherwise indicated, thermal (resp. fast) cumulative fission yields at infinite irradiation time are used to compute the 235 U, 239 Pu and 241 Pu (resp. 238 U) actinide fission spectra. In other words, neither off-equilibrium corrections nor the impact of different fission neutron energies are taken into account. Assessing the impact of both these effects and including them into the present calculations is beyond the scope of this article. However, off-equilibrium effects have already been studied in the past [17,91]. They typically yield a small 2% negative correction to the 1.8-3 MeV portion of a νe spectrum emitted by commercial pressurized water reactors burning fuel over a typical 12to 18-month cycle. The variation of the fission product yields with neutron energies has recently been studied by simulating and comparing low-enriched and highlyenriched reactor core designs [92]. The corresponding actinide fission νe spectra showed small O(1%) changes below ∼5 MeV and up to ∼10% beyond. Their corresponding IBD yields were found to change by less than 1%. Both these corrections then fall well within the final uncertainty budget of an actinide fission spectrum calculation (see figure 12). They also give an IBD yield correction negligible with respect to its typical uncertainty (see Table VIII).

The latest FY evaluations from the JEFF, ENDF and JENDL librairies are here considered. The FY data are extracted respectively from the JEFF-3.3 [START_REF] Plompen | [END_REF], the ENDF/B-VIII.0 [93] and the JENDL-5 [94] releases. Following the prescription of [95], some of the FY data of the ENDF/B-VIII.0 release have been corrected from an erroneous evaluation, mostly leading to anomalously large uncertainties. As opposed to past releases, the JENDL-5 release includes now for the first time FY information which do not rely anymore on the ENDF evaluation [90]. Although sharing common experimental data and similarities in the underlying models to describe the fission process, the FY evaluation of these three database releases were conducted independently, making then worth a comparison for the calculation of reactor νe spectra. Any of these libraries ever included information about the correlations between the fission product yield information, whether these correlations come from the experimental data or the evaluation method. The complexity of determining FY correlations is a worldwide recognized problem and has been investigated over the last decade by many groups [96][97][98]. A set of matrices estimating the covariances of the independent and cumulative FY for the fission of the 235 U, 238 U, 239 Pu and 241 Pu actinides and exclusively sourcing from the evaluation process have been recently estimated in [99] for the JEFF-3.3 and ENDF/B-VIII.0 releases. Covariance matrices originating from the FY evaluation process are also included in the JENDL-5 library [90]. They are used here as an attempt to assess the impact of FY correlations among fission products.

Figure 7 (a) compares the νe spectra from the thermal fission of 235 U using cumulative FY from these 3 libraries. The most notable differences are driven by a limited number of fission fragments having important yields which can differ up to an order of magnitude. For instance, the decrease in the ENDF/JEFF ratio around 4 MeV is mostly due to 102 Tc. A list of the most relevant isotopes having notable and significant FY differences among these libraries is shown in Table IV. These differences, combined with a slightly different list of FPs in these three libraries, induce isotopic IBD yield variations at the level of ∼0.5% for 235 U, ∼2% for 239 Pu and 241 Pu, and ∼3% for 238 U. These IBD yield variations indicate a small tension between these three librairies when considering the uncertainty budget associated to FY as shown in Table VIII. Moreover, the different fission νe spectra are not always consistent with each other within their 1σ uncertainty bars, as seen in Figure 7 (a). The relative uncertainty on the 235 U νe spectrum induced by the cumulative FY uncertainty is illustrated in Figure 7 (b) for the three libraries. For each library, the FY uncertainties are propagated either by assuming the FY uncorrelated or by using the covariance matrices as described above. FY-to-FY correlations, induced by the evaluation process in the JEFF-3.3 and JENDL-5 libraries, decrease the 235 U fractional uncertainty in the 2-7 MeV region. A 1.1 → 0.8% decrease in the corresponding IBD yield uncertainty is observed. Similar effects are also found for the 238 U, 239 Pu and 241 Pu νe spectra and IBD yields. In the opposite way, the FY uncertainties and covariances evaluated in the ENDF/B.VIII.0 library increase both the fractional uncertainty budget and the IBD yield (2.2 → 4.0%) of the 235 U spectrum. This behavior has not been further investigated at the present stage. It should however be noted that the ENDF/B.VIII.0 FY evaluation dates back from the 2000s, and then may well be outdated with respect to the more recent JEFF-3.3 and JENDL-5 evaluations.

To ease any future comparison, the FY evaluation from the JEFF-3.3 library is here chosen, as it is the the most commonly used among the past and current summation predictions [17,41,42]. Furthermore, the following summation calculations conservatively keep uncorrelated FY information among the fission fragments because (i) the previously mentioned set of covariance matrices are incomplete as they do not include any correlations sourcing from the experimental data the FY evaluation are based on and (ii) the use of these matrices results in a decrease of the uncertainties both on the isotopic IBD yields and the actinide fission νe spectra. Beyond the existing correlations between individual fission product yields, correlations between different fissioning isotopes such as 235 U, 238 U, 239 Pu and 241 Pu in the present study, are also expected [97]. No evaluation could however be found neither among the existing libraries nor in the literature. They are then disregarded in the present work.

B. Nuclear structure and decay data

As shown by eq. 4, 5 and 11, the nuclear structure and decay information of all the fission products making FIG. 7. Impact of the evaluated fission yield libraries on the computation of the νe spectrum from the thermal fission of 235 U. (a) Ratio of the νe spectra using the ENDF/B-VIII.0 and JENDL-5 libraries over the one using the JEFF-3.3 library. Uncertainties are propagated using correlations originating from the fission yield evaluation process (see text for more details). (b) Systematic uncertainty on the 235 U νe spectrum induced by the FY uncertainties for these three libraries. The dashed lines correspond to an uncertainty propagation procedure using correlations originating from the FY evaluation process, as opposed to the solid lines which assume the FY information to be uncorrelated. IV. List of fragments having relevant cumulative fission yield differences for 235 U among the JEFF-3.3, the ENDF/B-VIII.0 and JENDL-5 libraries. These fission fragments are mainly responsible for the differences seen in the associated ratios in Figure 7 (a).

up a reactor νe spectrum are necessary. In the present work, they are mostly extracted from the 2020, June 29 th release of ENSDF [START_REF] Nndc | From ENSDF database as of June 2020[END_REF], which includes nuclear structure and decay data for over 3000 nuclides. An ENSDF evaluation can be found for approximately 70% of the fission products entering the composition of a reactor spectrum, which corresponds to more than 600 isotopes and several thousands of transitions. The nuclear level properties and β transition information are extracted using the ENSDF++ program [100]. The following sections present the different types of evaluated nuclear and decay data necessary for the computation of a fission product β/ν e spectrum, and how their associated uncertainties are treated. Most importantly, many fission products have missing or incomplete decay information. When necessary, the procedure used to circumvent this lack of data is also described.

Branching ratio and β -intensity

Branching ratios (BR) and β -intensities along with their respective uncertainties are read from ENSDF. In the computation of a fission fragment spectrum, the sum of BRs among all listed transitions is normalized to its β - intensity. Corresponding uncertainties are propagated using a Monte Carlo method. All transitions having complete BR information are treated as Gaussian distributed random variables. The same goes for transitions missing a BR uncertainty, which are assigned a relative 10% uncertainty. This choice corresponds to the median value of the relative BR uncertainty distribution among all known transitions in ENSDF. Finally, transitions having only an upper (resp. a lower limit) as information about their corresponding BR are treated as random variables following an upper bounded (resp. lower bounded) uniform distribution. The total uncertainty of an isotope β/ν e flux originating from the BR information must equal the reported β -intensity uncertainty. Because any information on BR correlations are available in ENSDF, artificial correlations are introduced to meet this constraint. During this step, the phase space of possible correlations is randomly probed to pick the set of correlation parameters maximizing the corresponding isotope IBD yield uncertainty. Further details about the BR and β -uncertainty propagation method are given in [START_REF] Perisse | Modeling of reactor antineutrino spectra[END_REF]. Figure 12 illustrates the contribution of this uncertainty source to the fractional uncertainty of the 235 U νe fission spectrum, showing that it is among the least important with a O(1%) contribution below 7 MeV. These uncertainties combine to a ∼0.4% uncertainty once propagated to the calculation of the corresponding IBD yield.

Endpoint energy, nuclear level spin and parity information

The endpoint energy of a transition is estimated using eq. 5, and needs the total β decay energy Q β , the parent isotope isomeric state level E IS and the daughter isotope level E lvl information. Fission product Q β information are extracted from the 2020 release of the atomic mass evaluation (AME-2020) database [101], as the proposed evaluation is based on a least square analysis using multiple experimental data, and therefore seems more robust than a single measurement. As a result, the corresponding uncertainty is slightly reduced compared to the analytical Q β calculation based on the mass difference. The metastable parent nucleus energy level E IS and the corresponding uncertainty are taken from ENSDF. If missing, this information is then retrieved from the 2020 release of the NUBASE database [102]. Finally, the different daughter β-feeding state energies along with their uncertainty are extracted from ENSDF. The transition endpoint energy uncertainty is then estimated quadratically summing the Q β , E IS and E lvl uncertainties. Spin and parity of the parent and daughter nuclear levels are necessary input information to determine the type of a transition. These information are also extracted from the ENSDF database. Should they may be missing or not fully determined, the following choices are made for the computation of a transition spectrum. When a transition exhibits multiple spin/parity combinations, the corresponding β/ν e spectrum is computed as the spectrum average of all associated forbiddeness degree (FD) unique branches. If the spin of either the parent or daughter nu-cleus is missing, the transition spectrum is defaulty computed as an allowed one. The uncertainties sourcing from endpoint energy, nuclear level spin and parity information are simultaneously propagated using a Monte Carlo method in the computation of a fission fragment spectrum. Therefore, they are combined all together into a single covariance matrix. Endpoint energy is treated as a Gaussian distributed random variable, and is constrained to yield only positive values. Because the endpoint energy uncertainty of a transition is usually dominated by the fission fragment Q β uncertainty, and also because the Q β value is used to compute the endpoint energies of all fission fragment transitions, endpoint energy information is considered to be fully correlated among those transitions. After an endpoint energy is sampled, incompleteness in the spin/parity information of a transition is then considered. When many spin/parity combinations are possible, all resulting unique FD are evenly sampled. If a transition misses a spin information, the transition type is randomly sampled between allowed, 1 st , 2 nd and 3 rd unique forbidden transitions. Similarly to the BR and β -intensity information, uncertainties associated to endpoint energy, nuclear level spin and parity information are found to negligibly contribute to the total uncertainty budget of an actinide fission β/ν e spectrum (see Figure 12). They typically induce a ∼0.1% uncertainty on the 235 U, 238 U, 239 Pu and 241 Pu isotopic IBD yields (see Table VIII).

C. Correction for the Pandemonium effect

The decay scheme of a parent radionuclide is usually inferred by measuring the intensity and energy of the γrays emitted in the deexcitation cascade of the daughter nucleus. In the past decades, experimental apparatuses widely used High Purity Germanium detectors (HPGe) for the detection and measurement of these γ rays in coincidence with the β particles, because of their excellent energy reconstruction performances. However, these devices mostly suffered both from limitations and from an incorrect characterization of their detection efficiency at high energies, often leading to an underestimate of the daughter nucleus level density at high excitation energies and hence biasing the parent nuclide decay scheme. This effect, first pointed out in [103], is called the Pandemonium effect and is known to be widely present in modern evaluated nuclear databases. At the level of a β decaying isotope, the Pandemonium effect underestimates (resp. overestimates) the β and νe spectra at low energies (resp. high energies). Beyond the summation prediction of reactor νe spectra, having reliable decay data is also important for nuclear reactor operation and safety considerations [104]. Therefore, several experimental efforts are conducted to correct the nuclear structure information and the decay scheme of the most important nuclides known to be Pandemonium-affected. The following sections discuss the different sources of Pandemonium-corrected decay information, the extraction of the relevant information and how they are applied to correct the library of nuclear structure and decay information used in the present actinide fission νe spectrum calculations.

Total Absorption Gamma-ray Spectrometry data

The most reliable Pandemonium-free sources of decay data come from Total Absorption Gamma ray Spectroscopy (TAGS) measurements. The TAGS experimental technique generally uses an arrangement of highefficiency γ-ray detectors (typically NaI or BaF 2 scintillating crystals) with a nearly 4π coverage able to fully reconstruct the γ cascade following a nuclear β decay. As opposed to HPGe devices, such γ detectors have a modest energy resolution, thus requiring the use of deconvolution techniques to properly assess the decay scheme of the parent nucleus. In this work, and whenever it is possible, TAGS data are prioritized over the ENSDF data of any known Pandemonium-affected nuclide. A first campaign of TAGS measurements was conducted in 1997 at the INEL ISOL facility using NaI(Tl) scintillation detectors by the group of Greenwood et al. [105]. Using these data, the decay schemes of 49 short-lived fission products extracted from ENSDF were corrected.

More recently, the emerging and increasingly pressing needs of reliable decay information for the predictions of νe fluxes and spectra emitted at nuclear reactors further accelerated the experimental efforts to measure (or remeasure) a selection of Pandemonium-affected radionuclides with the TAGS technique [106,107]. In this work, the decay schemes of 45 radionuclides have been then retrieved and corrected when necessary using these recent TAGS data. They are listed in Table V. Among these, 6 radionuclides ( 89 Rb, 90 Rb, 90 Rb m , 91 Rb, 94 Sr, 140 Cs) were already measured by the group of Greenwood et al., which the more recent TAGS data have then been prioritized over. Moreover, the 137 Xe ENSDF data have been validated by the TAGS measurement reported in [108]. Finally, the ENSDF decay data of 99 Zr and 87 Kr remained uncorrected in the present work since their corresponding TAGS data were either incomplete or not available in a usable format. For information, 99 Zr ( 87 Kr) represents 1.0%, 0.7%, 1.3% and 0.9% (0.3%, 0.1%, 0.2% and 0.1%) of respectively the 235 U, 238 U, 239 Pu and 241 Pu expected IBD yields.

All selected TAGS data sets have been here updated using the AME-2020 database [101] for the estimate of the Q β energies, and with the NUBASE-2020 database [102] for the determination of the metastable level energies and β -intensities. In addition, TAGS data alone often miss spin and parity information. Therefore, correspondences between the reported nuclear level energies in the TAGS data with those documented in ENSDF have been systematically searched for to preserve these information when using the decay scheme of a TAGSmeasured isotope. If no such correspondence could be identified, the corresponding transition was arbitrarily treated as allowed.

Table VI breaks down the calculation of each actinide fission νe spectrum and flux by source of decay information, especially showing that ∼60% (resp. ∼45%) of the 235 U and 239 Pu (resp. 238 U and 241 Pu) IBD yield is calculated using TAGS-corrected nuclear decay information. In a complementary way, Figure 8 shows the contribution of the TAGS source of decay information to the calculation of the 235 U and 239 Pu fission νe spectra, showing that they amount to ∼60% of the spectrum in the 2-8 MeV energy range. The impact of correcting the ENSDF decay data from the Pandemonium effect is illustrated on Figure 9 (a) for the calculation of the 235 U νe spectrum. A very similar behavior is observed for the corresponding β spectrum. The use of the data from Greenwood et al. and from the recent TAGS measurements decreases (resp. increases) the νe spectrum above (resp. below) ∼2 MeV by a few percents. This trend directly stems from the fact that Pandemonium-affected fission fragments have their corresponding β/ν e spectra underestimated (resp. overestimated) at low (resp. high) energies. As a result, the corresponding 235 U IBD yield decreases by 2.4% after including the TAGS data from Greenwood et al., and by another 5.4% after adding the most recent TAGS data as listed in Table V. The Pandemonium effect impacts the other actinide fission β/ν e spectrum and flux in a very similar way. Including and using the TAGS data as described previously result in a ∼6-8% decrease of their respective IBD yield.

Integral β spectrum measurements

Besides TAGS data, another valuable set of Pandemonium-free data comes from integral β spectrum measurements. Following the work of [17], measurements of the continuous β and γ-ray spectra emitted after the decay of 111 fission products by the group of Tengblad et al. [127,128] were here considered to further correct nuclear decay data extracted from ENSDF. These 111 fission products are mostly short-lived radionuclides with large Q β energies, and make up for ∼90% of a typical nuclear reactor νe flux above 6 MeV [127]. Among these 111 measured β spectra, 44 were found to be con- sistent with calculations using ENSDF decay data. The corresponding integral β spectrum data were hence disregarded. Among the remaining 67 isotopes, 23 isotopes have also recently been measured with the TAGS method as described in the previous section. As opposed to the TAGS technique, integral β measurements only give access to the full β spectrum of a radionuclide without any information about the underlying β decay scheme. Therefore, the integral β spectrum data set associated to these 23 radionuclides has also been disregarded. The integral β spectrum measurements of 44 isotopes finally remain after this selection. They are listed in Table VII. They have then been used to model their νe contribution instead of directly using the corresponding ENSDF extracted decay information.

Isotope 80 Ga, 81 Ga, 82 Ga, 83 Ge, 79 As, 81 As, 82 As, 83 As, 85 As, 86 As, 83 Se, 83 Se*, 89 Br, 90 Br, 87 Kr, 96 Rb, 97 Sr, 97 Y*, 98 Y, 99 Y, 99 Nb, 130 Sn, 130 Sn*, 131 Sn, 133 Sn, 131 Sb, 133 Sb, 134 Sb, 135 Sb, 136 Sb, 137 Sb, 136 Te, 137 Te, 135 The νe counterpart of each of these 44 fission products has been modeled by applying the so-called conversion method, in which the associated β spectrum is adjusted by a set of 2 to 6 virtual allowed branches, each with a branching ratio and an endpoint energy as free parameters [17]. The corresponding νe spectrum is then obtained by applying energy conservation to the adjusted virtual β branches, i.e. by substituting the β particle kinetic energy E with the νe energy E ν = E 0 -E, where E 0 is the adjusted endpoint energy. Applying this procedure, these 44 fission fragments are found to make a 10-17% contribution to each isotopic IBD yield (see Table VI) and a ∼10% contribution to an νe actinide fission spectrum (see Figure 8). Although the conversion procedure reproduces each isotope experimental β spectrum to less than a percent over the whole energy range, nothing ensures that the associated νe spectrum is accurately described by the set of adjusted virtual branches. As an example, Figure 10 depicts the ratio of νe spectra computed using the conversion procedure on integral β spectrum measurements over the same spectra using TAGS data for the previously mentioned 23 fission products both sharing these sources of data. Local excesses exceeding ∼50% are clearly visible. An uncertainty accounting for a potential bias originating from the conversion procedure must therefore be derived and applied to the 44 fission product νe spectra converted from the integral β measurements. The construction of the associated covariance matrix uses the 23 fission products both sharing integral β spectrum and TAGS measurements, and proceeds as follow: (i) each of these 23 fission products is assigned an individual covariance matrix computed using the relative difference between spectra expressed in normalized νe energy E ν /Q β from converted integral β measurements and from using TAGS data (see Figure 10) (ii) a single covariance matrix is computed by averaging out these 23 individual covariance matrices. The obtained covariance matrix encompasses then the average relative difference between a νe spectrum derived from integral β measurement conversion and from using TAGS data. It is then applied to model the uncertainty of each of the 44 νe spectra converted from integral β measurements. As an example, the red band in Figure 10 shows the 1σ uncertainty derived for 145 La. Because the distribution of the relative difference between the νe spectrum converted from integral β spectrum measurements and the νe spectrum calculated using TAGS data does not show any particular pattern, the uncertainty applied to the converted integral β spectra is treated as uncorrelated among the 44 associated isotopes. Using this approach, the uncertainty derived on an individual isotope IBD yield is typically ∼20%. At the isotopic IBD yield level, the associated uncertainty totals ∼1.5% of the uncertainty budget, making it the second most important source of uncertainty (see Table VIII). Figure 12 displays the associated fractional uncertainty for the calculation of the 235 U fission νe spectrum. It shows that the conversion of the integral β spectrum measurements is the dominant source of uncertainties above 5 MeV.

Finally, the impact of incorporating the data from the integral β measurements of Tengblad et al. [127] in the calculation of the 235 U actinide fission νe fission spectrum is illustarted by the blue curve in Figure 9 (a). These data brings a significant >10% decrease to the 4 MeV portion of the spectrum, especially because as stated previously, they correspond to radionuclides with large Q β energies. They bring a slightly larger correction to the calculation of the uranium (4.9% and 6.2% decrease for 235 U and 238 U, respectively) than the plutonium (2.9% and 3.7% decrease for 239 Pu and 241 Pu, respectively) isotopic IBD yields.

Uncorrected decay information

Although the decay information of the most important fission products entering the computation of a reactor νe spectrum have been reassessed and corrected when necessary, a substantial amount of the remaining fission products may still be potentially affected and thus be uncorrected from the Pandemonium effect. Such radionuclides are usually identified by comparing the energy E lvl max of the highest recorded level of the daughter nucleus to the total energy Q β available for the β decay. If a signifi-FIG. 10. Miscalculation of νe spectra resulting from the conversion procedure applied to integral β spectrum measurements of a set of 24 fission fragments also having TAGS data (blue lines). The x-axis is expressed in normalized νe kinetic energy Eν /Q β . The red line is the ratio associated to 145 La. The red band represents the 1σ uncertainty modeled to cover calculation errors in the case of 145 La (see text for more details).

cant difference between E lvl max and Q β is observed, a radionuclide is then highly suspected to be Pandemoniumaffected. Using a similar criterion, 29 isotopes present in ENSDF have been identified and priorily selected by the Working Party on International Evaluation Co-operation of the Nuclear Energy Agency (WPEC-25) for a new measurement of their corresponding β decay scheme with the TAGS technique [106,107]. As shown by Table VI, these isotopes represent 10-12% of an isotopic IBD yield. For curiosity, this selection was extended using a very loose criterion, i.e. searching the full ENSDF database for fission products with complete decay schemes and having at least a 20% difference between E lvl max and Q β . This criterion selects 172, 160, 202 and 190 additional isotopes for 235 U, 238 U, 239 Pu and 241 Pu respectively, meaning that in a worst case scenario, 17 to 20% of an isotopic IBD yield is impacted by a residual Pandemonium effect. The contribution of these isotopes to the calculation of the 235 U and 239 Pu νe fission spectra are also depicted on Figure 8. The blue short-dashed (resp. long-dashed) lines shows the WPEC-25 (resp. extended) selection, showing that the residual Pandemonium effect would mostly impact the 5 MeV portion of a νe spectrum.

Because no other source than ENSDF is yet available, directly correcting these isotopes for the Pandemonium effect is impossible. The followed strategy is instead to apply an average correction, which is constructed by considering the average impact of the Pandemonium effect among a set of 81 radionuclides showing different decay information between their ENSDF and TAGS records. The average ratio of spectra computed with TAGS data over spectra computed with ENSDF data is taken as the correction. A covariance matrix is derived by considering the dispersion of the 81 ratios with respect to the average ratio. Using such approach, the IBD yield uncertainty for a fission fragment is ∼20% and safely covers the amplitude of the associated residual Pandemonium correction, which on average amounts to 18% for the 81 radionuclides discussed previously. At the actinide fission spectrum calculation stage, the uncertainty derived for the residual Pandemonium effect is treated as fully correlated among the relevant fission products. This choice is motivated by the fact that the Pandemonium effect systematically impacts the concerned fission fragment β/ν e spectra in the same way. The robustness of the residual Pandemonium correction construction has been checked by comparing the computation of a 235 U fission spectrum (i) correcting the 81 fission fragments with their TAGS data and (ii) correcting these same fission fragments with the residual Pandemonium correction instead. Results are displayed in Figure 9 (b) and demonstrate that the residual Pandemonium correction reproduces fairly well the (true) correction of the spectrum using TAGS data to better than 3% below 8 MeV. Furthermore, the difference is largely covered by its associated uncertainty, as depicted by the red band on that same figure. Because the rather loose selection discussed above might misidentify a fair amount of fission fragments having a wrong decay scheme, the residual Pandemonium correction and associated covariance matrix is here only applied to the 29 fission fragments selected by the WPEC-25 group. Applied this way, it decreases the 235 U, 238 U, 239 Pu and 241 Pu IBD yields by (2.3±2.6)%, (1.9±2.0)%, (2.0 ± 2.3)% and (1.8 ± 2.1)%, respectively. The amplitude of the residual Pandemonium correction to these isotopic IBD yields is then largely covered by the associated uncertainty. This is not surprising given the way the correction was constructed, but reassuring though given that this correction may not be perfectly suited to these 29 fission fragments. The residual Pandemonium correction uncertainty amounts to ∼2.5% when propagated to the calculation of the 235 U, 238 U, 239 Pu and 241 Pu IBD yields, making it by far dominant among all considered sources of uncertainty (see Table VIII). The corresponding fractional uncertainty obtained in the computation of the 235 U fission νe spectrum is also illustrated in Figure 12, showing that it prevails below 5 MeV.

D. Nuclides with no data

The crossing of the FP list extracted from the JEFF-3.3 evaluated FY database [START_REF] Plompen | [END_REF] with ENSDF [START_REF] Nndc | From ENSDF database as of June 2020[END_REF] leaves hundreds of emitters with neither nuclear structure information nor decay data. Table VI shows the estimated contribution of these nuclides with no data (NND) to the calculation of the 235 U, 238 U, 239 Pu and 241 Pu νe flux. The NND are usually short-lived and exhibit high Q β energies far above the IBD threshold. They can therefore play a non-negligible role in the calculation of a reactor νe spectrum, with a typical ∼5-10 % contribution to the IBD flux. The modeling of the NND component in past summation calculations was usually done either using the Gross Theory of β-decay [129,130] or the socalled Q β effective modeling [14,17] or a mix of both approaches [41,42]. None of these effective modelings is expected to accurately predict this contribution. Most importantly, no uncertainty treatment relative to this contribution has ever been proposed in the literature.

In the present work, another approach to compute the contribution of the NND has been attempted. In a nutshell, this approach computes and averages the summation spectra associated to several pools of fission fragments to estimate the NND contribution, and uses the dispersion among these spectra to build an associated covariance matrix. In a first step, each NND is attributed a pool of several nuclides having a complete decay scheme and chosen such that their total β-decay energy Q β ranges within ±10% of the NND Q β . The resulting pool size typically varies from a few tens of nuclides up to a hundred. Then, each pool has a reference β/ν e spectrum estimated by averaging out the individual nuclide spectra. A covariance matrix is also built by using the dispersion of the individual nuclide spectra around the pool average spectrum. Summing each of these NND emulated spectra with their corresponding cumulative FY and β -intensities then gives an estimate of their contribution to the calculation of an actinide fission νe spectrum. The uncertainty derived with this pool modeling on each of the NND IBD yield is typically at the level of ∼60%. No correlations among the calculated NND spectra are considered when propagating these uncertainties. At the fission actinide level, the NND contribution uncertainty to the IBD yield typically amounts to ∼1%, except for 235 U. As shown by Table VI, NND have a smaller contribution of only 3.2% to the 235 U IBD yield, thus leading to a total uncertainty of 0.6% (see Table VIII). The associated fractional uncertainty obtained in the case of the 235 U actinide fission νe spectrum is displayed in Figure 12, with a rising trend following the contribution of the NND to the νe flux (see left panel of Figure 8). The NND source of uncertainty starts to significantly contribute in the high-energy portion ( 6 MeV) of the spectrum.

The impact of the pool method for modeling the NND component in the summation calculation of the 235 U and 239 Pu fission νe spectra is investigated by comparing a modeling using the Gross Theory of allowed β-decay and a Q β effective modeling using 3 evenly distributed transitions having the same branching ratios. As illustrated on Figure 11, notable changes are visible in the high en- ergy part of the spectra above 6-7 MeV. In this energy regime, the Gross Theory (resp. Q β effective modeling) predicts smaller (resp. larger) fluxes than the previously described modeling of the NND component. The corresponding 1σ uncertainty, pictured by the blue band on Figure 11, pretty much (resp. hardly) covers the reported differences with respect to the Q β effective modeling (resp. the Gross Theory calculation). Both the Gross Theory and the Q β effective modeling bring a negative 1-3% change to the 238 U, 239 Pu and 241 Pu IBD yields with respect to the present calculation of the NND component. The corresponding changes to the 235 U IBD yield are however smaller than 1%. This can be traced to a smaller contribution of the NND for this actinide (see e.g. Table VI or Figure 8). As shown by Table VIII, the present pool modeling of the NND component induces an uncertainty of about 1% on each isotopic IBD yield.

V. RESULTS AND DISCUSSION

All along the previous sections, the summation calculation of the νe spectrum and flux following the fission of the 235 U, 238 U, 239 Pu and 241 Pu actinides has been methodically dissected and revised, keeping as one of the main objective the construction of a detailed and realistic uncertainty budget. The main results following this work are manyfold. They are discussed here after and challenged against state-of-the-art predictions and measurements.

Figure 6 and Figure 12 summarize how the uncertainty budget of the 235 U actinide fission νe spectrum respectively breaks down according to β-decay modeling and evaluated nuclear data sources of uncertainty. The main lesson learned is that the quality and the incompleteness of the evaluated nuclear databases dominate by far the uncertainties. This finding equally applies to the other actinides. Although this fact was already known for decades, it has here been rigorously assessed and quantified. In more details, the β decay scheme correction of the known Pandemonium-affected fission fragments significantly changes the spectrum, with an approximate -4%/MeV linear correction between 2 and 7 MeV (see Figure 9 (a)). This correction combines different sets of nuclear data, and exhibits an uncertainty budget which is largely dominated by the conversion of the integral β spectrum measurements of Tengblad et al. [127]. This uncertainty especially prevails the total uncertainty budget of the 235 U νe spectrum above ∼5 MeV. The very likely possibility that many other fission fragments are left with an incorrect decay scheme led to the ad hoc construction of the residual Pandemonium correction, as described in section IV C 3. Although the magnitude of this correction is far smaller than the correction arising from the combined sets of TAGS and integral β spectrum data (see Figure 9 (a)), its associated uncertainty has been constructed in a conservative way, and therefore dominates the final uncertainty budget of the present calculations below 5 MeV. Finally, the incompleteness of the fission fragment decay information, mostly represented here by the contribution of the nuclides without any decay data (NND, see section IV D), only impacts the present calculations above 6-8 MeV depending on the considered fission actinide (see e.g. Figure 8). Consequently, the associated uncertainty prevails the combined uncertainty budget in that same energy regime. Nuclear decay data wise, a sizeable 20-25% portion of an actinide fission spectrum is then still weakly modeled and points out the necessity of priorily measuring both the Tengblad et al. (see Table VII) and WPEC-25 (see section IV C 3) selections of fission fragments with the TAGS technique in order to improve the summation method. In the other hand, the NND component of the actinide fission spectra, which emerges in the very high energy portion of the spectrum, could well be better understood and tackled taking advantage of high energy measurements of reactor νe , such as the one recently released by the Daya Bay experiment [131]. Summation calculations are also sensitive to the details of the fission yield evaluation. Small tensions between the most recent evaluations have been observed (see e.g. Figure 7), and therefore calls for a careful review of the used nuclear input data. The lack of a complete covariance matrix estimate for any of the available evaluations especially makes the uncertainty budget associated to FY information incomplete. Furthermore, the yields populating either the ground state or the isomeric state of a fission product are also weakly evaluated information, when available [132]. They can significantly differ from one library to another, and can therefore considerably impact the calculations, especially in the high energy portion of the spectrum. Beta decay wise, pushing the branch modeling to a high level of refinement has been found to give a smaller impact (see section III). In particular, using advanced nuclear structure calcula-tion to realistically model the 23 main non-unique forbidden transitions, totaling all together ∼25% of the expected IBD yield, gives a modest 5% change to the νe fission spectra when compared to the widely used ξapproximation treatment (see Figure 5). This change is especially visible in the 4 MeV portion of the spectrum, where the contribution of these non-unique forbidden transitions takes over. Overall, the improved treatment of the non-unique transitions presented in this work gives a 0.5-5% uncertainty in the computation of the νe fission spectra above 2 MeV, which is about 2-3 times less than the uncertainties arising from the quality and incompleteness of the evaluated nuclear decay data discussed above.

These conclusions hold for the three other actinides, and can be further appreciated while computing their corresponding IBD yields. Table VIII presents the four isotopic IBD yields as obtained in the present work, along with a detailed breakdown of their respective uncertainty budget. They are found to have a total ∼3% relative uncertainty, which is slightly larger than the uncertainty budget corresponding to the latest conversion predictions [17,18]. Figure 13 (a) compares the 235 U and 239 Pu IBD yields against the Estienne-Fallot (EF) summation prediction [41,42], the Huber-Mueller (HM) [17,18] and the Kurchatov Institute (KI) [START_REF] Kopeikin | [END_REF] conversion predictions, as well as a selection of measurements achieved by the Daya Bay [133], NEOS-II [START_REF] Kim | The 30 th International Conference on Neutrino Physics and Astrophysics (NEUTRINO 2022)[END_REF], RENO [START_REF] Bak | [END_REF], Double Chooz [START_REF] De Kerret | Double Chooz Collaboration)[END_REF], Bugey-4 [51,136] and STEREO [137] IBD experiments. The EF and HM IBD yields are here evaluated by considering fission spectra estimated at 450 days of irradiation time. They are expected to negligibly differ with respect to a spectrum computed under full equilibrium conditions. The EF fission spectra are directly taken from [42], while off-equilibrium corrections at 450 days estimated in [17] are applied to the HM prediction. Moreover, HM fission spectra are completed below 1.825 MeV and above 8.125 MeV using the present summation calculations, which gives a ∼0.5% correction to the corresponding isotopic IBD yields. The comparison of the present IBD yield calculations to the EF prediction does not show any significant discrepancies within the reported uncertainties. Although these two summation predictions differ in e.g. the β decay modeling details or the treatment of the NND component, this result is not surprising because they use very similar sets of Pandemonium corrected nuclear decay data. As stated earlier, the correction of the evaluated nuclear decay data from the Pandemonium effect is the most impactful to the summation calculations of reactor νe spectra. A (7.5 ± 3.9)% difference in the 235 U IBD yield is observed with respect to the HM prediction, while 239 Pu IBD yields show a very good agreement. This difference is not highly significant, but is in line with the recent 235 U IBD yield measurements conducted by the LEU and HEU experiments cited above. Figure 13 (b) further illustrates this observation. Together with the latest STEREO measurement of the 235 U IBD yield [137], it expresses in the To investigate a step further the IBD yield differences with respect to the HM model, Figure 14 (a) compares the present summation calculations to the aggregate β spectra measured in the 1980s at the Institut Laue-Langevin (ILL) with the BILL magnetic spectrometer, for the thermal fission of 235 U and 239 Pu respectively. Off-equilibrium effects at low energies caused by the short neutron irradiation time of the ILL measurements are taken into account by using the FISPACT-II numerical code [138] together with IFYs from the JEFF-3.3 database [START_REF] Plompen | [END_REF]. The FP activities were calculated after a 12-h and a 36-h irradiation time for 235 U and 239 Pu, respectively. As a first approximation, their associated uncertainties were taken as those coming from the CFY evaluation. Significant discrepancies, linearly ranging from -15% up to +5% for 235 U and -10% up to +20% for 239 Pu can be observed between 1 and 5 MeV. Figure 14 (b) plots the S 5 /S 9 aggregate β spectrum ratios constructed from the ILL data and from the present summation calculations. Although a (7.0 ± 3.0)% mean offset is still present, these S 5 /S 9 ratios are in a closer agreement. The linear deviations observed in Figure 14 (a) then hint at a possible systematic effect both present in the 235 U and 239 Pu β spectra, and partially compensating when constructing the S 5 /S 9 ratio. At the present stage, nothing indicates whether this systematic effect comes from the ILL data or from the present calculations. For further comparison, data points from the KI (S 5 /S 9 ) measurement are also superimposed on Figure 14 (b). They exhibit a better agreement than the ILL data when compared to the present summation calculations. Still, significant deviations can be seen, especially at high energies.

The IBD yield comparison discussed above is further examined in Figure 15, which displays the ratio of both the 235 U and 239 Pu νe spectra as predicted with BESTI-OLE in the present work to either the HM or the EF model. Unsurprisingly, similar deviations than those observed in the direct comparison to the ILL data (see Figure 14 (a)) are also observed here in the comparison to the HM prediction. The 235 U and 239 Pu νe spectra from the present work and the EF prediction agree to within ± 5% in the energy range below ∼7 MeV. The large differences observed at higher energies is suspected to mostly come from a different treatment of the NND component. The EF summation prediction mostly uses the gross theory of β decay, which gives smaller νe fluxes at high energy than the pool modeling proposed in the present work (see Figure 11). Last but not least, using FY evaluation from the JEFF-3.1.1 library as in the original EF prediction has been found to significantly improve the agreement [137] is pictured by the orange vertical line. The light and dark shaded bands are respectively the 68% CL and 95% CL associated uncertainty. The solid line (dotted line) ellipses correspond to 95% CL (99% CL) contours from a global analysis using fuel evolution and absolute rate measurements at LEU and HEU reactors [30]. The dashed blue line corresponds to the (S 5 /S 9 ) aggregate beta spectrum ratio measured at the Kurchatov Institute. The blue triangle lying on this line corresponds to the KI prediction. The corresponding 68% and 95% CL contours are not displayed not to overload the figure. They are exactly the same than those of the HM prediction.

between the corresponding actinide fission νe spectra, es-pecially for the Plutonium isotopes. This last point again ) ratio measurement performed at the Kurchatov Institute (KI) [START_REF] Kopeikin | [END_REF]. The corresponding (S 5 /S 9 ) BESTIOLE /(S 5 /S 9 ) BILL,KI double ratios are displayed in the bottom panel.

demonstrates the importance of a robust evaluation of the fission fragment yields for more accurate summation calculations.

The last point of comparison focuses on the shape of the predicted fission νe spectra.

IBD spectrum measurements extracted from the combination of the PROSPECT data together with either the STEREO [139] (here denoted SP) or the Daya Bay [140] data (here denoted DBP) are here used as benchmarks. Figure 16 shows how the present summation calculations compare to the unfolded 235 U and 239 Pu νe experimental spectra. In details, both the experimental spectra and the summation prediction are area-normalized to perform a shape-only comparison. To account for residual effects in the experimental data unfolding process, the summation prediction is filtered using published smearing matrices as prescribed in the supplementary materials of [139,140]. Given the uncertainties both coming from the experimental measurements and those estimated in the present work, an overall good agreement between data and prediction is observed. A χ 2 value was computed and respectively gave χ 2 /ndf=19.1/21, χ 2 /ndf=25.8/23 and χ 2 /ndf=13.8/22 for the SP 235 U, DBP 235 U and 239 Pu data, demonstrating an overall good agreement with the present summation prediction in the 1.8-7.5 MeV energy range.

In light of the spectrum deviations recently reported in the 5-7 MeV energy regime with respect to the HM [139] and of the unfolded 235 U and 239 Pu spectra from the Daya Bay and PROSPECT joint measurement [140] to the summation prediction (BESTIOLE). The red curves correspond to the bestfit Gaussian distortions (see text for more details), whose respective parameters are displayed in Table IX.

model, a Gaussian distortion was searched for:

M(E ν ) = K • S k (E ν ) 1 + A e -(Eν -µ) 2 2σ 2 , ( 25
)
where K is a global normalisation parameter allowing for a shape-only comparison, and S k (E ν ) is the actinide fission spectrum as defined in eq. 3 and eq. 1. The results are summarized in Table IX. They are also superimposed (red solid lines) to the displayed ratios in Figure 16. As anticipated and shown by the χ 2 /ndf values reported in Table IX, none of the experimental reactor νe datasets significantly favors the Gaussian distortion hypothesis in this energy regime. The DBP 235 U measurement, which shows the largest deviation to the summation prediction, prefers the Gaussian distortion hypothesis at the 2.3σ level.

VI. CONCLUSIONS

In this article, the summation method was deeply improved using an advanced β-decay formalism together with recent evaluated nuclear data. For the first time, a complete uncertainty budget accounting for all known effects likely to impact the calculation of reactor νe spectra is proposed. The modeling of all β transitions known to contribute to a reactor νe spectrum has been greatly refined using the Behrens and Burhing formalism. A majority of these transitions are of the allowed and unique forbidden type. They are now computed accurately through the implementation of the main electro-magnetic corrections to the Fermi theory and an exact calculation of their corresponding shape factors. Furthermore, the treatment of the non-unique forbidden transitions, which were for long ago anticipated to play an important role, was realistically tackled using nuclear structure calculations. In particular, the ξ-approximation, which is abusively used to model these transitions, has been demonstrated to have a modest percent level impact on the prediction of isotopic IBD yields. The quality of the evaluated nuclear data was also extensively assessed. The fission yield evaluations available in the most recent JEFF-3.3, ENDF/B-VIII.0 and JENDL-5 libraries were used and compared, showing significant discrepancies in the resulting actinide fission νe spectra at 5 MeV energies. These discrepancies together with a lack of a complete set of covariance matrices correctly describing the correlations between the fission fragment yields make the uncertainty budget associated to the summation calculations incomplete in that respect. Furthermore, the reliability of the fission fragment decay data was scrutinized using the most up-to-date Pandemonium free measurements available in the literature and in online databases. Although many experimental efforts are being conducted to improve the quality and the reliability of those data, the present summation model remains potentially affected by a residual Pandemonium effect still not corrected for in the present day nuclear databases. An ad hoc correction was then constructed to take it into account. This effect currently dominates the associated uncertainty budget, making it necessary to further measure or remeasure the potentially impacted fission fragments with the TAGS technique. Finally, a new approach to estimate both the contribution and uncertainty of all known fission fragments left with no decay data has been developed and validated against the usual Q β effective and Gross Theory modeling. These fission fragments only impact the high energy portion of the actinide fission νe spectra, and play for now a secondary role. Following all these improvements, this newly revised summation model was then extensively compared against a selection of state-of-the-art predictions and measurements. General good agreement is achieved with measured IBD yields at LEU and HEU reactors, especially favoring the RAA to be mostly caused by an overestimate of the 235 U flux in the HM model. While significant discrepancies with the 235 U and 239 Pu aggregate β spectra measured at the ILL are observed, the present work shows closer agreement with the recently measured S 5 /S 9 aggregate spectrum ratio at the Kurchatov Institute. The summation calculations achieved 25) in the STEREO and PROSPECT (SP) 235 U and Daya Bay and PROSPECT (DBP) 235 U and 239 Pu unfolded νe spectra. The last two columns, denoted H0 and H1, report the χ 2 /ndf values obtained for the "no Gaussian distortion" and "Gaussian distortion" hypotheses, respectively. in this work also demonstrate good shape agreement with a set of recently measured νe spectra from the Daya Bay, PROSPECT and STEREO experiments. All of these comparison studies show that the present summation calculations pretty well describe the most recent IBD flux and spectrum measurements. They are however unable to reconcile them with the original ILL data the HM conversion prediction relies on, further casting doubts on the reliability of these data. The present work does however not allow to favour any particular scenario for understanding the origin of the RAA. Further improvements both in the future summation modeling and measurements of reactor νe flux and spectra are necessary to firmly conclude.

Finally, the uncertainty budget constructed in the present work had to address in a realistic and conservative way the many systematic effects arising from gaps in the modeling of several fission fragments. These deficiencies either come from knowingly non-valid approximations (e.g. ξ-approximation), unreliable information (e.g. Pandemonium effect) or even missing information. The uncertainties and/or corrections associated to these effects were all estimated following a similar strategy, which uses subsets of fission fragments having known and reliable information as proxies. Implicit in this strategy is the assumption that the uncertainties and/or corrections derived from using these proxies are fully representative of those that would apply to fragments actually having these gaps. Although the validity of this assumption could be questioned, the observed good agreement of the present summation calculations with the latest IBD flux and spectrum measurements indicates no significant bias in either the modeling of these corrections or the modeling of these uncertainties.

To conclude, the summation method made an important step toward becoming a reference tool for the prediction of reactor νe fluxes and spectra. In this respect, the calculated νe (β) fission spectrum data points for each of the four actinides are provided under the form of 25 keV binned histograms in the supplementary materials of this article. Moreover, a total covariance matrix for the major four νe (β) actinide fission spectra, including each spectrum covariance matrix as well as cross-covariance matrices, is provided. The νe spectra and the associated covariance matrices have been extended below the IBD energy threshold to 0 MeV and at high energies up to 12.5 MeV, both allowing for a proper comparison with the forthcoming CEvNS experiments and future IBD ex-periments sensitive to the highest reactor νe energies. Finally and for completeness purposes, the spectra and covariance matrix of the most relevant activation products created by the neutron irradiation of reactor fuel and structural materials, namely 239 U, 239 Np, 28 Al, 56 Mn, 6 He, 52 V, are also provided in the supplementary materials under the exact same format.

FIG. 2 .

 2 FIG. 2. Comparison of allowed and forbidden unique transitions in the Behrens and Bühring β decay formalism. (a) νe spectra computed with accurate λ k Coulomb functions for a set of fictitious transitions of different forbiddenness degree with Z = 46, A = 117, and E0 = 10 MeV. (b) Impact of accurate λ k Coulomb functions compared to the λ k = 1 approximation on the calculation of shape factors for this same set of transitions. Solid lines include the nuclear finite-size effect into the computation of Coulomb functions, while dashed lines additionally include the atomic screening effect (see section III B for more details).

FIG. 5 .

 5 FIG. 5. Impact of the ξ-approximation to the calculation of the 235 U fission νe spectrum. Panel (a) (resp. panel (b)) plots the ratio of the β (resp. νe) spectrum including the 23 non-unique forbidden branches from TableIIas obtained from the nuclear structure calculations (NSC) described in section III C 1, over the same spectrum where all non-unique transitions are treated under the ξ-approximation. The red area represents the 1σ uncertainty band resulting from the application of the ξ-approximation to these 23 non-unique forbidden branches. Black error bars represent the uncertainty derived from the nuclear structure calculations if applied to these same 23 transitions instead.

FIG. 6 .

 6 FIG. 6. Breakdown of fractional uncertainties associated to the modeling of β branches in the summation calculation of the 235 U fission νe spectrum. Uncertainties associated to the treatment of non-unique forbidden transitions are further broken down between the 23 branches computed with nuclear structure and described in Table II (blue dashed line) and the rest being ξ-approximated (blue dotted line).

  MeV] 102 Tc (4.29 ± 0.06) × 10 -2 (2.24 ± 0.36) × 10 -3 (4.29 ± 0.04) × 10 -2 4.53 102 Tc * (4.18 ± 1.51) × 10 -5 (4.29 ± 0.47) × 10 -2 (2.16 ± 0.78) × 10 -5 4.58 97 Y (2.16 ± 0.25) × 10 -2 (4.89 ± 1.12) × 10 -2 (1.90 ± 0.26) × 10 -2 6.82 97 Y * (2.52 ± 0.20) × 10 -2 0 (3.06 ± 0.27) × 10 -2 7.49 100 Nb (5.54 ± 0.10) × 10 -2 (3.11 ± 1.00) × 10 -2 (5.71 ± 0.13) × 10 -2 6.40 100 Nb * (6.18 ± 1.95) × 10 -3 (3.11 ± 1.00) × 10 -2 (4.87 ± 1.25) × 10 -3 6.71 96 Y (4.66 ± 0.15) × 10 -2 (6.00 ± 0.96) × 10 -2 (5.73 ± 0.14) × 10 -2 7.11 86 Ge (8.67 ± 3.01) × 10 -6 (6.29 ± 1.01) × 10 -3 (3.45 ± 1.20) × 10 -7 9.56 88 As (1.57 ± 0.61) × 10 -5 (1.24 ± 0.56) × 10 -3 (2.56 ± 0.99) × 10 -5 13.43 92 Rb (4.37 ± 0.19) × 10 -2 (4.82 ± 0.07) × -2 (5.01 ± 0.22) × 10 -2 8.09 TABLE

FIG. 8 .

 8 FIG. 8. Contributions of the different sources of decay information used to compute the 235 U and 239 Pu fission νe spectra. Contributions associated to solid lines add up to 100%. Contributions from isotopes possibly affected by a residual Pandemonium effect (dotted or dashed blue lines) belongs to the ENSDF data contribution (solid blue line). The RP WPEC-25 case corresponds to the contribution of 29 radionuclides identified by the WPEC-25 group potentially having a Pandemonium effect while the RP extended case additionally includes hundreds of isotopes based on a broader selection (see section IV C 3 for further details).

FIG. 9 .

 9 FIG. 9. Correction of the Pandemonium effect in the calculation of the 235 U νe fission spectrum. (a) Impact of different sources of Pandemonium-free decay data. The ratios of spectra computed with different nuclear decay data sources over a spectrum calculated using ENSDF extracted decay data only are displayed. The Greenwood, Recent TAGS, Tengblad and RP correction labels respectively correspond to TAGS data coming from [105], to TAGS-corrected radionuclides as listed in Table V, to integral β spectrum measurements from [127] and to the 29 isotopes which the residual Pandemonium correction has been applied to. The yellow band shows the 1σ uncertainty resulting from applying the residual Pandemonium correction. (b) Coherence of the residual Pandemonium correction, as derived in section IV C 3. The black curve compares a 235 U spectrum computed with TAGS-corrected data (Greenwood and Recent TAGS) to a spectrum computed using ENSDF extracted decay data only but applying the residual Pandemonium correction to these same isotopes instead of using their TAGS data. The red band is the 1σ uncertainty resulting from the residual Pandemonium correction.

FIG. 11 .

 11 FIG.11. Impact of different approaches to compute the contribution of the fission fragments with no decay information to the summation calculation of the 235 U and 239 Pu fission νe spectra. The comparison of the Gross Theory and Q β effective modeling is made against the pool modeling described in section IV D. The blue area centered on unity represents the 1σ uncertainty band derived for the pool modeling.

FIG. 12 .

 12 FIG. 12. Uncertainties in the summation calculation of the 235 U fission νe spectrum, broken down by source and type of evaluated nuclear data.

FIG. 13 .

 13 FIG. 13. Comparison of the IBD yields as obtained with BESTIOLE in the present work to a selection of state-of-the-art predictions and measurements. (a) Comparison of the isotopic IBD yields for 235 U, 239 Pu and the combination of 235 U, 238 U, 239 Pu and 241 Pu as measured at LEU commercial reactors. The shaded areas correspond to the 1σ uncertainty band estimated from the present summation calculations. The EF IBD yield predictions miss an uncertainty, because those were not evaluated in [41, 42]. (b) Comparison of IBD yields expressed relatively to the HM prediction in the ( 235 U, 239 Pu) plane. The HM conversion prediction is pictured by the blue cross. The green dot and red inverted pyramid respectively correspond to the present summation calculations from BESTIOLE and from the EF prediction. The dark (light) shades are the 68% CL (95% CL) contours for the BESTIOLE summation calculations (green) and the HM prediction (blue). The latest STEREO measurement of the 235 U IBD yield[137] is pictured by the orange vertical line. The light and dark shaded bands are respectively the 68% CL and 95% CL associated uncertainty. The solid line (dotted line) ellipses correspond to 95% CL (99% CL) contours from a global analysis using fuel evolution and absolute rate measurements at LEU and HEU reactors[30]. The dashed blue line corresponds to the (S 5 /S 9 ) aggregate beta spectrum ratio measured at the Kurchatov Institute. The blue triangle lying on this line corresponds to the KI prediction. The corresponding 68% and 95% CL contours are not displayed not to overload the figure. They are exactly the same than those of the HM prediction.

FIG. 14 .

 14 FIG. 14. Comparison of the BESTIOLE summation prediction to the ILL aggregate β spectrum data [5, 6]. Prediction uncertainties are represented by the blue band. Black error bars correspond to experimental uncertainties. (a) Ratios to the 235 U (top) and 239 Pu (bottom) ILL β spectra. (b)The top panel shows the (S 5 /S 9 ) β spectrum ratios both coming from the summation prediction (BESTIOLE), from the ILL data (BILL), and from the recent (S 5 /S 9 ) ratio measurement performed at the Kurchatov Institute (KI)[START_REF] Kopeikin | [END_REF]. The corresponding (S 5 /S 9 ) BESTIOLE /(S 5 /S 9 ) BILL,KI double ratios are displayed in the bottom panel.

FIG. 15 .

 15 FIG. 15. Comparison of the BESTIOLE summation calculations to the EF and HM 235 U and 239 Pu νe fission spectra. Prediction uncertainties from BESTIOLE are represented by the blue band. Black uncertainty bars correspond to HM uncertainties. No uncertainties are available for the EF prediction.

TABLE I .

 I Contributions of the different types of β decay transitions to the summation calculation of the νe fluxes (top) and IBD yields (bottom) from the fission of the 235 U, 238 U, 239 Pu and 241 Pu actinides. Computations are done using cumulative fission yields from JEFF-3.3

		235 U	238 U	239 Pu	241 Pu
	νe flux [%]
	Allowed	57.9	62.9	64.8	66.9
	1 st non-unique forbidden	26.1 (6.7) 22.8 (5.9) 22.3 (4.6) 21.1 (4.7)
	1 st unique forbidden	11.4	7.8	8.0	6.4
	Other non-unique forbidden 3.4	2.8	3.7	3.1
	Other unique forbidden	0.1	0.1	0.1	0.1
	Nuclides with no data	1.1	3.6	1.1	2.4
	IBD yield [%]
	Allowed	43.9	50.1	52.8	55.4
	1 st non-unique forbidden	41.9 (30.6) 33.4 (20.9) 33.2 (23.4) 30.0 (19.2)
	1 st unique forbidden	10.0	5.7	8.4	5.4
	Other non-unique forbidden 0.9	1.0	0.9	1.1
	Other unique forbidden	0.1	0.1	0.1	0.1
	Nuclides with no data	3.2	9.7	4.6	8.0

TABLE V .

 V 87 Br, 88 Br,94 Rb[112] 89 Kr, 89 Rb,90 Kr,90 Rb,90 Rb m , 92 Rb [113] 100 Nb, 100 Nb m , 102 Nb, 102 Nb m [120] 101 Nb, 105 Mo, 106 Tc, 107 Tc [121] List of fission fragments whose decay information are corrected following recent TAGS measurements.

	Isotope	Reference
	76 Ga	[109]
	84 Br, 85 Br	[110]
	86 Br, 91 Rb	[111]
	93 Rb, 139 Xe	[113]
	94 Kr	[114]
	94 Sr	[115]
	95 Rb, 137 I	[116]
	96 Y, 96 Y m	[117]
	98 Nb	[118]
	142 Cs	[119]
	100 Tc	[122]
	102 Tc, 104 Tc, 105 Tc	[123]
	101 Zr, 102 Zr, 109 Tc	[124]
	103 Mo, 103 Tc, 140 Cs	[125]
	103 Nb, 104 Nb m	[126]
	137 Xe	[108]

TABLE VIII .

 VIII 235 U 238 U 239 Pu 241 Pu IBD yield [10 -43 cm 2 /fission] Isotopic IBD yields and their corresponding uncertainties, broken down according to uncertainty sources from evaluated nuclear data and modeling of β branches.( 235 U, 239 Pu) plane the result of a global rate analysis using IBD data from HEU and LEU experiments relatively to the HM prediction[30], and shows that both are in tension with the latter. The IBD yields obtained through the summation method with BESTIOLE in the present work and in the EF prediction are also shown for comparison. They come in a very good agreement with these experimental results, and all together suggest that the RAA is mosty caused by an overestimate of the 235 U νe flux. The dashed blue line, which indicates how the 239 Pu/ 235 U IBD yield ratio would scale according to the recent KI measurement of the (S 5 /S9

			6.25 10.01 4.48 6.58
	Uncertainty [%]
		Endpoint + Spin-parity	0.1 0.1 0.1	0.1
		Branching ratio + β -intensity 0.4 0.3 0.4	0.3
	Data	Residual Pandemonium Tengblad	2.6 2.1 2.3 1.4 1.7 1.6	2.1 1.4
		Nuclides with no data	0.6 1.2 0.9	1.3
		Fission yield	1.1 1.2 1.3	1.3
		Weak magnetism	0.3 0.3 0.2	0.2
		Radiative correction	0.1 0.1 0.1	0.1
	Modeling	Nuclear structure calculation 0.3 0.2 0.2 ξ-approximation 0.3 0.5 0.3	0.2 0.4
		IBD cross-section	0.1 0.1 0.1	0.1
		Total	3.3 3.2 3.3	3.2
		) aggregate beta
	spectrum ratio, further supports this interpretation. In
	particular, the KI prediction assumes that the (5.4 ±
	0.2)% offset measured in this ratio is entirely caused by
	a wrong normalisation of the ILL original 235 U aggregate
	beta spectrum measurement. As shown by Figure 13
	(b), it exhibits a much better agreement with the IBD
	yield experimental data than the HM prediction does.
	Finally, the experimental uncertainties reported in the
	recent measurements of the 235 U and 239 Pu IBD yields
	are of the same order of magnitude than those obtained in
	the present summation prediction. As can be seen in Fig-
	ure 13, they cannot yet fully conclude about the origin
	of the RAA, demonstrating that further improvements
	both in future reactor νe flux predictions and measure-
	ments are necessary.

  235 U) 0.098 ± 0.050 5.74 ± 0.30 0.55 ± 0.20 19.1/21 15.1/18 DBP ( 235 U) 0.169 ± 0.061 6.24 ± 0.30 0.76 ± 0.29 25.8/23 16.1/20 DBP ( 239 Pu) 0.105 ± 0.088 5.60 ± 0.21 0.44 ± 0.10 13.8/22 12.3/19 TABLE IX. Best-fit parameters resulting from the search of a Gaussian distortion (see eq.

	A	µ [MeV]	σ [MeV]	H0	H1
	SP (

ACKNOWLEDGMENTS

The authors thank Dr. Chikako Ishizuka for kindly providing the fission yield evaluation of the latest Japanese Evaluated Nuclear Data Library release. Dr. Muriel Fallot along with the Valencia-Nantes collaboration for TAGS evaluation, and Dr. Ãleksandra Fija lkowska along with the ORNL-Warsaw group for TAGS evaluation are also graciously thanked for providing TAGS data. This work was done in the framework of the NENuFAR project, which was supported by the direction of the cross-disciplinary programs at Commissariat à l'énergie atomique et aux énergies alternative (CEA).