The antineutrino energy structure in reactor experiments

The recent observation of an energy structure in the reactor antineutrino spectrum is reviewed. The reactor experiments Daya Bay, Double Chooz and RENO have reported a consistent excess of antineutrinos deviating from the flux predictions, with a local significance of about 4$\sigma$ between 4 and 6 MeV of the positron energy spectrum. The possible causes of the structure are analyzed in this work, along with the different experimental approaches developed to identify its origin. Considering the available data and results from the three experiments, the most likely explanation concerns the reactor flux predictions and the associated uncertainties. Therefore, the different current models are described and compared. The possible sources of incompleteness or inaccuracy of such models are discussed, as well as the experimental data required to improve their precision.


Introduction
In the last two decades, several neutrino oscillation experiments [1] have demonstrated that neutrinos are massive particles. Thus, neutrinos have become a main probe to explore physics beyond the Standard Model. Within the three neutrino paradigm, the neutrino oscillation probability can be described by three mixing angles (θ 12 , θ 23 , θ 13 ), two independent mass square differences (∆m 2 21 , ∆m 2 31 ), and one phase δ CP responsible for the CP -violation in the leptonic sector. While the dominant oscillations driven by θ 12 and θ 23 have been measured by different experiments in the so-called solar and atmospheric sectors, the third mixing angle θ 13 remained unrevealed until very recently. The first direct indications of a non-zero value of this angle has come from the accelerator-based experiments MINOS [2] and T2K [3]. However, the current accelerator neutrino experiments cannot measure θ 13 independently of other oscillation parameters. Complementing the role of accelerator-based facilities, reactor neutrino experiments stand as the direct way to provide an accurate value of θ 13 . In a two flavors scheme and for short baselines (L∼2km), the survival probability of a reactor electron anti-neutrinoν e with energy E ν can be described as: P (ν e →ν e ) ∼ = 1 − sin 2 2θ 13 sin 2 1.27∆m 2 31 (eV 2 )L(m) The value of θ 13 can be measured directly from the oscillation amplitude, inferred from an energy-dependent deficit in the number of observed neutrinos.
After a series of short and medium-baseline (∼100−1000 m) reactor neutrino experiments carried out in the 80s and 90s, a new generation of experiments has been operating for the past few years. Three different collaborations (Daya Bay [4], Double Chooz [5] and RENO [6]) have reported very precise measurements of the mixing angle θ 13 . As these experiments rely partially on the comparison of the observed antineutrino flux with respect to the expected one, a revision of the relatively old reactor flux models was performed in [7] and [8], becoming the new references and reducing the uncertainties at the 3% level. This re-evaluation of the flux led to the so-called reactor antineutrino anomaly [9], pointing at a possible short-baseline oscillation that would imply the existence of at least one sterile neutrino. While this suggests a possible underestimation of the reactor flux uncertainties, it does not impact the determination of θ 13 . Beyond this anomaly, Daya Bay, Double Chooz and RENO have reported very recently [10,5,11] an energy distortion around 5 MeV, which deviates from the expectation at about 3-4σ. Apart from reinforcing the idea of an underestimation of the flux errors, this experimental result has induced a world-wide effort in trying to understand the origin of this discrepancy.
This work reviews the observation of such a 5 MeV energy structure by the three current reactor experiments. The possible causes are described as well as the different experimental approaches carried out to identify its origin. The incompleteness of the reactor flux predictions is presented as the most likely explanation, so the different models developed so far are reviewed. Within those models, a number of possible sources of biases or error underestimations are listed, thus pointing at possible experimental ways to improve our current knowledge. This review is organized as follows: Sec. 2 describes the currentν e reactor experiments, relying on the flux predictions presented in Sec. 3; the observation of the energy distortion is reported in Sec. 4, followed by a critical analysis of its possible origin in Sec. 5; the reactor flux models are revisited in Sec. 6 as the most likely cause; Finally, Sec. 7 summarizes the state-of-the-art and discusses about the experimental data required to gain further knowledge on theν e reactor flux and the origin of the energy structure.

A new generation of reactor neutrino experiments
The most common way of detecting reactor neutrinos is via the inverse beta decay (IBD):ν e + p → n + e + . When this reaction takes place in liquid scintillator doped with ∼1% of Gadolinium, it produces two signals separated by about ∼ 30 µs: the first one due to the e + and its annihilation (prompt signal), and the second one due to the n capture in a Gd nucleus (delayed signal). This characteristic signature yields a very efficient background rejection. The prompt energy deposition (E e ) relates directly to the interacting antineutrino energy (Eν ): where M p , M n and M e are the proton, neutron and electron masses, respectively. The observed energy spectrum is the convolution of the reactor ν e flux and the IBD cross section. As shown in the left panel of Fig. 1, the mean energy of theν e spectrum is around 4 MeV, which corresponds to a prompt energy E e of ∼3 MeV. For this energy, the oscillation effect due to θ 13 starts arising at L ∼ 0.5 km and reaches the first maximum around 2 km, where the effect of θ 12 is still negligible as can be seen in the right panel of Fig. 1. Therefore, neutrino reactor experiments with short baselines offer a clean laboratory to search for θ 13 .  In spite of its characteristic signature, the IBD signal can be mimicked by the accidental and correlated backgrounds. The accidental background stands for the random coincidence of a positron-like signal coming from natural radioactivity, and the capture in the detector of a neutron created by cosmic muon spallation in the surrounding materials. The correlated background consists of events which may mimic both the prompt and the delayed signals of the IBD. Along with the stopping muons, the fast neutrons and cosmogenic isotopes, both generated in muon interactions, are the main sources of this background. Fast neutrons entering the detector lead to proton recoils, thus faking a prompt signal, before being captured. Muons crossing the detector can produce long-lived β-n decay isotopes, like 9 Li and 8 He. As the half-life of such cosmogenic isotopes is ∼100 ms, their decay cannot be associated to the muon interaction.
The sensitivity to the θ 13 -driven oscillation is optimized by detecting a deficit in the expected neutrino events around 2 km away from the nuclear power plant (far detector), as shown in right panel of Fig. 1. However, some of the largest systematic errors in reactor experiments arise from the uncertainties in the originalν e fluxes. In order to reduce them, a relative comparison between two or more identical detectors located at different distances from the reactors becomes critical. As originally proposed in [12], a near detector placed a few hundred meters away can measure the fluxes before any oscillation takes place. The comparison between the far and near detectors leads to a breakthrough in the sensitivity to θ 13 , as all the fully correlated systematic uncertainties cancel out. Further steps in the sensitivity optimization rely on reducing the relative detection efficiency and energy scale uncertainties of the detectors, as well as on minimizing the backgrounds.
Following the above ideas, a new generation of reactor experiments are running since 2010. In China, the Daya Bay experiment [4] has built a far site and two near sites meant to measure theν e fluxes from the 6 cores (17.4 GW th in total) of the three power plants existing in the area. The Double Chooz experiment [5] operates two identical detectors located 400 m (since 2014) and 1050 m away from the two 4.25 GW th reactor cores of the CHOOZ nuclear plant in France. RENO [6] also consists of two identical detectors measuring the antineutrino fluxes generated at the 6 cores (17.3 GW th in total) of the Youngwang nuclear plant in South Korea. Although there are some differences in the detector designs of the tree experiments, they all rely on the same principles and technology. The detectors are divided into three concentric volumes: the target (the inner-most volume filled with Gd-doped liquid scintillator), the γ-catcher (filled with undoped liquid scintillator) and the buffer (filled with mineral oil). The light produced by interactions in the liquid scintillator is read out by a number of photomultiplier tubes (PMTs) located in the buffer walls.

Antineutrino flux prediction
In a nuclear reactor, about 6 antineutrinos from β-decays are generated per fission, releasing an average energy of about 200 MeV. As the unstable fission products are rich in neutrons, they undergo β − decays generating a nearly pure electron antineutrino flux. Only four isotopes, whose fission products can produceν e with energies above the IBD threshold (1.8 MeV), contribute to more than 99% of the flux: 235 U, 239 Pu, 238 U and 241 Pu. However, such a flux consists of a superposition of thousands of β-decay branches. A fraction of the neutrons produced in the 235 U fissions is captured by 238 U, giving place to mostly 239 Pu. Thus, the core burns 235 U while accumulating 239 Pu as it is operated, in the so-called burnup process. Apart from these two main isotopes which make up about 90% of the flux, the 241 Pu and 238 U fissions contribute to the remaining 10%. From a practical point of view, this implies that an accurate reactor flux prediction relies on two main aspects: 1) the simulation of the time evolution of the core fuel composition (i.e, the contributions of each one of the four main isotopes), and 2) the knowledge of the β spectra associated to the decay chains of the fission products.
In order to compute the expected neutrino flux in a reactor experiment like Daya Bay, Double Chooz or RENO, three main ingredients need to be taken into account: 1) the detector-related normalization terms, 2) the reactor flux as a function of time, and 3) the IBD cross section. In absence of oscillations, the number of expected antineutrinos from a nuclear core can be described as: where ǫ is the detection efficiency, N p is the number of protons in the target, L is the distance to the center of the reactor, and P th is the thermal power. E f is the mean energy released per fission: where α k is the fractional fission rate of the k th isotope (k = 235 U, 239 Pu, 238 U, 241 Pu). The mean cross section per fission σ f k is defined as: where S k (E) is the reference spectrum of the k th isotope and σ IBD is the inverse beta decay cross section. The three variables P th , E f and σ f are time dependent, with E f and σ f depending on the evolution of the fuel composition in the reactor and P th depending on the operation of the reactor.
The current reactor experiments have used the reference spectra S k (E) from [7,8] as an input to their oscillation analyses. However, as far as the determination of these spectra is concerned, S k (E) can be expressed as the sum of the contributions from all the fission products (N f ): where A f is the activity of the fission product and the spectrum S f (E) of each fission product is in turn a sum of N b β-branches connecting the ground state (or an isomeric state) of the parent nucleus to different excited levels of the daughter nucleus: being BR b f and E b 0f the branching ratio and the endpoint energy of the b branch of the f fission product, respectively, and Z f and A f the charge and atomic number of the parent nucleus. It is worth noticing that Equations (5) and (6) are valid for both electron and antineutrino spectra. The beta decay spectrum S b f for a single transition in a nucleus with endpoint energy E b 0f = E e − E ν is then: where S 0 is a normalization constant taking into account the phase space [7,13], F (E e , Z f , A f ) is the Fermi function accounting for the Coulomb interaction of the outgoing electron with the charge of the daughter nucleus, and C(E e ) is a shape factor [14] for forbidden transitions due to additional lepton momentum terms (C(E) = 1 for allowed transitions). Beyond these terms, some additional effects need to be considered for precision studies: this is the role of the δ(E e , Z, A) factor. It accounts for the radiative (R), finite size (FS) and weak magnetism (WM) corrections: The R corrections are due to the emission of virtual and real photons by the charged particles present in the β-decay, and it is computed in [15]. The FS correction accounts for the finite size of the nucleus, as the electric charge and hypercharge are not point-like [16,17]. The WM term refers to the induced current yielding the largest contribution to the shape of the β spectrum [16,18]. Finally, the simplified form from Vogel and Beacom [19] can be used to describe the IBD cross section: where and m e and E e + are the positron mass and energy. The variables m n and m p are the masses of the neutron and proton with ∆ = m n − m p . The constant K is inversely proportional to the neutron lifetime. Using the MAMBO-II measurement of the neutron lifetime [20] leads to K = 0.961 × 10 −43 cm 2 MeV −2 .

Reactor flux models
In order to predict the reference spectra S k (E), two different approaches are developed. The ab initio or summation method takes advantage of the available information on the β decays of each fission fragment (nuclear databases), summing over each nuclide's individual spectrum to obtain an aggregate spectra. On the other hand, the so-called conversion method exploits the aggregate β spectra measured in the Institut Laue-Langevin (ILL) [21,22,23,24], fitting the data to a set of virtual β branches and converting the result into the corresponding antineutrino spectra (e.g. [25]). While it is worth noticing that both methods rely on measured β spectra, the conversion approach yields the most precise results since the uncertainties are constrained by the ILL measurements. Although this is reviewed in this work, the errors associated to the conversion method have been claimed to be at the level of 2-3%. As the nuclear databases are known to suffer from a lack of relevant data (concerning both β decays and fission yields) and from the need of more precise measurements, the summation method provides typically an envelope error of about 10-20%. Recently, there have been improvements in both the conversion and summation techniques [7,26], being one of the main goals to optimize the results from the current reactor experiments. Given the limitations of the ab initio approach, the reactor antineutrino spectra have been estimated historically relying on the total electron spectra associated with the beta decays of all fission products of 235 U, 239 Pu, and 241 Pu. Such β-spectra were obtained at ILL by irradiating thin target foils of these isotopes with thermal neutrons. As 238 U nuclei undergo fission with fast neutrons, the associated spectrum could not be measured at that time and therefore its prediction has been typically based on the summation method. Same applies to all spectra above 8 MeV, as the ILL measurements were performed only up to this energy. In [7], a mixed approach has been developed combining the precise reference of the electron spectra from ILL with the physical distribution of beta branches of all fission products provided by the nuclear databases. This new analysis has provided a better handle on the systematic errors of the conversion, and a new set of reference spectra for 235 U, 239 Pu, 241 Pu and 238 U (although the latter is still based on a purely summation technique). While the shapes of the spectra and their uncertainties are found to be comparable to that of the previous analysis of the ILL data, the normalization is shifted by about +3% on average, thus leading to the reactor neutrino anomaly. The re-evaluation of short-baseline reactor data in the light of these new reference spectra reveals a deficit in the number of observed antineutrinos, which might be explained in terms of sterile neutrino oscillations. One of the main reasons for this normalization shift is the treatment of the corrections δ(E e , Z p , A p ) in Eq. 7, and in particular the WM term. These corrections have been further investigated in [8], deriving a consistent set of reference spectra in both shape and normalization. To complete the picture of the state-of-the-art conversion method, it must be noticed that the cumulative β spectrum of the fission products of 238 U has been finally measured in [27], in the range from 2.875 MeV to 7.625 MeV.
Despite the larger uncertainties, the summation method is still a powerful tool to predict reactor antineutrino fluxes. To start with, this approach provides estimations of the reference spectra which are independent from the measurements at ILL. As these measurements are unique, a cross check based on nuclear databases is specially valuable. The summation method is also the only way to predict the antineutrino energy spectra beyond 8.0 MeV for 235 U, 239 Pu and 241 Pu, and beyond 7.6 MeV for 238 U. While it is true that some authors (e.g. [7,8]) have provided polynomial parameterizations of the spectra that can be used to extrapolate the predictions above 8 MeV (as typically done by reactor experiments), such an extrapolation is not physically motivated and the associated error cannot be estimated in a robust way. Furthermore, the binning of 250 keV in which the ILL spectra were originally published is large enough to hide possible structures coming from some contributions to the reactor fluxes. Taking into account all these considerations, the summation technique has become a main tool to shed light on the reactor antineutrino spectra and the nuclear databases have been recently improved by new β-decay measurements. In particular, a new set of reference energy spectra has been obtained in [26] taking into account the new measurements of the 102;104;105;106;107 Tc, 105 Mo and 101 Nb nuclei. These measurements are taken with the Total Absorption Technique (TAS), insensitive to the Pandemonium effect [28] which typically affects the γ spectrometry with Ge detectors. Beyond the relevant improvement in the summation-based reference spectra, the work in [26] highlights the need of new TAS measurements.
In order to compare the state-of-art reference spectra, the ratio of the summation spectra derived in [26] to the conversion-based predictions ( 235 U, 239 Pu and 241 Pu from [8] and 238 U from [27]) is shown in Fig. 2. Both set of predictions agree within 1σ, although no strong conclusions can be settle Pu 241 Figure 2: Ratio of summation-derived spectra (from [26]) to the start-of-art conversion-derived spectra (from [8] and [27]), as a function of theν e energy. The shadowed band shows the 1σ error.
given the large error bands.

Observation of an energy structure around 5 MeV
In the Neutrino 2012 conference, the RENO collaboration mentioned the observation of an excess in the number of antineutrinos between 4 and 6 MeV of the positron energy spectrum [29], with respect the to flux prediction. In Neutrino 2014, both RENO and Double Chooz collaborations reported and quantified such an excess [11,30]. Daya Bay also presented a similar energy structure at the ICHEP 2014 and at NuTel 2015 [31] conferences. The first journal publication of this effect has been provided by Double Chooz in [5]. Using the data from the far detector, a best fit value of the flux normalization of 9±2% (with respect to the central value prediction) is quoted between 4.25 and 6 MeV. This translates into a significance of 3σ. The energy distortion is consistent with the ones presented in previous Double Chooz publications [32,33,34], where a significant measurement of the data-prediction discrepancy was not possible due to the limited statistics and the non-optimized detector energy response. Daya Bay has also released a paper [10] on the reactor neutrino flux measurement, concerning both the normalization and the energy spectrum. The measured prompt energy spectrum in the near detectors shows a deviation from the reactor models with a significance beyond 2σ over the full energy range, and around 4σ between 4 and 6 MeV. The excess in this energy range has been estimated to be about 1% of all events in both the near and far detectors. Once corrected by the θ 13 oscillation effect, the energy spectra measured in the near and far detectors are consistent. Two reactor predictions have been considered in Daya Bay, one based on the conventional ILL models and another one based on the re-evaluations provided in [7,8]. The disagreement between the data and the prediction arises in both cases, being the significance of the deviations very similar. The significance has also been computed adopting two different approaches, one relying on the contribution of the χ 2 of each energy bin and another one on the p-values within local energy windows, yielding consistent results. In addition, the RENO collaboration has shown in the the proceedings of Neutrino 2014 [11] the energy structure in the 4-6 MeV energy range for both the near and far detectors. The observed excess ofν e is consistent among the two, and the significance of the excess at the near detector is estimated to be 3.5σ. The excess ofν e with respect to the total expected flux is quoted as 2.3±0.4 (data) ± 0.5 (prediction)% for the near detector, and 1.8±0.7 (data) ± 0.5 (prediction)% for the far detector.
The ratio of the background-subtractedν e candidates spectrum to the non-oscillation prediction is shown in Fig. 3 for the Daya Bay and RENO near detectors and for the Double Chooz far detector. The excess reported by the three collaborations amounts to about 10% over the expected number ofν e in this energy range, and both RENO and Daya Bay have observed consistent structures in their near and far detectors. It might be argued that the excess observed in RENO is larger than in Daya Bay and Double Chooz. However, given the discrepancy between RENO and Daya Bay in the 2-4 MeV range (see Fig. 3), such a difference in the amplitude of the distortion might be due to differences in the flux predictions beyond the 4-6 MeV window. In order to compare the observed distortions in a robust and quantitative way, theν e prediction of the three reactor experiments should be based on the same reactor model (which comprises not only the reference spectra S k (E), but also the simulation of the reactor core evolution, the treatment of the spent fuel, etc). It is also worth noticing that the significances of the excess quoted by the three experiments cannot be directly compared, as they are computed in different ways. Daya Bay normalizes the predicted spectrum to the observed number of events, thus evaluating the discrepancy in terms of the energy spectrum between 4 and 6 MeV, and not the total rate. On the other hand, Double Chooz performs an evaluation based on the total predicted and expected rates in the 4.25-6.00 MeV energy window. Independently of how the distortion significance is estimated, it is limited by the uncertainties in the flux prediction, which are at the level of 2-3% for both the rate and the spectral shape.
As the three experiments detect the reactorν e in the same way and with Best-fit sin Figure 3: Ratio of background-subtractedν e candidates to non-oscillation prediction, as a function of the IBD prompt energy. Left: Daya Bay and RENO ratios forν e candidates observed at the near sites. Right: Double Chooz ratio forν e candidates measured at the far detector. The shadowed region represents the typical reactor error derived from the [8] and [27] referenceν e spectra, which is dominant for the considered energy region. very similar detectors, the possible causes of this energy structure are common and might include, in principle, detector and/or background issues. An explanation in terms of the reactor flux prediction (incompleteness of the model or underestimation of the uncertainties) would be also correlated among the three experiments, as they all rely on the conversion method to obtain the 235 U, 239 Pu and 241 Pu antineutrino spectra. There are however two differences. Firstly, Double Chooz uses in [5] the measurement in [27] to derive the 238 Uν e spectrum, while Daya Bay and RENO take the summation-based spectrum from [7]. Secondly, Double Chooz constrains the flux normalization to the measurement in Bugey4 [35], taken 15 m away from the core. This is why Double Chooz quotes a flux normalization error of 1.7%, thus reducing the errors in Daya Bay and RENO (2.7% and 2.0%, respectively). Beyond this, Daya Bay has also explored some variations to the reference model but concluded that the structure still remains. In particular, the local distortion around 5 MeV cannot be described extending the reactor model with a single β-branch or a mono-energetic line. While this disagreement between data and flux models needs to be investigated, it must be noticed that the impact on the θ 13 mixing angle is negligible. As demonstrated in [5] for Double Chooz, even with only the far detector being used for the oscillation analysis, the θ 13 measurement is not affected by the energy distortion. This can be easily understood since the amplitude of the θ 13 -driven oscillation is vanishing around 5 MeV. In the case of Daya Bay and RENO, whose analyses involve both near and far detectors, the impact of the energy structure is even smaller as the role of the flux prediction is not as relevant due to the inter-detector comparison.

Possible sources of the energy structure
The spectral shape of the energy structure at 5 MeV cannot be produced by any standard neutrino oscillations scenario, even considering sterile neutrinos. In particular, it has been observed by Daya Bay and RENO at two different baselines (the near and far detectors). Therefore, it can be assumed that the discrepancy between the data and the flux models might be due to one of these reasons: 1) the existence of non-standard IBD interactions, 2) a detection issue distorting the energy scale, 3) an unaccounted background source, and 4) missing contributions to the reactor models. Being the reactorν e spectrum uncertainty of the order of 2-3% and the maximum deviation between data and prediction around 10%, the current reactor experiments cannot establish this discrepancy beyond a significance of ∼4σ. However, the available data allow for a dedicated analysis on the possibles causes of the discrepancy. The current reactor experiments have been capable of reinforcing the case for a reactor-model explanation, while disfavoring other possible causes, namely the misinterpretation of the detector response and the incompleteness of the background model. The dedicated studies addressing the possible sources of the prompt energy spectrum distortion are described below.

Antineutrino interactions
An unaccounted or non-standard neutrino interaction in the detectors of the reactor experiments might lead to an excess of observed neutrinos. In the energy range of reactorν e (below 10 MeV), the typical cross section of charged and neutral current neutrino interactions follow an increasing pattern with energy after a given threshold. This kind of trend can hardly explain a bump-like excess around 5 MeV in the positron energy spectrum. Within the target volume of the detectors, the antineutrinos can interact basically with H, C and Gd. In the γ-catcher, only with H and C. As Double Chooz has observed the energy structure using neutrons captures in Gd and in H [34], an unaccounted interaction with Gd can be excluded. The antineutrinos might interact with some C isotope with enough energy to separate one neutron, remaining the final nuclei in excited state. The combination of the de-excitation γ and the neutron might mimic the IBD signal. However, the rate of such process should be rather small (as there are not empirical evidences in the current experiments) and could not explain the ∼10% excess.

Energy scale
A detector-related issue affecting the energy scale might also explain, in principle, the energy structure. However, some studies performed by Daya Bay and Double Chooz rule out this possibility. In [5], the accuracy of the energy scale around 5 MeV has been confirmed by spallation neutrons captured on carbon, which occur predominantly in the γ-catcher volume as the capture cross section is smaller than on Gd. The C captures result in an energy peak at 5 MeV, whose agreement between data and MC simulation has been found to be within 0.5 %. Along the same lines, Daya Bay, Double Chooz and RENO have also checked the energy reconstruction by means of the β decays of 12 B collected in data, showing no energy distortion when compared to the corresponding simulation. This is consistent with the fact that any nonlinear effect, due to the scintillator properties or the electronics response, is observed for energies above 4 MeV. Furthermore, the energy resolution estimated with the collected data is also in good agreement with that of the Monte Carlo.

Background model
The events found in the 4-6 MeV range fulfill all the IBD characteristics, in particular concerning the neutron capture time and distance distribution, and the spatial distribution of the prompt signals. Thus, the three collaborations disfavor the hypothesis of an unaccounted background contribution. In addition, the reactor off data taken in Double Chooz allows for an independent and inclusive background measurement, thus accounting even for possible unknown sources [36]. The measured total rate in [5] according to the candidates selection cuts is computed to be 0.75±0.37 events/day. While keeping the independence with the background model, this background measurement is slightly modified by means of a Reactor Rate Modulation (RRM) [37] fit: 0.90 +0. 43 −0.36 events/day. These total background measurements are lower than the sum of the individual background sources accounted for in the background model (accidental coincidences, fast-neutrons/stopping-muons, and cosmogenic isotopes): 1.6 +0. 41 −0.17 events/day (1.7σ discrepancy with respect to the reactor off measurement). Therefore, the existence of an unaccounted background source, leading to the excess around 5 MeV, is strongly disfavored. Beyond this comparison with the inclusive background measurement, dedicated studies with reactor on and reactor off data have been developed to look for direct indications of unknown background sources. No significant evidences have been found.

Reactor flux model
The remaining possible cause of the energy distortion is that one of an additional reactorν e component beyond the current model. In particular, if the excess around 5 MeV is due to an unaccounted reactor contribution, it must be correlated to the reactor power. On the other hand, if it is due to an unknown background, the rate of the excess should be independent of the power. Such a correlation has been demonstrated by Daya Baya, Double Chooz and RENO, by estimating the excess for different reactor powers. Daya Bay has shown in ICHEP 2015 the time stability of the prompt energy spectrum and the time distribution of events for two different energy windows (4.5-5.5 MeV and 3.0-4.0 MeV), proving that the structure remains the same over time, and thus for different conditions of the reactors operation. In [5], the Double Chooz collaboration shows the correlation of the excess with the reactor power in a flux-model independent way, by parameterizing the spectrum and measuring an effective excess for different reactor conditions. Consistent results are found forν e candidates obtained with neutron captures in Gd and H. RENO has also reported in Neutrino 2015 and NuTel 2015 such a correlation by means of the measurement of the excess for different reactor powers, as shown in left panel of Fig. 4.
As the Double Chooz RRM analysis utilizes the correlation between the observed rate and the thermal power to derive both the mixing angle θ 13 and the total background rate, it can be used to test the hypothesis of a bias in the flux prediction. In particular, it can confront the data to the background model and the flux model at the same time, thus providing indications about the most likely cause of the energy structure. In [5], five independent RRM fits have been carried out in different energy regions, constraining sin 2 2θ 13 to the best fit value in [4] while leaving as free parameters both the total background rate and a flux normalization term (with respect to the central value of the flux model). The best fit values of the background rate are fully consistent with both the background model and the reactor off measurement, while the best fit values for the flux normalization deviates (2σ) from the prediction in the 4.25-6.00 MeV window, as shown in the right panel of Fig. 4. This result is consistent with the reported correlation between the excess and the thermal power, thus reinforcing the case for a flux model bias and disfavoring again the background model as the source of the energy distortion. If one constrains the total background rate to the background model, the discrepancy between the flux model and the RRM best fit value is increased to 3.0σ.

Reviewing the reactor flux predictions
As discussed above, the most likely explanation for the energy structure is that one of unaccounted contributions in the reactor flux models. Hereafter, this work assumes this is the unique source of the data-prediction disagreement around 5 MeV. Under this well motivated assumption, the different approaches to the flux estimations need to be reviewed from a critical point  [11]). Right: RRM best fit values of the reactor flux normalization (with respect to the central value prediction) in the far detector of Double Chooz (data from [5]). Results with and without the background (BG) model constraint are shown with empty squares and solid dots, respectively. of view. To start with, the errors of the conversion-based reference spectra must be re-evaluated somehow, as the current quoted uncertainties do not cover the observed energy distortion. This is obviously related to the identification of the possible missing pieces in the reactor models. Such contributions in the conversion-based flux might be related to: 1) the aggregate β measurements, 2) the conversion procedure itself, and 3) the nuclear corrections, mostly related to the forbidden decays. To shed light on all these possibilities, the summation-based predictions might play a major role, but the associated limitations need to be taken into account. All these aspects are discussed below.

The aggregate spectra
As the 235 U, 239 Pu and 241 Pu reference spectra obtained by means of the conversion method rely on the ILL data, any issue affecting the ILL spectrometer would propagate to the reactor flux predictions. In principle, biases in both the overall normalization and the energy reconstruction (or the associated errors) might be possible, thus giving rise to the reactorν e anomaly and the energy structure around 5 MeV, respectively. As pointed out in [38] following a summation approach, the presence of a bump between 5-7 MeV in both the calculated electron and antineutrino spectra might be an indication of an artifact in the original ILL measurements rather than an effect of the conversion method.
Beyond the ILL data, the summation-derivedν e spectra from 238 U (used as reference in Daya Bay, RENO and first publications by Double Chooz) might be considered a candidate to explain the energy distortion, given the large associated errors. Because of the different experimental setups, RENO reports that about 12% of the fissions are due to 238 U, while Daya Baya quotes only 7.6%. As the energy structure in RENO is about 50% larger than in Daya Bay, this might indicate that the 238 U fissions are contributing to it. This isotope is indeed responsible for about 20% of theν e flux between 4 and 6 MeV. In [39], it has been reported that two different databases predict, within the summation scheme, bumps in the region of interest. However, the amplitude is not large enough to cover the structure observed in reactor experiments. Furthermore, the analysis in [39] has not considered the work in [27], where the 238 U aggregate spectrum has been measured. On the contrary, Double Chooz has used [27] to predict theν e flux and has found that the energy bump remains, with roughly the same amplitude.

The conversion procedure
The conversion technique has been reviewed extensively in the literature since the first antineutrino predictions based on the ILL data [25]. As discussed previously, the method was recently improved in [7,8]. While the different approximations have been performed to fit the data to a number of virtual β branches, the results have been consistent as far as energy shape and error budget are concerned (this is not the case of the overall normalization). This leaves small room for a possible issue in the technique itself. However, the ILL reactor is different to that ones currently used in the reactor experiments. The neutron flux spectra at typical pressurized water reactors (like the ones in Daya bay, Double Chooz and RENO) are harder in energy than the thermal spectrum of the ILL reactor. As highlighted in [39], this opens the possibility of epithermal neutron contributions to the 235 U, 239 Pu, 241 Pu and 238 U fissions, resulting in a shoulder at 5 MeV in theν e spectrum. However, since there are not fission yield measurements for the nuclei that dominate that energy region, this hypothesis is hard to demonstrate or refute.

Nuclear corrections
The uncertainties quoted in [7] have been revisited in [13], since they lead to a significance of the reactor neutrino anomaly of about 3σ. As described in Sec. 3, an antineutrino spectrum can be estimated from a beta spectrum if the linear combination of operators involved in the decay, the endpoint energy and the nuclear charge are known. However, the fission β spectra involve about 6000 decays, being forbidden about 30% of them. This im-plies that some assumptions are needed when deriving the reactorν e flux, given the limited knowledge on the structure of the forbidden transitions. Such assumptions affect eventually the error budget of the predictions in both the conversion and summation methods. In [13], it has been noticed that different treatments of the forbidden transitions (and the associated δ corrections) can lead to antineutrino spectra that differ both in shape and magnitude at about 4%. In particular, if all forbidden decays are treated as allowed transitions, the antineutrino spectra are increased, thus leading to the reactor neutrino anomaly reported in [9]. However, this is not always the case if different approaches for the forbidden transitions are adopted. It is concluded that uncertainties in theν e predictions are about 4%, implying an increase of roughly a factor 2 with respect to the estimations in [7,8]. Along the same lines, the effect of first forbidden transitions on the β-decay neutrino spectra is analyzed in [40] by performing microscopic nuclear structure calculations. The authors conclude that these decays may be responsible for a fraction of the deficit of neutrinos observed in the reactor experiments. Although the works in [13,40] are addressing the issue of the reactor neutrino anomaly, the conclusions apply also to the energy structure observed around 5 MeV: the disagreement between the data and the flux prediction might be due to the forbidden transitions contributing to that energy region.
Concerning the nuclear corrections to be applied to the forbidden transitions, three relevant points have been highlighted in [39]. First, it has been noticed that several of the β decays contributing to the bump region have a total angular momentum and parity which involve no WM correction. This increases the flux predictions with respect to [7,8,26], where this fact has not been taken into account. Second, the shape factor C(E e ) is not the same for all forbidden transitions: in [7], the C(E e ) corresponding to a unique forbidden transition has been assumed for all the cases. Finally, the FS corrections for these transitions applied in the literature are always approximated. Despite these considerations, a more accurate treatment of the forbidden decays performed in [39] cannot account for a significant fraction of the energy structure.

The limitations of the summation method
In order to get insights on the above topics, the summation method is a powerful handle as it provides an ILL-independent set of reference spectra and a tool to estimate the effect of the different nuclear corrections and assumptions. However, it is worth remarking that the errors associated to this technique are too large to establish any conclusion. Furthermore, the use of different databases can lead to somehow different conclusions. As an example, the ENDF/B.VII.1 compiled nuclear data [41] has been used in [38] to derive theν e spectra, yielding an energy bump in the antineutrino energy in the 5-7 MeV region (E e =4-6 MeV). However, the ENDF/B.VII.1 data used in the analysis is not taking into account the new TAS measurements described in [26]. By using the reference spectra from [26] or the ENDF/B.VII.1 database combined with the new TAS measurements (hereafter updated ENDF/B.VII.1) [42], the bump is significantly reduced. A detailed comparison of the updated ENDF/B.VII.1 and JEFF-3.11 [43] decay libraries has been presented in [39]: in the case of the JEFF-3.11-based results, the bump is totally removed. Finally, it must be noticed that the same kind of limitations of the available nuclear data arise when considering the reactor neutrino anomaly (i.e, the flux normalization). As an example, Daya Bay has measured an absoluteν e rate in good agreement with the current world average and with the ENDF/B-VII.1 prediction (i.e., no indication of anomaly). However, comparing the Daya Bay data with the JEFF-3.1.1 estimations yields a deficit in theν e rate, thus suggesting the anomaly.

Summary and discussion
The three current antineutrino disappearance reactor experiments (Daya Bay, Double Chooz and RENO) have observed an energy distortion around 5 MeV in the prompt energy spectrum (∼6 MeV inν e spectrum), deviating from the predictions at a ∼4σ level as an excess in the number ofν e . The structure is observed in the near detectors of Daya Bay and RENO (with baselines around 300 m) and in the far detectors of the three experiments (with baselines around 1-2 km), so it cannot be explained in terms of any standard neutrino oscillations. Given the correlation of the excess with the reactor power, the three collaborations conclude that the most likely explanation is an incompleteness or bias in the reactor flux models. Although the origin of the energy structure might not be related to the reactor antineutrino anomaly (that might be described in terms of sterile neutrino oscillations), both features reinforce the case for a revision of the current reactor flux predictions. A possible underestimation of the error budget in the reactor models, due to unknown or not well described contributions, has to be considered.
There are two general methods to estimate the reactor fluxes as a composition of theν e spectra from the main fissile isotopes. The most precise one, and the state-of-the-art reference for reactor experiments, is the conversion method: it relies on the β aggregate spectra measured in the ILL, fitting the data to a set of virtual branches and converting the results into the correspondingν e spectra. Currently, this method quotes an uncertainty of 2-3%. The second approach is the so-called summation method: it builds theν e spectra as the sum of each nuclide's individual β spectrum, according to the available information in the nuclear databases. This technique typi-cally yields an error envelope of 10-20%. While uncertainty associated to the summation-based spectra covers the ∼10% deviation observed in the experimental data around 5 MeV of the prompt energy spectrum, the error of the conversion-based spectra does not. The later method might be affected by: 1) an issue in the ILL data, 2) intrinsic limitations of the technique, and 3) nuclear effects or uncertainties not accounted for. The forbidden transitions might play a major role in the later aspect.
The available reactor data can be used to shed some light on the puzzle of the flux predictions. In particular, one can perform fits to the observed ν e spectrum for different sets of parameters and assumptions used in the predictions (like the number of β branches or the treatment of forbidden transitions). The fit results can help to identify or rule out possible contributions to the prompt energy shoulder at ∼5 MeV. It is also possible to take theν e spectra from Daya Bay, Double Chooz and RENO, and deconvolute them back to the corresponding β spectra, which in turn can be confronted to the ILL data. The observation of the same structure in the β spectra would be an indication of a bias in the ILL data. Finally, the summation method can be used to review the error on the predicted spectra by analyzing the effect of the different approximations concerning the nuclear corrections and the forbidden transitions. If it is concluded from this kind of studies than the error in theν e spectra is larger than currently assumed (as suggested by some authors), the discrepancy between the observed data and the models would need to be reevaluated.
Despite the above considerations, the available reactor data is not enough to find out the actual origin of the energy distortion. The same applies to the current nuclear data concerning β decays. However, the situation might improve once the near-future campaign of very short-baseline reactor experiments (meant to explore the possibility of sterile neutrino oscillations) starts delivering data. In order to rule out the possibility of an issue in the original ILL measurements, a new aggregate β spectrum is the most direct approach. While this new measurement would be valuable, it might not be the ultimate solution to the origin of the discrepancy between the data and the models. Beyond cross-checking the ILL measurement, it is also worth exploring the limitations of the related conversion procedure, specially regarding the neutron spectra in different types of reactors. The comparison of aggregate β spectra measured in a very thermal reactor and in a reactor with a harder neutron spectrum, as proposed in [39], would suffice to quantify the impact on the predictedν e flux.
The available nuclear data is neither capable of identifying the origin of the energy structure by means of the summation method. In order to improve the precision and accuracy of the summation-derived predictions, new β decay measurements are needed. In particular, improving the knowledge on the forbidden β transitions is crucial. The summation-method itself can be used to define the list of most relevant transitions to be measured, by tagging the nuclei that contribute the most to the energy region of interest. As a matter of fact, the main contributors have been identified in works like [38,42]. As most of the relevant transitions are first forbidden, the β spectrum needs to be measured with high precision so the shape correction factor can be explored. As already demonstrated in the literature, the TAS technique has arose as the best option for several transitions. A new campaign of measurements will boost the capabilities of the summation method, thus becoming a major tool to resolve the nature of the reactorν e energy structure.

Note added in proof
After the submission of this manuscript, the RENO collaboration has released a paper [44] where the observation of the energy structure is described. Although with 500 live days of data instead of 800, the paper accounts for the results presented in the Neutrino 2014 and NuTel 2015 conferences (cited in this review), including the figures of the prompt energy spectrum at both far and near detectors.