Improved STEREO simulation with a new gamma ray spectrum of excited gadolinium isotopes using FIFRELIN

The STEREO experiment measures the electron antineutrino spectrum emitted in a research reactor using the inverse beta decay reaction on H nuclei in a gadolinium loaded liquid scintillator. The detection is based on a signal coincidence of a prompt positron and a delayed neutron capture event. The simulated response of the neutron capture on gadolinium is crucial for the comparison with data, in particular in the case of the detection efficiency. Among all stable isotopes, 155Gd and 157Gd have the highest cross sections for thermal neutron capture. The excited nuclei after the neutron capture emit gamma rays with a total energy of about 8MeV. The complex level schemes of 156Gd and 158Gd are a challenge for the modeling and prediction of the deexcitation spectrum, especially for compact detectors where gamma rays can escape the active volume. With a new description of the Gd (n,γ\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$ \gamma$\end{document}) cascades obtained using the FIFRELIN code, the agreement between simulation and measurements with a neutron calibration source was significantly improved in the STEREO experiment. A database of ten millions of deexcitation cascades for each isotope has been generated and is now available for the user.

The Stereo experiment [1] detects electron antineutrinos produced in a compact research reactor at the Institut Laue-Langevin (ILL) in Grenoble, France. With the Stereo data, it is possible to test the hypothetical existence of a sterile neutrino in the eV mass range [2]. Antineutrinos are detected via the inverse beta decay reaction (IBD)ν e + p → e + + n in a gadolinium (Gd) loaded organic liquid scintillator (LS) [3]. The positron ionization gives rise to a prompt signal -related to the antineutrino energy-, while the neutron thermalizes and diffuses in the liquid. The gamma emission following the radiative capture of the neutron creates a delayed signal after typical coincidence time of a few tens of μs.
The capture time can be strongly reduced by the addition of Gd with its very high cross section for thermal neutron capture (∼ 10 5 barn for some of the isotopes) to the LS. At a Gd-concentration of 0.2 wt.% in Stereo, the average capture time is 18 μs, one order of magnitude lower than for an unloaded scintillator where capture is mainly by hydrogen (H). More than 80% of the neutrons are captured on Gd. The resulting excited Gdnuclei decay to the ground state by emitting on average four gammas with a total energy of about 8 MeV. This is well above the typical energies of the natural radioactivity backgrounds and the 2.2 MeV line of H(n, γ), providing a clean detection channel. Therefore, several past, present, and upcoming neutrino detectors take advantage of the Gd-loaded LS technology. All of these experiments rely on a precise knowledge of the Gd-deexcitation process, where the gamma multiplicity and single energies are especially relevant for segmented and/or compact detectors such as Stereo. These quantities play a crucial role in the reconstructed spectrum of the experiments, since high energy gammas are more likely to escape the detector or leave only a part of their energy via Compton scattering, and populate thereby the tail of the reconstructed peak. The DANSS experiment, for example, reported some tension between calibration and simulation data in the low energy part of the Gd spectrum [4]. To determine the neutron detection efficiency in the Daya Bay experiment, four different models are used to estimate the gamma energy and multiplicity distributions [5]. An accurate modeling of these distributions is also of major relevance in water Cherenkov detectors with Gd-loading [6], since the light production is very sensitive to the energies of the single gammas due to the Cherenkov threshold. Therefore, a good understanding of the cascade for the most relevant isotopes 156 Gd and 158 Gd is of primary importance to reduce the systematic uncertainties. However, none of these isotopes have complete experimentally known nuclear level schemes and branching ratios up to 8 MeV. The cause of small discrepancies between the Stereo data and simulation was identified to be largely due to an inaccurate description of the gamma deexcitation of the Gd nuclei. By using 156 Gd * and 158 Gd * gamma cascades from the nucleus deexcitation code Fifrelin developed at CEA-Cadarache (France) these discrepancies are reduced.
The Fifrelin code is developed for the evaluation of fission data and has already proven its ability to make accurate predictions regarding neutron and gamma properties [7][8][9]. It provides important information on crucial parameters like the gamma multiplicity and the particle energies. The code is made of two parts: one is the assignment of the initial state of the fission fragment and the other is the deexcitation process. For the Stereo experiment, only the deexcitation part is used. Once the initial states are given, a set of nucleus level schemes is sam- pled allowing to take into account nuclear structure uncertainties. The deexcitation processes are then performed within a Monte-Carlo Hauser-Feschbach framework, based on Bečvář's algorithm [10], and extended to the n/γ emission by Régnier [11]. After a thermal neutron capture (E n = 25meV), the 156 Gd * and 158 Gd * nuclei have an excitation energy approximately equal to the compound nucleus neutron separation energy S n (E n S n ), equal to 8.536 MeV and 7.937 MeV, respectively. Knowing the 155 Gd, 157 Gd ground state spin-parity J π = 3/2 − , the selection rules give to the excited nuclei a parity −1 and two allowed spins of 1h and 2h. Following the latest nuclear data evaluations [12][13][14] a 2h spin is assigned to both nuclei, corresponding to their first resonance spin. In Fifrelin, a realization of the nuclear level scheme ( fig. 1) uses all the experimental knowledge and the missing information comes from nuclear models: -for E ≤ E RIPL , all the energy levels are known and are retrieved from the Reference Input Parameter Library (RIPL-3) database [15]. If a spin and/or a parity are missing, the code samples them from the momentum and parity theoretical distributions detailed below. -for E RIPL < E ≤ E limit , only a few levels are experimentally determined. Additional discrete levels are then sampled until the level number matches the theoretical level density. This is done until E limit , defined by a nuclear level density set to 5 · 10 4 MeV −1 (default value). -for E > E limit , the number of levels is innumerable and corresponds to the continuum. Therefore, levels are gathered in energy bins (dE = 10 keV by default) having a specific J π given by a model.
The values of E RIPL , E limit and E M for the relevant Gd isotopes are given in table 1. The theoretical nuclear level density ρ(E, J, π) used to complete the level scheme is with ρ tot (E) the total nuclear level density, P (J|E) the energy dependent angular momentum distribution and P (π) the parity distribution. The positive and negative parities are assumed to be equiprobable, P (π = ±1) = 0.5. The angular momentum distribution is with σ 2 (E) the spin cut-off parameter defining the dispersion of the nucleus angular momentum. More information on its origin and its parametrization can be found in [15,16]. The total nuclear level density follows the Composite Gilbert Cameron Model (CGCM) [17], using the Constant Temperature Model (CTM) at low energy and the Fermi Gas Model (FGM) at high energy with E M the energy where the CTM and FGM nuclear level densities match along with their derivatives. The CGCM parametrizations used here can be found in [15]. During a deexcitation step, all the transition probabilities Γ p (i → f, α) to go from a given initial state (i) to a final state (f ) in emitting a particle p with given properties (α) are computed. Then, one transition is sampled among all of them. Generally, models only give access to the average partial width Γ p (I → F, α) associated to a transition from an initial set of levels ( Finally, the partial width to go from I to F is where p is the transition energy, and with δ(α, J π i , J π f ) accounting for spin and parity selection rules, y(I → F ) the Porter-Thomas factor simulating the transition probability fluctuation described by a [18] and ρ(E f , J π f )dE the number of levels in F . At S n , neutron emission is unlikely. Conversion electrons are taken into account using the BrIcc code [19]. For gamma emission, the average partial width (Γ γ ( γ , XL)) depends on the emitted gamma energy ( γ ), its type X (electric or magnetic), its multipolarity L and the radiative strength function model (f XL ): The E1 transition is described by the Enhanced Generalized Lorentzian Model (EGLO) which is a Lorentzian function where an asymptotic term is added to better reproduce low energy experimental data [15,20]. The other XL transitions are best described by a Standard Lorentzian Model (SLO) [21] defined by where B n is the neutron binding energy and R a nucleus mass dependent ratio. More details can be found in [15]. The Stereo simulation is based on Geant4 libraries [22] and includes the detailed geometry of the detector, the description of its response with special emphasis on light emission and collection [1]. Neutron transportation is handled by the NeutronHP libraries, in which microscopic interaction cross-sections are from the ENDF/B-VII.1 evaluation. Standard deexcitation processes in Geant4 do not offer a satisfactory treatment on an event-by-event basis regarding energy conservation. As a consequence, when a neutron is captured on a Gd isotope, standard NeutronHP processes are bypassed and an user-defined process is used in the Stereo simulation. An empirical gamma-cascade treatment was initially performed using an additional support for the GLG4sim package [23], developed specifically for neutrino detection in LS. This implementation gave satisfactory results in larger detectors for well contained energy depositions. In the new Stereo simulation, the deexcitation cascades from Fifrelin are directly used. In both cases, the deexcitation products are generated isotropically. For natural Gd, the gamma multiplicity per cascade is about 4, and differs only by a few percents between GLG4sim and Fifrelin. The major difference between the codes can be seen in the energy distribution of single deexcitation products, presented in fig. 2. About 15% of the gammas generated in the Fifrelin simulation have an energy higher than 3.5 MeV, while they only account for 7% in the GLG4sim modeling. Recently, independent measurements [6] have shown that these high-energy gammas are needed for an accurate description of the cascade. This is of primary importance for the Stereo detector since at such energies around 5 MeV, the mean free path of a gamma in the LS is 40 cm, comparable to the characteristic size of a Stereo cell. Conversion electrons are present in about 70% of the cascades with a most probable energy of 70 keV. Due to very low emission probability, electrons of more than 200 keV represent less than 1% of the sample.
In the Stereo experiment, the neutron responsecharacterizing the delayed signal-is monitored using an americium-beryllium (AmBe) source deployed regularly in 5 of the 6 identical 91 cm high target cells at 5 different vertical levels (10, 30, 45, 60 and 80 cm from the bottom). Neutrons are produced at a 15 · 10 3 s −1 rate through the reaction: α+ 9 Be → 12 C+n. In about 60% of the cases [24], the neutron emission is accompanied by a 4.4 MeV gamma from carbon deexcitation. A coincidence selection is then applied to isolate the neutrons from these gammas and to get a clean and pure neutron capture sample without background: delayed signals are searched in a time window of 100 μs after a prompt 4.4 MeV gamma signal, and contributions from random coincidences (accidental background) are statistically subtracted. The resulting delayed energy spectrum is presented in fig. 3, where both the 2.2 MeV peak from H(n, γ) and the ∼ 8 MeV from Gd(n, γ) are visible, along with the simulations. The shape of the tail over all the energy range is greatly improved with the Fifrelin description, for central positions, as well as for border positions, more sensitive to escaping gammas. Performing a Kolmogorov-Smirnov test on the tail from 3 to 7 MeV, we find an agreement with a probability of 11%, providing no indication for a systematic effect in the description at the central position, whereas the test showed clear incompatibility between data and the GLG4sim spectrum at more than 5 standard deviations. As expected, the presence of higher energy gammas tends to correct the balance between low and high energy events in the tail. The Stereo energy scale being anchored to the low energy gamma of 54 Mn [1], the mean positions of the reconstructed peaks are artificially higher than literature values, due to quenching effects. These non-linearity effects are calibrated and taken into account in the Stereo simulation. In order to assess the improvement of this new simulation without considerations on any residual linear systematic effect on the energy scale, the reconstructed energy of both simulations is scaled such that the mean position of the H(n, γ) peak from the simulation matches the data. In this way, the agreement for the Gd(n, γ) peak is evaluated relatively to the H capture peak at 2.2 MeV, and an agreement for the Gd mean peak position at the sub-percent level is achieved with the Fifrelin simulation.
The neutrino detection uncertainty in Stereo is dominated by the systematic uncertainty on the delayed neutron detection efficiency. Beyond selection cuts related to the event topology [2], the delayed event of an IBD candidate is required in this article to be within a time window of (0.5-70) μs after the prompt signal. The largest inefficiency is coming from the (4.5-10) MeV energy cut on the delayed event set in order to select only Gd events. The lower threshold is chosen to maximize the signal-tobackground ratio and to minimize the systematic uncertainties, and it is clear from fig. 3 that a significant part of the Gd events -with large energy leakage-is not included. The correct description of the spectrum over the full energy range is therefore essential.
To quantify the impact of the event selection cuts, the neutron capture spectra are divided in two parts: a H window (1.5-3) MeV and a Gd window (3)(4)(5)(6)(7)(8)(9)(10) MeV. The ratio of events in the Gd window (N Gd ) and the total sum (N Gd + N H ) is defined as the Gd-fraction (ε Gd ) The impact on the delayed selection cuts (time and energy) are evaluated by defining the IBD efficiency ε IBD , fraction of events in the (N Gd ) passing the tighter delayed selection used in the neutrino analysis: ε IBD = N (4.5-10) MeV and (0.5 μs <ΔT <70 μs) N Gd .
The total delayed detection efficiency ε tot is then the product of these two terms: The numbers in table 2 illustrate that the Stereo data favor the Gd spectrum from the Fifrelin events as compared to the GLG4sim simulation. Using GLG4sim, the ratio R = ε Data tot /ε MC tot quantifying the agreement between data and MC for the total efficiency was found to be 0.9537 at the most central calibration point as shown in table 2. In the new simulation, ε MC tot matches the data within 1% (R = 0.9953). For ε Gd , small discrepancies remain, mainly in the border positions. The data/MC ratio of ε Gd is very sensitive to the treatment of the neutron propagation. In particular the modeling of neutron scattering and thermal diffusion in the detector as well as the neutron detection cross-section could induce an additional mismatch. The border regions are more sensitive to such inaccuracies. Therefore, as an extreme case, the calibration data at the top of cell 1 was investigated, for which the source was located only 12 cm (8 cm) from the target wall (cell top). Overall very good data/MC agreement is achieved for ε IBD due to the improvements related to the new simulation input of Fifrelin (see fig. 3).
In summary, the correct description of the deexcitation process of the Gd nuclei after neutron capture is crucial for neutrino detection experiments using Gd in general. In particular, this is the case for small detectors sensitive to gamma escape such as Stereo. Since nuclear level schemes are not completely experimentally known, nuclear models are needed. The Stereo description of the Gd cascade was greatly improved using the Fifrelin nuclear code, making benefit of the most updated nuclear databases and user feedback on nuclear evaluations.
Open access funding provided by Max Planck Society. This work is supported by the French National Research Agency (ANR) within the Project No. ANR-13-BS05-0007 and the programs P2IO LabEx (ANR-10-LABX-0038) and ENIGMASS LabEx (ANR-11-LABX-0012). We acknowledge the support of the CEA, CNRS/IN2P3, the ILL and the Max-Planck-Gesellschaft.

Data Availability Statement
This manuscript has associated data in a data repository. [Authors' comment: We make available ten millions of deexcitation cascades for each isotope at https://doi.org/10.5281/zenodo.2653786, since other running and upcoming projects might profit from these data as well.] Publisher's Note The EPJ Publishers remain neutral with regard to jurisdictional claims in published maps and institutional affiliations.
Open Access This is an open access article distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/4.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.