Gamma Ray Spectrum from Thermal Neutron Capture on Gadolinium-157

We have measured the $\gamma$-ray energy spectrum from the thermal neutron capture, ${}^{157}$Gd$(n,\gamma){}^{158}$Gd, on an enriched $^{157}$Gd target (Gd$_{2}$O$_{3}$) in the energy range from 0.11 MeV up to about 8 MeV. The target was placed inside the germanium spectrometer of the ANNRI detector at J-PARC and exposed to a neutron beam from the Japan Spallation Neutron Source (JSNS). Radioactive sources ($^{60}$Co, $^{137}$Cs, and $^{152}$Eu) and the reaction $^{35}$Cl($n$,$\gamma$) were used to determine the spectrometer's detection efficiency for $\gamma$ rays at energies from 0.3 to 8.5 MeV. Using a Geant4-based Monte Carlo simulation of the detector and based on our data, we have developed a model to describe the $\gamma$-ray spectrum from the thermal ${}^{157}$Gd($n$,$\gamma$) reaction. While we include the strength information of 15 prominent peaks above 5 MeV and associated peaks below 1.6 MeV from our data directly into the model, we rely on the theoretical inputs of nuclear level density and the photon strength function of ${}^{158}$Gd to describe the continuum $\gamma$-ray spectrum from the ${}^{157}$Gd($n$,$\gamma$) reaction. Our model combines these two components. The results of the comparison between the observed $\gamma$-ray spectra from the reaction and the model are reported in detail.

We have measured the γ-ray energy spectrum from the thermal neutron capture, 157 Gd(n, γ) 158 Gd, on an enriched 157 Gd target (Gd 2 O 3 ) in the energy range from 0. 11 MeV up to about 8 MeV. The target was placed inside the germanium spectrometer of the ANNRI detector at J-PARC and exposed to a neutron beam from the Japan Spallation Neutron Source (JSNS). Radioactive sources ( 60 Co, 137 Cs, and 152 Eu) and the reaction 35 Cl(n,γ) were used to determine the spectrometer's detection efficiency for γ rays at energies from 0.3 to 8.5 MeV. Using a Geant4-based Monte Carlo simulation of the detector and based on our data, we have developed a model to describe the γ-ray spectrum from the thermal 157 Gd(n,γ) reaction. While we include the strength information of 15 prominent peaks above 5 MeV and associated peaks below 1.6 MeV from our data directly into the model, we rely on the theoretical inputs of nuclear level density and the photon strength function of 158 Gd to describe the continuum γ-ray spectrum from the 157 Gd(n,γ) reaction. Our model combines these two components. The results of the comparison between the observed γ-ray spectra from the reaction and the model are reported in detail.

Introduction
Gadolinium, A 64 Gd, is a rare earth element. Its natural composition ( nat Gd) includes isotopes with the atomic mass numbers A = 152,154-158 and 160. The element features the largest capture cross-section for thermal neutrons among all stable elements: ∼ 49000 b. This is due to the contributions of the isotopes 155 Gd (60900 b [46]) and especially 157 Gd (254000 b [46]).
In nuclear physics, gadolinium isotopes have been studied in neutron capture γ-ray spectroscopy and photoabsorption measurements to obtain information on its nuclear structure and properties [10,15,20,21,29,31,41,42,44,54,[58][59][60]. The spectroscopic (n, γ) measurements allow to catalogue neutron capture resonances and to probe the high density of nuclear energy levels around the neutron separation energy S n in the product nucleus A+1 Gd. Moreover, they allow to identify discrete nuclear states between the ground state and S n of A+1 Gd. Together with the inverse reaction (γ, n) in photoabsorption measurements, the neutron capture γ-ray spectroscopy allows to determine the nuclear level density and the photon strength function of Gd.
Recently, natural gadolinium also plays a role in experimental neutrino physics through the identification of the electron anti-neutrino (ν e ) interactions. The presence of gadolinium enhances the tagging of the neutron produced from the inverse beta decay (IBD) reaction of the MeVν e on a free proton: ν e + p → n + e + . Until now, the element has been used as neutron absorber only by scintillator-based detectors [1,3,8,49]. However, the addition of gadolinium to water Cherenkov neutrino detectors is studied and will soon be applied on a large scale in Super-Kamiokande (SK) [55,61].
One important property of the A Gd(n, γ) reaction is that the deexcitation of the compound nucleus A+1 Gd * proceeds not necessarily by one but by a cascade of on average four γ-ray emissions [21]. Due to the Cherenkov threshold, the variable number of γ rays and their energy distributions within the cascades effectively decreases the mean visible energy release from the neutron capture to below the Q-value. As a consequence, a reliable assessment of neutron tagging efficiencies in Cherenkov detectors with the help of Monte Carlo (MC) simulations strongly depends on a precise model for the full γ-ray energy spectrum from the thermal Gd(n, γ) reaction. More seriously, such a model is important for non-hermetic ν e monitors [49], where an accurate assessment of their neutron detection efficiency strongly depends on a precise model for the γ-ray energy spectrum from Gd(n, γ).
There have been several publications on measured γ-ray spectra from Gd(n, γ) reactions for neutron energies ranging from meV to MeV [4,15,31,58,60]. Recently, the Detector for Advanced Neutron Capture Experiments (DANCE) at the Los Alamos Neutron Science Center (LANSCE) has extensively studied the γ-ray energy spectra from the radiative neutron capture reaction at various multiplicities in the neutron kinetic energy range from 1 to 300 eV for 152,154,155,156,157,158 Gd targets [10,21,42]. Their comparison of the data to MC simulations with the DICEBOX package [14] showed fair agreement. There are some publications [4,30] measuring prompt prominent γ rays with limited acceptance, but there 2/27 have been few measurements of the prompt γ rays covering almost the full spectrum from 0.1 MeV to 9 MeV from the capture reaction on 157 Gd at thermal neutron energies, which enable us to compare them with the modeling in Monte Carlo simulation.
In the following, we report on a measurement of the γ-ray energy spectrum from the radiative thermal neutron capture on an enriched 157 Gd sample with excellent γ-ray energy resolution, high statistics and low background. It was performed with the germanium (Ge) spectrometer of the Accurate Neutron-Nucleus Reaction Measurement Instrument (ANNRI) [32,[37][38][39][40] that was driven by a pulsed neutron beam from the Japan Spallation Neutron Source (JSNS) at the Material and Life Science Experimental Facility (MLF) of the Japan Proton Accelerator Research Complex (J-PARC) [47]. Using the time-of-flight (TOF) method, capture reactions of neutrons in the energy range from 4 to 100 meV could be accurately selected for the analysis. The obtained data covers the entire spectrum from 0.11 MeV to about 8 MeV with observed γ-ray multiplicities one to three. Based on our data and a Geant4 [2,6] detector simulation of our setup, we have developed a model to generate the full γ-ray spectrum from the thermal 157 Gd(n,γ) reaction. This constitutes an important step towards a corresponding model for the nat Gd(n,γ) reaction, which is ultimately relevant for neutrino detectors with gadolinium loading.

Physics Motivation
It is a common technique for ν e detection in the MeV regime to search for the delayed coincidence signals from the products of the IBD reaction ν e + p → n + e + , which has a threshold energy of about 1.8 MeV [51,52]: The "prompt signal" occurs a few nanoseconds after the interaction and originates from the energy loss and the annihilation of the emitted positron. At low energies, when the invisible recoil energy of the neutron can be neglected, one can reconstruct the ν e energy from the prompt event's visible energy E prompt as E ν = E prompt + 0.782 MeV [12].
The "delayed signal" stems from the γ-ray emission following the capture of the thermalized neutron on a nucleus of the detector's neutrino target material. Neutrons produced by neutrinos in the MeV regime via the IBD reaction typically have kinetic energies up to several tens of keV and interact between ten to twenty times via elastic scattering with hydrogen before they are thermalized [7,25]. The mean timescale τ cap for the neutron capture depends on the concentrations n i and the thermal neutron capture cross-sections σ cap,i of the nuclei i in the detector material as well as on the mean velocity v n of the produced neutrons: τ cap ∝ 1/(n i σ cap,i v n ). With hydrogen, carbon and oxygen nuclei naturally being present in common low-energy ν e detectors, e.g., organic liquid scintillator and water Cherenkov detectors, the mean neutron capture time is usually on the order of a few tens to hundreds of microseconds. Table 1 summarizes thermal neutron capture cross-sections and Q-values for the most abundant isotopes of these elements.
Recently, it has become a common technique to add a mass fraction of 0.1-0.2% of gadolinium into the neutrino targets of organic liquid scintillator [1,3,8] and water Cherenkov [11,22,55,61] detectors in order to enhance the neutron tagging efficiency for IBD events. This basic technique was first demonstrated in the discovery of the neutrinos with a cadmium-loaded liquid scintillator in 1956 [51,52]. On the multi-kiloton scale, ν e Table 1 Cross-sections [46] and Q-values [19] Table 2 Relative abundances of gadolinium isotopes in natural gadolinium [53] and their radiative thermal neutron capture cross-sections [46]. detection with gadolinium-enhanced neutron tagging will first be done by SK. A corresponding project, SK-Gd, will start soon, after EGADS successfully demonstrated the sustainable gadolinium loading of water [55,61]. The demonstrated feasibility to load common neutrino target materials with gadolinium is based on two positive properties: the large capture cross-section for thermal neutrons, especially of 157 Gd, and the high Q-value, 7937 keV [19] for 157 Gd(n, γ), compared to the values listed in Table 1. The reason for the large cross-section of the gadolinium isotope is an s-wave neutron capture resonance state in the thermal energy region with a resonance energy of 31.4 meV for 157 Gd [50]. A list with the thermal neutron capture cross-sections of all the gadolinium isotopes in natural gadolinium, which defines the composition of how gadolinium is commonly loaded to neutrino target materials, is given in Table 2. 1 The about 8 MeV excitation energy from the Gd(n, γ) reaction is released in several γ rays. Due to the calorimetric measurement, liquid scintillator detectors simply need to look for this energy deposition, assuming that all the γ rays are fully contained inside the active volume. A water Cherenkov detector, however, detects only a part of it due to the above-mentioned energy threshold. Therefore, good understanding of the multiplicities of γ rays from Gd(n, γ) reactions and their energy distributions in the range 0.1-8 MeV is an important prerequisite to properly predict neutron tagging efficiencies in gadolinium-doped water Cherenkov detectors based on MC simulations.

Experiment
We performed our measurements of the thermal neutron capture on gadolinium with an enriched 157 Gd target inside the Ge spectrometer of ANNRI [32,[37][38][39][40] at JSNS of J-PARC in December 2014. The JSNS complex provides neutrons with energies up to 100 keV. Its beam is one of the most intense pulsed neutron beams for precise neutron TOF experiments in the world, especially in the thermal energy region. The ANNRI detector, located at Beam Line No. 4 [32] of the MLF, is dedicated to measure cross-sections and γ-ray spectra of neutron-nucleus interactions with excellent energy resolution compared to other γ-ray spectrometers.

Detector Setup
During our measurements, the JSNS was powered by a 300 kW beam of 3 GeV protons in "double-bunch mode" that hit a mercury target at a repetition rate of 25 Hz. This created a double of 100 ns wide neutron beam bunches with 600 ns spacing every 40 ms. At the target position inside the ANNRI spectrometer, which is located 21.5 m from the neutron beam source, the neutron beam delivered an energy-integrated neutron intensity of about 1.5 × 10 7 /cm 2 /s.
The ANNRI spectrometer consists of two Ge cluster detectors with anti-coincidence shields made of bismuth Ge oxide (BGO) and eight co-axial Ge detectors. Since the co-axial detectors were still in repair after the Tohoku earthquake on 11 March 2011, we only used the two Ge cluster detectors shown in Fig. 1(a) in the present analysis. The clusters are placed perpendicular to the aluminum beam pipe ( Fig. 1(b)), with the front faces 13.4 cm above and below the target position. They provide a combined solid angle coverage of 22% with respect to this point. As shown in Fig. 1(c) each of the 7 crystals in the cluster has its hexagonal surface facing the target. The dimensions of a Ge crystal are shown in Fig. 1(d) and (e).
The BGO anti-coincidence shield for one Ge cluster (see Fig. 1(a)) consists of a cylindrical BGO counter, which is separated into twenty readout blocks: twelve around a cluster and eight covering its rear side. The shields provide a total solid angle coverage of 55% with respect to the target.
In order to reduce background γ rays from the neutron capture by the aluminium layer on the beam pipe, the inner face of the pipe is lined with a layer of lithium fluoride of ∼1 cm thickness. Moreover, shields made of lithium fluoride and lithium hydride are located between the pipe and the Ge clusters to protect the crystals from the impinging neutrons. The remaining γ-ray background was measured directly by placing only the empty target holder, whose dimensions are shown in Fig. 1(f), inside the neutron beam.

Data Acquisition
The data acquisition (DAQ) system [36] was triggered when at least one of the fourteen Ge crystals had a collected charge equivalent of more than 100 keV. All further energy depositions in the crystals within a time window of 560 ns (smaller than the double-bunch spacing) after the trigger were combined with the initial deposition to form an event. Within this event, we only considered crystals with a collected charge corresponding to more than 100 keV as hit. The crystal hits of cluster were accepted if none of the 20 surrounding BGO blocks had an energy deposition greater than 100 keV within the same time window. The data stored per event included the neutron TOF, given by the time difference between the first detected hit of a crystal (trigger time) and a signal from the JSNS, as well as the collected charge (energy deposition) and the hit time delay with respect to the trigger time of every hit crystal.
For the purpose of dead time correction, signals from a random pulse generator with an average rate of 570 Hz were fed into the pre-amplifier of every Ge crystal and simultaneously counted by a fast counter. The amplitudes were set to be about an energy equivalent to 9.5 MeV. The ratio r L,i (=N r,i /N s ) of the number of pulses N r,i recorded by the i-th crystal to the number of pulses N s corrects the absolute elapsed time of the experiment T for the dead time of the crystal's DAQ system after a trigger, giving the crystal's effective live time as T L,i = T · r L,i . On average, r L,i is about 94%. The dead time correction is important for calibration and background subtraction.

Event classification
We assigned a multiplicity value M and a hit value H to each recorded event. We defined the multiplicity M as the combined number of isolated sub-clusters of hit Ge crystals at the upper and the lower clusters. A sub-cluster is formed by the neighboring hit Ge crystals and can be of size ≥ 1. The hit value H describes the total number of Ge crystals hit in the event. The multiplicity M represents the number of γ rays and the hit value H represents the lateral spread of γ rays. Figure 2 shows some examples (right) together with the numbering scheme used to reference individual Ge crystals (left).
Since we assume that M is the number of detected γ rays, this implies that sub-clusters with sizes greater than one are mainly due to scattering of one γ-ray between neighboring crystals and not due to multiple γ rays.

Detector simulation
Based on the geometry and material specifications for ANNRI (see Fig. 1), we have developed a detailed detector simulation using version 9.6 patch 04 of the Geant4 toolkit. It uses the Geant4 implementation G4EmPenelopePhysics of physics models for low-energy photon, electron-positron interactions developed for the PENELOPE (PENetration and Energy LOss of Positrons and Electrons) code version 2001 [26].
With the MC simulation we evaluated the detector response to the simultaneous propagation of one or more γ rays with specified energies through the setup. During the simulation of an event, each Ge crystal accumulated the energy depositions from charged particles. This information was then used to realize the trigger and veto scheme described in Sect. 3.2. We validated the MC simulation by comparing its outcomes for the fractions of different event classes (Sect. 3.3) and the energy spectra observed by the single crystals to the data taken with calibration sources (Sect. 3.5): Radioactive 60 Co dominantly emits two γ rays, 1173 keV and 1332 keV, after its β − decay to 60 Ni. We used the lower cluster of ANNRI to tag one of the γ rays in a single crystal and looked at the crystal hit configuration created by the other γ ray in the upper cluster. This resulting hit configuration was classified with multiplicities 7/27 To study the hit configurations at higher energy, we used the single γ ray of 8579 keV from the thermal neutron capture 35 Cl(n, γ) reaction, which is produced via direct M1/E2 transition from 8579 keV to the ground state (2 + → 2 + ) [46]. Table 3 summarizes the fractions of the different event classes created by the γ rays of different energies in our experimental data and our MC simulation. We only selected events with M = 1, i.e., with one sub-cluster of hit crystals. Using the MC simulation, we estimated the background contribution which comes mainly from 6 prominent two-step deexcitation γ rays from 8579 keV using the MC simulation [24] to be about 9% for M1H2 case and 24% for M1H3 case. The table lists the values after subtracting these contributions. The systematic errors to the numbers given in the table due to this overlap effect are negligible.
As one can see from Table 3, the agreement between data and MC for the two 60 Co lines is very good. Despite the errors for the experimental data on the 8579 keV line from 36 Cl due to the above corrections, the agreement with the MC simulation is also good. Overall, the summary shows that energy migration to neighboring and distant crystals within a Ge cluster, which increases with increasing γ-ray energy and arises from Compton scattering of the γ ray or γ-induced e + e − pair production, is correctly reproduced within our MC simulation.
Moreover, Fig. 3 shows the energy spectra for 60 Co (left) and 137 Cs (right) from M1H1 events observed in our calibration data and corresponding MC. One can see that, in addition to the multiplicities, also the spectral shapes are very well reproduced by our detector simulation.

Background and calibration data
In order to measure the background for the experiment, which originates mostly from γ rays from the interactions of the beam neutrons with materials other than the target, we placed the empty target holder for 6 hours into the neutron beam. Figure 4 shows the background energy spectrum observed by one of the crystals for M1H1 events together with the spectrum observed in the measurement with the enriched 157 Gd sample before background subtraction. The background spectrum, after processing in the same way as the data and the live-time normalization, contributes only ∼0.06% to the gadolinium data spectrum. The energy calibration of the ANNRI Ge crystals was done with known γ-ray lines from the radioactive sources 60 Co, 137 Cs and 152 Eu as well as from the deexcitation of 36 Cl after the thermal 35 Cl(n,γ) reaction in a sodium chloride (NaCl) target. Table 4 summarizes the measurement time and number of observed events for the different sources and targets.
The energy resolutions (σ(E)) of all the 14 crystals were measured over the energy from 0.3 to 8 MeV and they are expressed as σ(E)(keV) = 1.8 + 0.00041E(keV). With the known activities β of our 60 Co, 137 Cs and 152 Eu sources, we estimated absolute single photopeak efficiencies ε i (E γ ) at different energies E γ for each crystal (i) as

9/27
where N i (E γ ) is the number of detected γ rays in the ±3σ region of a Gaussian fitted to the photopeak observed by the i-th crystal at energy E γ , BR γ is the branching ratio for the decay branch emitting the γ ray of energy E γ , and T L,i is the corrected livetime. The single photopeak efficiency values at various energies from the measurements with the radioactive sources and the NaCl target cover the range from 344 keV to 8579 keV for each crystal. The values for one of the crystals are depicted in Fig. 5. The relative efficiency values for the NaCl target were normalized with respect to the dominant 7414 keV line, which itself was normalized with our MC simulation.
The corresponding prediction for each crystal (i) was calculated using the MC simulation as where N i (E γ ) is the number of γ rays and N (E γ ) denotes the total number of generated γ rays with energy E γ . The data points and the MC simulation are in good agreement. The data from the 60 Co and 137 Cs calibration sources also allowed us to check the uniformity of the detector. For this purpose, we compared the nominal value of the source's activity to the value measured by each Ge crystal. The ratios of data to nominal value are shown in Fig. 6.
Taking the error bars into account, the spread of the single 60 Co ( 137 Cs) ratios with respect to the mean ratio shows a uniformity of the detector response over the solid angle of the detector at the 8% (14%) level. In other words, the detection efficiency is well understood over all crystals at this level of uniformity. Further details are described in Appendix A.

Gadolinium data
For the measurement with gadolinium we attached the enriched 157 Gd target in the form of gadolinium oxide powder (Gd 2 O 3 ) in a teflon sheet to the designated holder within the neutron beam line at the center of the ANNRI detector. Taking into account the isotopic composition of the commercial gadolinium sample (Table 5) and the dominant cross-section of 157 Gd (see Table 2), the target is essentially a pure 157 Gd target for thermal neutrons. A total of 1.81 × 10 9 events were collected with this target in about 44 hours of data taking.    where m n is the neutron mass and L is the 21.5 m distance between neutron source and target. The resulting neutron energy spectrum is shown in Fig. 7. In order to study the γ-ray spectrum solely from thermal neutron capture on 157 Gd, we only selected events from neutrons in the kinetic energy range [4,100] meV for the present analysis. After the neutron energy selection and the subtraction of the background, the resulting event sample was divided into sub-samples based on the multiplicity M and hit value H of the events. Figure 8 shows the energy spectra observed by the crystal 6 for different multiplicity values M (M1H1, M2H2 and M3H3). We mainly show the spectra from the hit configurations M1H1, M2H2 and M3H3, since they are the majority of the events among each multiplicity value (M=1, 2 and 3) and they are less subject to the overlap with multiple γ rays.
The observed spectra are consistent within about 7% for the dominant M1H1 events for all 14 detectors. The observed energy spectra are dominated by the γ rays from the thermal 12/27 157 Gd(n, γ) reaction, especially when we selected M1H1 and M2H2 events, since a clean single hit on one crystal suppresses the effect of Compton scattering. At low energy, the spectra are slightly distorted by the effect of the Compton scattering.
The M1H1 spectra in Fig. 8 exhibit two components: discrete peaks, very well visible below 1.6 MeV and above 4.8 MeV, 2 and a continuum, most prominent between the previous energy regions. The origins and features of these components and how we implemented them in our spectrum model will be discussed in the following sections.
4. Gamma-rays from thermal 157 Gd(n, γ) reaction: Emission scheme and model Our approach to model the γ-ray spectrum is a separate description of the continuum component and the discrete peaks visible in Fig. 8. We followed the strategy of the GLG4sim package [57] for Geant4 to which we compare our results in Sect. 5.

Emission scheme
After the thermal neutron capture on 157 Gd, the remaining 158 Gd * compound nucleus is in an s-wave neutron capture resonance state with an excitation energy of 7937 keV and spin-parity J π = 2 − [46]. It deexcites via a cascade of on average four γ-ray emissions [21] to the ground state of 158 Gd with J π = 0 + .
As illustrated on either side of Fig. 9, the density of nuclear levels increases with increasing excitation energy from the domain of well separated (discrete) levels, where spin and parity of the states are known, to a quasicontinuum where individual states and energy levels cannot be resolved. Since the two regions are connected smoothly, there is no obvious, sharp boundary between them. For the purpose of modeling, an arbitrary transition point is commonly defined at an excitation energy up to which supposedly complete information on discrete levels is available, e.g., 2.1 MeV in Ref. [21].
As depicted on the left of Fig. 9, the continuum component of the γ-ray spectrum from the thermal 157 Gd(n, γ) reaction stems from multi-step deexcitations of 158 Gd * . Such intermediate transitions from the neutron capture state down towards the ground state can occur between (unresolvable) levels in the quasicontinuum (dashed lines), within the domain of discrete levels (solid lines) and between two levels from each of these smoothly connected regions. Both the number and energy values of the emitted γ rays (i.e., the intermediate levels) are random.
The discrete peaks on top of the continuum mainly originate from the transition from the neutron capture state to the low-lying levels as illustrated on the right of Fig. 9 and they are studied in the previous publications [4,15]. While the discrete peaks in our model are based on the previous publications and their intensities are adjusted to agree with our own data, we employ a statistical approach to describe the continuum component in the γ-ray energy spectrum of 158 Gd * , which dominates with a contribution of (93.06 ± 0.01)% to our data. The approach is to follow Fermi's Golden Rule [25], which states that the transition probability per unit time is proportional to the product of the transition matrix element squared between the initial and the final states and the state density at the final state: Left: Illustration of a multi-step γ-ray cascade from the neutron capture (n-capture) state down towards the ground state via many intermediate levels in the deexcitation of 158 Gd * after the thermal 157 Gd(n, γ) reaction. Right: Example for a two-step cascade that proceeds via a low-lying level and includes the emission of a high-energy γ ray. This cascade type contributes to the creation of discrete peaks in the high-energy part of the γ-ray spectrum.
Starting from an excited state with energy E a , the differential probability dP (E a , E b )/dE that the nucleus undergoes a transition to a state with energy E b (< E a ) and emits a γ ray of multipolarity XL (X = E, M and L = 1, 2, . . . for electric, magnetic and angular momentum, respectively) with energy E γ = E a − E b is expressed as The first factor, ρ(E b ), is the nuclear level density (NLD) at the final state (E b ) [18]. The second factor is the sum over the transmission coefficients T XL (E γ ) for the different multipolarities, each of them depending on the corresponding photon strength function (PSF), f XL (E γ ), for electromagnetic decay.
Since ρ(E b ) increases exponentially as E b increases, this factor favors transitions from E a to E b which is large and close to E a , and thus favors small The PSF describes the coupling of a photon with given energy and multipolarity to the excited nucleus. Electric dipole (E1), magnetic dipole (M1) and electric quadrupole (E2) radiation are the most relevant multipolaries [18]. The experimental photonuclear data [13,23] show that the photoabsorption cross section of statically deformed spheroidal nuclei like 158 Gd is well approximated by that of the Giant Dipole Resonance (GDR) as the superposition of two Lorenzian lines, corresponding to oscilations along each of the axes of the spheroid. The simplest model for PSF is thus called the Standard Lorentzian Model (SLO) having two Lorentzian forms [9,18]: where usually two sets (i = 2) of parameters are used to describe the two GDRs in terms of a resonance energy E i (in MeV), resonance width Γ i (in MeV) and a peak cross-section 14  Gd reaction and its composition in terms of contributions from the first (red), second (blue), third (orange), fourth (green) γ ray and other γ rays (gray). The distributions are normalized such that the total continuum spectrum actually is a binned probability distribution and that the relative contributions of the single components are properly reflected. Right: γ-ray multiplicity distribution obtained from the continuum part of our spectrum model. Five million events were generated for the distributions.
σ i (in mb). As shown in Fig. 12, this factor favors large E γ in the energy regime < 8 MeV and the factor E 2L+1 γ favors large E γ as well. As a result of two competing factors of the nuclear level density and the transmission coefficients, we obtain the broad peak structure in the continuum spectrum distribution of each γ ray. The distribution P (E a , E b ) for the first, second, third and the fourth γ ray (as in Fig. 9-left) is shown in Fig. 10. Their combination for the total continuum spectrum is shown as well. The dips in the distributions at 0.4 MeV and 7.5 MeV are due to a corresponding feature in the NLD model we employ. Fig.10-left shows a dominance of the continuum component above 5 MeV by the first γ ray. This naturally suggests that the discrete peaks above 5 MeV are generated from the first transition [4,15].
For completeness, Fig.10-right shows the γ-ray multiplicity distribution as generated by the continuum part of our spectrum model. Some more remarks on this distribution are given in Sect. 4.2.1.

The MC Model ("ANNRI-Gd model")
In line with the GLG4sim approach [57], our model for the γ-ray spectrum from the radiative thermal neutron capture reaction 157 Gd(n, γ) consists of two separate parts: the discrete peaks contribute (6.94 ± 0.01)%, while the continuum component dominates the remaining (93.06 ± 0.01)% of our data.
The model is written in C++ and is used through the program structure of our Geant4based detector simulation. Gd level density (HFB) 158 Fig. 11 Tabulated values [18,33] for the NLD of 158 Gd from computations based on the HFB method [27,28]. We used linear interpolation between the points, which have a spacing of 0.25 MeV (0.5 MeV) below (above) 5 MeV excitation energy.
E a − E b , we complete Eq. (4) with a proper normalization as where δE is a finite energy step in our computations. The E1 transmission coefficient is Note that, due to the normalization in Eq. (6), the absolute values of the NLD and the PSF do not matter for our purpose. The detailed comparison of our model with our data is presented in Sect. 5. Fig. 10-left shows the binned probability distribution P (E a , E b ) of the γ-ray energies generated by our continuum model part. It also shows the contributions of γ rays from different transition steps in a cascade. In Fig. 10-right one can see the corresponding γ-ray multiplicity distribution. The probability distributions P (E a , E b ) are tabulated for each energy E a for the continuum part. To prevent the generation of infinite cascades or small, negative γ-ray energies, we artificially force a cascade to end after a finite number of steps: If the remaining excitation E falls below a threshold value of E thr = 0.2 MeV, one last γ ray with low energy E is generated. Therefore, one "additional" γ ray is produced per cascade, effectively increasing the mean value by one. This procedure assures the total energy conservation.

Nuclear level density.
To describe the NLD of 158 Gd as a function of excitation energy, we used a microscopic combinatorial level density computed according to the Hatree-Fock-Bogoliubov (HFB) method [27,28]. Tabulated values are provided separately for positive and negative parity levels [18,33]. We point-wise summed the two values and used linear interpolation in our calculations. For the modeling of the continuum component we used the HFB model results over the entire excitation energy range (Fig. 11). 16 Fig. 12 The E1 PSF for 158 Gd, given as a function of the γ-ray energy, that we used in our model for the spectral continuum component from the thermal 157 Gd(n, γ) reaction. It follows the SLO model in Eq. (5) with the parameters given in Table 6.

E1 photon strength function.
We list four sets of values for the GDR parameters E i , Γ i and σ i in Table 6 [56]. For our present model, only the first two E1 PSFs are used. Fig. 12 shows the resulting PSF shape. Two prominent GDR peaks i = 1, 2 are clearly visible.

4.2.2.
Discrete peaks. The previously described model for the continuum component of the γ-ray spectrum from the thermal 157 Gd(n, γ) relies on a continuous NLD description and thus does not reproduce sharp γ-ray energy lines in the observed spectra. We separately added this spectral component on top of the continuum part and tuned the strengths of the discrete peaks to match our measurement results.
Using our data from all of the ANNRI detector's Ge crystals and selecting the dominant M1H1, M2H2 and M3H3 events, we identified 15 known [48] discrete γ-ray lines above 5 MeV with high intensity after careful exclusion of single and double escape peaks. Following the previous assumption that the peaks in the high-energy part of the spectrum arise from the first transition, we refer to the corresponding γ rays as 'primary' γ rays. We also identified γ rays from subsequent transitions ('secondary' γ rays) in detected multi-γ events (M > 1) by looking at observed γ-ray energies besides the primary one used to tag the event. Two examples for the primary γ rays (5903 keV and 6750 keV) are shown in Fig. 13. The energies of the identified primary and secondary transition γ rays as well as their relative intensities as obtained from our data are listed in Table 7. Note that the direct transition from the 17/27 neutron capture state (J π = 2 − ) to the 158 Gd ground state (J π = 0 + ) is much suppressed because it is of M2 type. The energies in the table were not determined from our data but taken from [48]. For cases where peaks obviously overlapped and could not be disentangled, we treated them combined in the intensity evaluation and list the mean primary γ-ray energy. The secondary γ rays in Table 7 were used as tag to identify further parts of the corresponding decay branches with information from Ref. [48].
A comparison between the mean intensities from our data and values documented in Ref. [48] is shown in Fig. 14 for primary γ rays (left) and the secondary γ rays (right) listed in Table 7.
By summing all the relative intensities in Table 7 we get (6.94 ± 0.01)% as an estimate for the fraction of neutron capture events that contribute to the discrete peaks. The remaining part, (93.06 ± 0.01)%, contributes to the continuum part of the γ-ray spectrum. We note that although our mean values for the discrete primary (secondary) γ-ray intensities in Table  7 are lower by about 20% than the values in the literature, the effect of this difference on the total spectra is very small since the contribution of the discrete component to the total spectra is less than 7%.
In order to implement the identified discrete peaks into our model, we converted the listed mean intensity values to probabilities summing to one and hard-coded them together with the γ-ray energies of the different cascades into our γ-ray generator. A particular cascade from Table 7 is then generated according to this probability distribution. We finally obtain MC γ-ray spectrum as the sum of the continuum and discrete peaks as shown in Fig. 15.

Model performance
To assess the performance of our model for the γ-ray spectrum from the thermal 157 Gd(n, γ) reaction, we compared its output to the measured data. With our model, we separately simulated 2 × 10 8 n-capture events for the discrete peaks and the continuum component, and then merged in corresponding proportion. The resulting spectrum is shown in Fig. 16, along with the data. Table 7 List of the 15 discrete peaks from primary γ rays we identified in our data. The stated energies are taken from Ref. [48], rounded to nearest keV. In four cases the table lists the unweighted mean energy of known peaks that overlap in our data: (i) 6001 keV combining 5995 keV and 6006 keV, (ii) 5669 keV combining 5661 keV and 5677 keV, (iii) 5595 keV combining 5582 keV, 5593 keV and 5.610 keV, as well as (iv) 5167 keV combining 5155 keV and 5199 keV. (Intensities listed w.r.t. total data) On the left side of Fig. 16 one can see the resulting energy spectra observed by one of the crystal for M1H1 events in our data and the simulations. The same for GLG4sim simulations is also shown. Both simulated spectra were normalized such that the total sum of counts in each spectrum's range from 0.8 to 8 MeV is equal to the corresponding sum in the data spectrum. The low-energy limit excludes strong shape deviations below 0.8 MeV due to an excess of low-energy γ rays in data. The shape of the energy spectrum in our data is significantly better reproduced by our model, as seen from the ratios "Data/MC" on the right of Fig. 16. For the presented spectrum with the 200 keV binning, the mean deviation of the single ratios from the mean ratio is about 17% for our model. We checked that the 19 Table 7 to the corresponding values listed in the CapGam data base [50] and ENSDF [48].  Fig. 15 The MC-generated total energy spectrum from thermal 157 Gd(n, γ) for the peripheral crystal 6 (black), and the comprising continuum (red) and discrete (blue) component of the MC energy spectra shown separately. results for our model shown in Fig. 16 are consistent among the 14 Ge crystals for different event classes. We also simulated the γ-ray spectrum with the generic photon evaporation model [5] from the standard tools of the Geant4 framework. Resulting shape deviations with respect to the data were found to be larger than for the GLG4sim model.
The comparison between our data and our model for M1H1 events in Fig. 16 shows some systematic discrepancies: From 0 to 1.2 MeV and from 2.2 to about 6 MeV the model underpredicts the number of γ rays. Between 1.2 and 2.2 MeV and above 6 MeV it makes overpredictions. The shape description of these regions is interconnected because of the conservation of the energy for every capture reaction.
Overall, our model significantly improves the description of the shape of the γ-ray energy spectrum from the thermal neutron capture on 157 Gd with respect to other available models usable with the Geant4 structure. A second aspect to evaluate is the detected γ-ray multiplicity. Our experiment allows to classify events in terms of a multiplicity M of sub-clusters of hit crystals and the total number H of hit crystals. Table 8 lists the fractions of events from the classes M = 1, 2, 3, 4, M ≤ H ≤ 4 in our data and in the MC sample from our model.
The frequency of the most dominant event classes, M1H1 and M2H2, is well reproduced within about 3%; the total agreement between MC and data for M = 1 is about 2%. For event classes with MxHx, x ≥ 2, which have lower frequency / statistics, the MC tends towards increasing underpredictions.

Conclusion
A good model for the γ-ray energy spectrum from the radiative thermal neutron capture on natural gadolinium is an important prerequisite for MC studies to evaluate the efficiency to tag neutrons from IBD events in gadolinium-enhanced ν e searches. This is especially true for water Cherenkov or segmented detectors, where the energy distribution within a γ-ray cascade from the neutron capture on gadolinium heavily impacts the detectability of this marker signal.
Using the Ge spectrometer of the ANNRI detector at MLF, J-PARC, we performed a measurement of the γ-ray energy spectrum from thermal neutron capture on 157 Gd.
Based on our data and a Geant4 simulation of the ANNRI detector, we have developed a model for the γ-ray energy spectrum from the thermal 157 Gd(n, γ) 158 Gd reaction. This marks an important step towards a model for natural gadolinium, which is of use for gadoliniumenhanced ν e measurements.
While the strength information of 15 discrete peaks above 5 MeV in our data is directly included into the spectrum model, the continuum component is modeled using a statistical approach. We used the Standard Lorentzian PSF model to describe the strength of the photon-nucleus coupling as a function of the γ-ray energy.
The measured spectrum agrees within ∼17% with that of our ANNRI-Gd model at 200 keV binning. We found this outcome to be a significant improvement compared to other spectrum descriptions, e.g., from the standard Geant4 or the GLG4sim package. The completeness of 21/27  primary γ ray of energy E γ , there is a chance that one of the secondary γ rays vetoes the primary γ ray hit by directly going into the BGO shield of the corresponding Ge cluster. This effectively reduces the photopeak efficiency compared to the case where solely the primary γ ray would be emitted. The γ rays from the thermal 35 Cl(n, γ) 36 Cl reaction do not allow the determination of absolute efficiency values since the number of emitted γ rays is unknown. Therefore, we computed efficiency values relative to the photopeak efficiency of the most intense line at 7414 keV among our selected lines. The normalization of the reference efficiency was obtained from our MC simulation.
We corrected the single photopeak efficiency for this trigger effect differently for 36 Cl / 152 Eu and 60 Co. From the complex decay and deexcitation schemes of 36 Cl and 152 Eu we only selected γ rays for the efficiency determination that are dominantly emitted alone or with just one additional γ ray in their particular decay channel: 5517 keV, 7414 keV, 7790 keV and 8579 keV for 36 Cl; 344 keV, 779 keV, 1112 keV and 1408 keV for 152 Eu. Relevant branching ratios can be found in Ref. [17,45]. This selection allowed for an easier estimation of the above described inefficiency in the two γ-ray cases by multiplying the raw photopeak efficiency value for a crystal with the correctition coming from our Geant4 MC simulation. It is calculated from the single photopeak MC efficiency ε MC i (E γ ) for the γ ray of interest with energy E γ and the corresponding single photopeak MC efficiency ε MC i,2γ (E γ ; E γ2 ) obtained when the second γ ray with E γ2 is simultaneously propagated through the detector.
For the 60 Co source, which essentially always emits two γ rays (E 1 = 1173 keV and E 2 = 1332 keV) [16], we determined the corrected single photopeak efficiency directly through a fit: We look at a pair of crystals (i, j), i = j, where each crystal is on a separate cluster of ANNRI. The number of observed M1H1 events where E k (k = 1, 2) is deposited in crystal (i) is N ik with its error σ ik . One expects this value to be N ik = βT r L,i ε ik (1 − C i ) with ε ik ≡ ε i (E k ), the live time T and the dead time correction factor r L,i . The efficiency correction (1 − C i ) is due to the inefficiency described above. For a given pair of crystals and the two γ rays this yields four combinations of crystal and γ ray. Moreover, we look at the coinciding detection of both γ rays in M2H2 events by the crystals (i, j). With E l = E k being the second γ ray, the observed number of coincidence events where E k (E l ) is detected in crystal i (j) is N ikjl . Its error is σ ikjl . The expected value is N ikjl = βT r L,ij ε ik ε jl W (θ ij ). Here, r L,ij is the dead time correction factor for the crystal pair (i, j), which typically is on the order of 90%. The factor W (θ ij ) accounts for the predicted angular correlation [43] of the γ rays from 60 Co with angle θ ij , which is given by the angle of the detector pair (i, j). A second combination, N iljk , simply follows from permuting the γ ray energies. With the in total six observables we minimized the expression 23/27 for 48 crystal pairs (i, j), one was excluded, to fit the four uncorrected single photopeak efficiencies, ε ik , ε il , ε jk and ε jl , and βT for different but fixed values of the constant C. The best agreement between the mean of the fitted values of βT and the nominal value was obtained for C = 0.225. Using this constant, we took the averages of the efficiency values per crystal and energy as final results. With this method we obtained a single photopeak efficiency of (1.3 ± 0.1)% at 1.3 MeV for all 14 Ge crystals combined. Figure A1 depicts the ratios of the single photopeak efficiencies from data and from MC at the single γ ray energies averaged over all 14 crystals (left) and for all 14 crystals averaged over the 11 data points (right). On both plots one can see that the weighted mean values of the ratios deviate by less than 10% from perfect agreement and maximum deviations are about 20%. The weighted sample standard deviation of the ratios for all crystals and data points is about 6%. From this study, we conclude that we understand the photopeak efficiency of each crystal not only over the energy range from 344 keV to 8579 keV but also uniformly over the entire solid angle of the detector and that we can reproduce the responce of each crystal very well by our Geant4 detector simulation. (2)) averaged over the 14 crystals at the fixed γ-ray energies (left) and averaged over the 11 data points for the single crystals (right). Linear interpolation between the points from the simulation was used to determine the MC efficiency at an intermediate energy. The calculations of the weighted mean values and the weighted sample standard deviations (error bars) take the errors of the data points (see Fig. 5) into account. Outer error bars indicate the extreme values of the ratios in the respective samples. The 35 Cl(n, γ) data point at 7414 keV is the reference for the normalization of the other data points of this reaction. It perfectly agrees with a ratio of one since it was normalized with the MC simulation.

B. Double/Triple Gamma Spectrum:
The M2H2 and the M3H3 classified data are the next dominant fraction of the data after the M1H1 discussed above. The corresponding spectra in data also agree with those generated by our model as shown in Fig. B1. With increased multiplicity, the slope of the spectrum grows softer, owing to reduced probability of emitting higher energy gamma rays, as also seen in Fig. 8 and independent checks from simulations as in the Fig. 10.