Nuclear deexcitation simulator for neutrino interactions and nucleon decays of 12 C and 16 O based on TALYS

Nuclear deexcitation associated with neutrino interactions and nucleon decays has a significant impact on a recent key observable for neutrino detectors: neutron multiplicity. However, most existing neutrino Monte Carlo event generators do not consider this process. For a comprehensive prediction of the neutron multiplicity, a nuclear deexcitation simulator based on TALYS named NucDeEx is developed. NucDeEx can be used for neutrino interactions and nucleon decayes for 12 C and 16 O. The author provides the NucDeEx codes and sample codes for applying the NucDeex to widely-used neutrino Monte Carlo event generators, NEUT, NuWro, and GENIE.


I. INTRODUCTION
Neutrino-nucleus interactions are frequently described by the impulse approximation, which divides the process into three parts: Initial nucleon state, primary interaction, and final state interactions (FSI).The initial nucleon state is usually parametrized by the Fermi gas model or spectral functions.The primary interaction is the process via the weak interaction of a neutrino with a single nucleon target, treating the other nucleons as mere spectators.The FSI is a cascade rescattering of particles in the nucleus.The description above is adopted in most of the existing neutrino Monte Carlo event generators.Considering actual interaction processes, a hole or holes are created after the FSI.Then, the residual nucleus frequently remains in the excited state and undergoes deexcitation, emitting particles such as protons, neutrons, α particles, and gamma rays.However, with few exceptions [1], this process was not simulated in the major generators until the recent studies described later.Despite this fact, this description has been used for a long time in various neutrino experiments with much success because the particles emitted by nuclear deexcitation have low energies of a few MeV and, except for gamma rays, are undetectable with most neutrino detectors.
Neutrino detectors can detect neutrons in three ways: Detection of recoil protons, deexcited gamma rays from recoiled nuclei, and gamma rays emitted by neutron capture, as shown in Fig. 1.Since the protons have a high Cherenkov threshold and short trajectories of a few centimeters, Cherenkov and tracker detectors are forced to impose high detection thresholds.Because of their low energy loss density, the gamma rays from deexcitation and neutron capture are generally difficult to detect with tracker detectors.The deexcited gamma rays also have a certain energy threshold for neutrons.On the other hand, gamma rays emitted by neutron capture do not require any detection threshold for neutrons.Neutrons thermalized by proton recoils remain in the liquid scintillator and water Cherenkov detectors for several hundred microseconds.In these detectors, most of them are cap-tured by protons emitting a 2.2 MeV gamma ray.Large liquid scintillator detectors such as KamLAND [2] and JUNO [3] are good at measuring these gamma rays because of their large light yield.On the other hand, detecting the 2.2 MeV gamma ray is challenging due to its small light yield for water Cherenkov detectors such as Super-Kamiokande [4].Ongoing detector updates are overcoming this issue by dissolving gadolinium (Gd) in water, such as Super-Kamiokande Gadolinium [5] and ANNIE experiment [6].Neutron capture on Gd has a much larger cross section than protons, carbons, and oxygens.In addition, it emits gamma rays of 8 MeV in total.Therefore, Gd-loaded water Cherekov detectors can detect neutron capture with high efficiency.These detectors are expected to improve the results of various physics analyses by using the measured neutron multiplicity to enhance flavor identification or signal-tobackground ratio.The following three physics targets are discussed: CP-phase measurements, diffused supernova neutrino background (DSNB) searches, and nucleon decay searches.In the CP-phase measurements, the charged-current quasielastic (CCQE) interaction (ν µ + n → µ − + p and νµ + p → µ + + n) is used as a signal.A comparison of neutrino and antineutrino oscillation is achieved by using accelerator neutrino to switch between neutrino and antineutrino dominant modes.However, some wrong-sign neutrinos contaminate the beam during the beam production process.The presence or absence of neutrons in the final state is expected to reduce the contamination in event selections.In the DSNB searches, the inverse beta decay (ν e + p → e + + n) is usually used as a signal, and the dominant background is atmospheric neutrinos.The background can be drastically reduced by requiring single neutron capture in the event selection while maintaining high signal efficiency [7].In the nucleon decay search, the signal usually does not accompany a neutron in the final state to violate the baryon number.If a nucleon bound in a nucleus decays, the residual nucleus is deexcited and frequently emits neutrons in the final state.The dominant background is again atmospheric neutrinos.Therefore, the signal-tobackground ratio can be increased by requiring the absence of neutron capture in the event selection [8].These analyses entirely rely on the neutrino (or nucleon decay) Monte Carlo event generators to predict how the particles emitted by neutrino interactions (or nucleon decays) respond in detectors.Therefore, the prediction accuracy of neutron multiplicity in the generators determines how much we can improve the physics results.
Returning to residual nuclear deexcitation, the fact that this process was not considered in most of the generators that employ the impulse approximation was a critical problem.This process was not simulated with few exceptions because of the small energy of the deexcited particles.However, neutron capture detection is susceptible to deexcitation since no detection threshold is imposed in neutron energies.Therefore, the highest priority is constructing a comprehensive description and simulator of neutrino-nucleus interactions, including nuclear deexcitation.NEUT [9] considers the deexcitation process only for 16 O based on Ref. [10], mainly for gamma rays.A more precise simulator applicable to 16 O and other nuclei is necessary.More simulation studies of the deexcitation process for application to neutrino interactions and nucleon decays have been discussed extensively in recent years [11][12][13][14][15][16][17].Notable points in these studies are briefly summarized below: Deexcitation simulation studies using TALYS [18], similar to this paper, are discussed in Refs.[13,14].However, since these were intended for use with liquid scintillator detector JUNO, these dealt only with carbon target and is a closed source, making it impossible to use with other detectors.A simulation study combining a neutrino generator NuWro [19] and FSI model INCL coupled with deexcitation code ABLA [20] was made focusing on carbon target in Ref. [15].ABLA is provided as an option to describe the deexcitation process in the INCL code.Another neutrino generator GENIE [21] can also simulate the deexcitation process using ABLA in principle since INCL is already implemented into it, although a detailed study has yet to emerge.ABLA lacks predictive power for gamma rays from low-energy discrete excited states, as commented in Ref. [15], and thus may be unsuitable for liquid scintillators and Cherenkov detectors.Another study using PEANUT in FLUKA [22] was performed by liquid argon time projection chambers, ArgoNeuT [17].The MeV-scale experimental data of ArgoNeuT shows good agreement with FLUKA prediction.This result indicates that PEANUT offers an excellent deexcitation model, though more verification in the GeV region and other nuclei, such as carbon and oxygen, is desired.In addition, incorporating PEANUT into the existing neutrino generators, which has yet to be achieved, might entail technical challenges.To summarize, there are several relevant studies, but many are limited in applicable nuclei, have limited validation studies, or are closed-source.The closed sources make their integration into the existing neutrino generators impossible and lead to limited use and validation of the deexcitation simulators in neutrino experiments.
Therefore, the author developed a dedicated software, NucDeEx, to simulate nuclear deexcitation.Two notable features of NucDeEx are as follows: It can be applied to 12 C and 16 O aiming for use in liquid scintillator and water Cherenkov detectors.Moreover, NucDeEx's code is provided as a package applicable to existing neutrino and nucleon decay Monte Carlo event generators.The former is because observing particles from deexcitation becomes a problem with these detectors.The latter is essential to be used in the running and future experiments, a crucial difference from similar studies that are not open source.This paper is organized as follows: Sec.II describes the overall simulation strategy; Sec.III then presents its validation using other prediction and non-neutrino beam experiments; Sec.IV finally shows the impacts of the deexcitation process on neutron multiplicity associated with neutrino interactions.In this section, three widely used neutrino Monte Carlo event generators are discussed: NEUT [9], NuWro [19], and GENIE [21].The NucDeEx codes and sample codes used with the generators are available at GitHub [23].

II. SIMULATION OVERVIEW
Neutrino Monte Carlo event generators are event-byevent simulators.The deexcitation simulator NucDeEx presented in this paper is also built as an event-by-event simulator for use with it.The basic concepts of a predecessor deexcitation simulator to NucDeEx were explained in Refs.[2,24].The predecessor was composed of two software: TALYS [18] and Geant4 [25].TALYS, an opensource package for nuclear reaction, calculates comprehensive branching ratios based on the Hauser-Feshbach model [26].Geant4, a widely used particle simulation package, calculates event-by-event kinematics using the branching ratios from TALYS.The predecessor was used with NuWro in atmospheric neutrino analysis at Kam-LAND [2].The analysis result showed excellent agreement in neutron multiplicity between predictions and observations.Note that the predecessor is not open-source and integrated into NuWro.
In order to make NucDeEx readily available for other experiments, it must be inter-operative with the existing neutrino Monte Carlo event generators.The predecessor requires Geant4, but its large library set may introduce complexity and compatibility issues.Aiming to minimize dependencies for implementation to the existing generators, the author modified the code to compute kinematics using ROOT [27] instead of Geant4.Since the generators already depend on ROOT, the fact that the deexcitation simulator depends on ROOT will not be a problem.The NucDeEx code has been opened with various other modifications, but the idea is similar to previous studies explained in Refs.[2,24].

A. Simulation procedure
Figure 2 shows a procedure of the deexcitation simulation.The information of the residual nucleus after the FSI is extracted from an event sample of neutrino interactions or nucleon decays.Particle kinematics is determined by the kinematics simulator using ROOT.The branching ratios and separation energies necessary to decide the kinematics are calculated with TALYS and tabulated in advance.Therefore, NucDeEx depends on TALYS regarding simulation outputs but is independent regarding software codes.The branching ratios and kinematics simulation details are described in Sec.II C and Sec.II D, respectively.
The key point of this design is that TALYS only calculates branching rations for nuclei, not simulating particle kinematics event-by-event and not giving excitation energy distribution.That is why the author prepares the particle simulation using Geant4 or ROOT.We need to know nuclear species, hole state, and excitation energy.These parameters are extracted from the outputs of neutrino or nucleon decay Monte Carlo event generators.
B. Branching ratios for p 3/2 -hole states Due to its low excitation energy, the deexcitation process from p 3/2 -hole states is simple.The branching ratios for p 3/2 -hole states were measured or predicted well in Refs.[11,28].NucDeEx refers to these data instead of TALYS.
In the simple shell model picture of 12 C, two nucleons lie in the s 1/2 shell, and four nucleons lie in the p 3/2 shell.No nucleon exists in the p 1/2 shell, which lies a few MeV above the p 3/2 shell.In this case, the p 3/2 -hole state always goes to the ground state (g.s.); no deexcitation happens.However, according to more precise shell model calculations, the p 1/2 shell is partially filled with nucleons From the event sample of neutrino interactions or nucleon decays, information on the residual nucleus is extracted and passed to the deexcitation simulator, NucDeEx.Using the branching ratios calculated with TALYS [18] (denoted as the orange part), particle kinematics is simulated (denoted as the blue part).Combining particle information, we get a comprehensive event sample, which includes the deexcitation process.
at a certain probability due to nucleon-nucleon correlation.As a result of this paring effect, about 20 ± 5% of the p 3/2 -hole states are expected to have excitation energies of a few MeV [12].Table .I shows a summary of the branching ratios for p 3/2 -hole states of 11 B * and 11 C * used in NucDeEx.The experimental data for 11 B * shown in this table agrees well with the predicted value mentioned above.Since no experimental data for 11 C * could be found, an analogy was assumed: While the same branching ratio is used, the excitation states are changed to ones that have the same spin-parity and similar excitation energy.The shell structures of neutrons and protons are identical except for the Coulomb potential of about 2-3 MeV.Especially in the low excitation energy region, there are excited states with the spin-parity at close energies for 11 B * and 11 C * .This fact suggests, to some extent, the validity of the analogy of assuming isospin symmetry.In reality, the symmetry is broken by the Coulomb potential, and the excitation energies are slightly different.Therefore, experiments on 11 C * to evaluate the effect of the symmetry breaking are desired.
As for 16 O, all shells including p 1/2 are filled with nucleons.The p 1/2 -hole states always go to the g.s. and p 3/2 -hole states go to excited states.Table.II shows summaries of the branching ratios for p 3/2 -hole states of 15 N * and 15 O * used in NucDeEx.It sometimes transitions to the excited states that emit a proton or multiple gamma rays.As well as 11 C * , no data for 15 O * could be found, so the same analogy is assumed.
From the p 3/2 -hole states, neither oxygen nor carbon is accompanied by a deexcited neutron.However, gamma rays affect the visible energy in neutrino detectors.The contributions from the gamma ray become relatively large and require prediction accuracy, especially in DSNB searches looking at low energy regions below 30 MeV [29] and proton decay searches [8].The probability of occurrence of p 3/2 -hole state is about twice larger than that of s 1/2 -hole states, which is discussed in Sec.II C. NucDeEx precisely considers the deexcitation from p 3/2 -hole states, although they do not affect neutron multiplicity.In the case of high excitation energies, such as s 1/2 -hole and multi-nucleon hole states, the deexcitation process becomes complicated and various particles are emitted.Branching ratios for γ, α, n, p, deuteron (d), triton (t), and 3 He are calculated with TALYS version 1.96.This calculation corresponds to the orange part in Fig. 2. A deexcitation simulation study using SMOKER code [30] was made in Ref. [12].SMOKER is another nuclear simulator that uses the Hauser-Feshbach model like TALYS, but it does not deal with deuteron, triton, and 3 He emissions.One of the features of TALYS is its ability to describe deexcitation, even for heavy ions.TALYS provides a global optical model potential to determine the transmission coefficients and several sophisticated level density models [31].The author configured TALYS to calculate the optical model for any compound nucleus.The default level density model is a combination of the constant temperature and a Fermi gas, which was adopted in previous studies [2,24].The author found that another model, the back-shifted Fermi gas model [32] gave better agreement with experimental data shown in Sec.III; therefore, the model is adopted in NucDeEx.
The calculated branching ratios of s 1/2 -hole states from 12 C and 16 O are shown in Fig. 3.The spin parity J π = 1/2 + is set for application to these states.At a typical excitation energy of 23 MeV (30 MeV) for 11 B * ( 15 N * ), neutron emission dominates about 50%.The results indicate that the deexcitation process has a nonnegligible impact on neutron multiplicity.Comparing the results of 11 B * and 11 C * , related to the 12 C target, there is an approximate isospin-symmetry in the branching ratios: For example, the branching ratio of n for 11 B * is similar to one of p for 11 C * .On the contrary, the symmetry is almost broken for 15 N * and 15 O * due to their larger Coulomb potential.A general trend, the branching ratios of p for neutron-hole nuclei are larger than ones of n for proton-hole nuclei, is seen through the results.In the simulation process, the emitted particles are first determined using the branching ratio shown in Fig. 3.The particle's energy depends on the daughter nucleus's excitation energy to be transited.This simulation step is determined as explained below.
When nuclei have high excitation energy, they rarely go to the ground state by a single-step decay but rather transition to an excited state and decay sequentially.Thus, we have to calculate the transition probabilities of all possible excited states of the daughter nuclei and the branching ratios from the states.The complete sets of information are also calculated and tabulated using TALYS.Figure 4 shows an example of the relative branching ratio of each particle as a function of the excitation energy of the corresponding daughter nucleus.After determining the emitted particles in the previous simulation step, the excited state of the daughter nucleus to be transited is determined using the relative branching ratio shown in Fig. 4. The relative branching ratios can also be interpreted as probability densities.The energies of the emitted particles are determined by this simulation step.
NucDeEx repeats these steps until the excitation energy of the daughter nucleus becomes zero.The branching ratios for all possible daughter nuclei (e.g., 10 B * produced via 11 B * → n + 10 B * ) required for these simulation steps are also calculated and tabulated in the same way.The details are omitted in this paper due to the large data size.
So far, we have discussed single nucleon hole states where we know the spin-parity exactly.However, multinucleon holes are produced by FSI with a probability of about 30% in GeV-order neutrino reactions.In the case of such multi-nucleon hole states, most neutrino Monte Carlo generators cannot determine the spin-parity of the residual nucleus because the Fermi gas model used in FSI does not consider shell levels.Therefore, branching ratios of multi-nucleon hole states are calculated using TALYS without specifying the spin-parity.The spin-parity is determined according to the level density model employed in TALYS, the back-shifted Fermi gas model.

D. Kinematics simulation
The kinematics simulation, corresponding to the blue part in Fig. 2, depends on ROOT libraries.It considers the separation energies and the four-dimensional momentum conservation.The deexcited particles are assumed to be emitted isotropically in the center-of-mass frame, namely the rest frame of the mother nucleus.The deexcitation decay process could be boosted because the residual nuclei tend to have about 200 MeV/c in momentum.The Lorentz boost from the center-of-mass frame to the laboratory frame is also considered in this simulation.The effect of the Lorentz boost is not so large but introduced to accurately reproduce the detection thresholds' effect when compared with experimental data in Sec.III.

III. RESULTS OF s 1/2 -HOLE STATE SIMULATION
Before showing the results of NucDeEx with neutrino interactions, its performance is discussed in this section.In order to validate the performance independently of neutrino interactions and generators, the deexcitation process from specific nuclear states is compared with ex-  15 N * reactions using proton beam at an energy of 392 MeV [10,33,34].They measured single-step and multistep decays of n (only for 15 N * ), p, d, t and α.The branching ratios of 3 He are negligibly small for these nuclei.Since both experiments have high incident energies, the impulse approximation is applicable, i.e., the types of incident particles do not affect the deexcitation process.Note that the single-step (multistep) decays in this paper are identical to the twobody decay (three-body decay and sequential decay) in Refs.[33,34].
The excitation energy must be specified instead of the neutrino Monte Carlo event generators in the comparisons.It is determined according to the spectral function (SF) by Benhar et al. [35], which is commonly used in generators, shown in Fig. 5.The SF provides probability distribution as a function of removal energy.The excitation energy is obtained by subtracting the separation energy.As shown in Fig. 5, this subtraction causes nonphysical negative excitation energies.We need more high-precision SF data considering discrete excited states to overcome this issue.There are two or three peaks in the excitation energy distribution, corresponding to s 1/2 -, p 3/2 -, and p 1/2 -hole states.In the following compar-isons, excitation energies above 16 or 20 MeV and below 35 or 40 MeV corresponding to the s 1/2 -hole states are selected.

A. Deexcitation from 11 B *
A comparison with the 12 C beam experiment of relative branching ratios for 11 B * with 16-35 MeV excitation energy range is shown in Fig. 6.All predictions shown here are based on TALYS.The differences come from the excitation energy distribution and level density model adopted in TALYS.The excitation energy distributions assumed in these predictions are different and slightly change the branching ratios: This work and results from Hu et al. used the SF, while the previous study (the green histogram) assumed the Lorentzian distribution with 23 MeV as the mean and 14 MeV as the FWHM.The level density model is the leading cause of the differences.As mentioned in Sec.II C, NucDeEx uses the back-shifted Fermi gas model, but the others use another level density model: A combination of the constant temperature and a Fermi gas model.The result of this work gives the best agreement with the experiment with uncertainty of about 15%.Since this experiment used 12 C beam, the center-of-mass frame was boosted, and particles emitted by deexcitation could be measured at very low thresholds in the laboratory frame.The agreement between the experimental data and this work indicates that NucDeEx predicts the branching ratios to low energies well.-hole: 19.4%  The magenta histograms show the results using NucDeEx, the green histograms show the previous study from [2], and the orange histograms show the predicted results by Hu et al. using TALYS [14].The experimental data by Panin et al. [28] is shown in the blue-hatched histograms.No experimental uncertainty is available.The branching ratios shown here account for only single-step decays.
Another comparison with the proton beam experiment of the absolute branching ratios with the same excitation energy range is shown in Fig. 7.We can compare branching ratios of both single-step and multistep decays from this data.The prediction accuracy of these decays depends not only on the branching ratios but also on the transition probabilities to the excited state, as shown in Fig. 4. A large discrepancy in t branching ratio is visible.The back-shifted Fermi gas model gives better agreement but does not solve the issue entirely; a notable difference remains.The authors of the experiment also discussed this issue, but their simulation (the blue histograms shown in Fig. 7) did not reproduce the data either.This discrepancy of t branching ratio is still an open question to be addressed in the future by performing validation experiments.The experiment applied the detection thresholds of 3.1-4.6MeV depending on particles [33], and these effects are considered in this simulation.According to NucDeEx prediction, these thresholds lead to considerable inefficiency, about 50%.To reduce this inefficiency and to better understand deexcited particles, the ion beam experiment such as Ref. [28] is preferable.The magenta histograms show the results using NucDeEx, the green histograms show the previous study from [2], and the orange histograms show the predicted results from Hu et al.
using TALYS [14].The black histograms denote experimental data from Yosoi et al. with statistical errors, and the authors provide the prediction using CASCADE code written with the blue histograms [33].The hatched or filled histograms represent the branching ratios for single-step decays, and the open histograms represent those for multistep decays.

B. Deexcitation from 15 N *
A comparison with the proton beam experiment of the branching ratios for 15 N * with 20-40 MeV excitation en-ergy range is shown in Fig. 8.Note that this experimental data was taken using the same beamline shown in Fig. 7, and the detection thresholds are also considered in this simulation.Owing to the efforts of installing neutron detectors, the experimental data of n branching ratios are available.The result of this work agrees with the experimental data in n branching ratio both of single-step and multistep decays within about 20%.Together with the results of n branching ratios for 11 B * , NucDeEx predictions show reasonable agreements with the experimental datasets.As well as 11 B * , there are large discrepancies for other particles, especially t.Since all predictions do not agree with the experimental data well and these experiments share the same beam line and detectors, validation experiments are desired.Conducting a new experiment using 16 O beam is also useful to measure particles with low energy thresholds, as demonstrated by Panin et al. introduced in Sec.III A.
Kobayashi et al. made an analysis focusing on the gamma ray from 15 N * using the same experimental data shown in Fig. 8 [10].The branching ratios for gamma rays above 3 MeV, which are detectable in water Cherenkov detectors, were investigated.Since multiple gamma rays can be emitted from these states, this measurement was of the total energy released by gamma rays E γ,tot .Table III compares gamma ray branching ratios between the experiment and a prediction using NucDeEx.The prediction agrees with the experiment for low-energy gamma rays below 6 MeV, but not above 6 MeV.The cause of this discrepancy is still under study because gamma rays above several MeV are important in searching for DSNB and nucleon decay using Super-Kamiokande.

IV. APPLICATION TO NEUTRINO INTERACTIONS
In this section, NucDeEx is evaluated with the existing neutrino Monte Carlo event generators, and its effect on neutron multiplicity is discussed.Three commonly used generators are used in this paper: NEUT version 5.6.3,NuWro version 21.09.02, and GENIE version 3.04.00with G18 10b 02 11a model configuration [36].To obtain accurate excitation energy distributions to be input to NucDeEx, the SF by Benhar et al. is adopted as the nuclear model in NEUT and NuWro.For these two generators, the excitation energy is calculated with the energies of incoming and outgoing particles obtained from the event sample.Since the SF has broad peaks and sometimes gives nonphysical negative excitation energies as shown in Fig. 5, exact energy conservation is not considered for discrete excited states, p 3/2 -and p 1/2hole states.We need more high-precision SF considering such discrete excited states to solve this problem.GE-NIE provides several sets of model configurations.This paper uses G18 10b 02 11a, which adopts the local Fermi gas and hN2018 FSI models.GENIE also offers an effective SF model as an alternative nuclear mode.However, it is not used because it did not show the peaks seen in Fig. 5.For GENIE, the excitation energy is randomly determined according to the SF by Benhar et al. irrespective of the kinematics of the event sample.CCQE interactions of ν µ and νµ with monochromatic energies of 1 GeV are investigated.The detailed treatment of excitation energy and hole states needed to connect NucDeEx and generators are explained in Sec.IV A. Then, the nominal outputs of the generators are described in Sec.IV B. At last, the results of applying NucDeEx are shown in Sec.IV C.

A. Excitation energies and hole states
In the case of single nucleon holes, the generators do not clarify the hole states while determining excitation energies.Table IV shows the excitation energy range and probability of each hole state.The hole states are determined only from the excitation energy in this application.This simple method is not ideal because the excitation energy distribution of each hole state has a finite width and is superposed, as shown in Fig. 5.A probability distribution for each hole state is needed to describe more accurately.NEUT assumes that contributions of nucleon-nucleon correlation with excitation energies off the peak of each shell as s 1/2 -hole states.It results in a larger probability of s 1/2 -hole states shown in Table IV.[35] shown in Fig. 5.The probability used in NEUT's original deexcitation process for 16 O is also listed as a Ref. [16].The excitation energy distribution in the multi-nucleon holes produced by FSI is poorly understood.In most generators, FSI considers separation energy but not excitation energy because it is based on the Fermi gas model, which does not consider shell levels.Hence, the generators cannot predict the excitation energy distribution in the multi-nucleon hole states.In applying NucDeEx to the generators, the effect of FSI is neglected, only considering the excitation energy produced by the target nucleon, even in the case of multi-nucleon hole states.It is worth mentioning that there are several studies that handle the FSI effect on excitation energy with rough approximation [14,15].This effect is from complex many-body systems of nuclei and is quite challenging to find a reasonable description.Exploring approaches like these studies and discussing more to reach a consensus is one of the essential issues to be addressed in the future.Note that the excitation energy and hole state are external input parameters to NucDeEx, and any improvements in the description will require modifications on the part of the generators.

B. Nominal output of generators
Neutron multiplicities of the nominal output of generators are shown in Fig. 9. NEUT implements the original deexcitation process only for 16 O based on Ref. [10], but it is disabled in these comparisons.The mean neutron multiplicities are shown in Table V.These results show significant variations in neutron multiplicity of about 20%, caused by different FSI implementations.This means that the FSI models need to be carefully studied to predict neutron multiplicity more precisely.NuWro (GENIE) predicts the smallest (largest) mean neutron multiplicity.This trend is consistent with nucleon transparency discussed in Ref. [37]: NuWro predicts larger transparency while GENIE gives smaller one.

C. After applying NucDeEx
Figure 10 shows neutron multiplicities after applying NucDeEx.The mean neutron multiplicities are summarized in Table V.The neutron multiplicities are increased by 20%-30% by deexcitation, which is equal to or greater than the differences of FSI models of generators as discussed in Sec.IV B. This result shows that the deexcitation processes significantly impact neutron multiplicity.It also indicates that although the uncertainty reduction of FSI is substantial, it is even more critical to correctly treat the deexcitation process in descriptions of neutron multiplicity associated with neutrino-nucleus interactions.
In the results of 16 O target, calculations from NEUT using its original deexcitation process are also listed.In the interaction of νµ + 16 O, the NEUT's calculation based on Ref. [10] gives a large discrepancy from the other calculations using NucDeEx.The difference comes from the probability of each hole state shown in Table IV and the treatment of multiple neutron emissions.The model implemented in NEUT does not consider multiple neutron emissions and overestimates single neutron emission instead.However, multiple neutron emission is not negligible in the case of high excited states with a neutron-rich nucleus, such as 15 N * s 1/2 -hole states created by CCQE interactions of νµ + 16 O → µ + + n + 15 N * .According to NucDeEx, single (multiple) neutrons are estimated to be emitted with a probability of about 55% (15%) from such states.The difference in the treatment of multiple neutron emissions results in a large difference in the probability of a single neutron emission.It leads to a large discrepancy in the number of events in a bin with a neutron multiplicity of two in Fig. 10.Since Super-Kamiokande Gadolinium has a high neutron detection efficiency of about 70%, the data is expected to verify the difference of deexcitation models in the future.
The increase in neutron multiplicity due to deexcitation depends on the generators, i.e., the FSI models, because deexcitation occurs after FSI, meaning the deexcitation depends on the residual nuclear states produced by FSI.Even if neutrino detectors measure neutron multiplicity, separating the deexcitation and FSI contributions would be quite difficult.The neutron travel distance of these two contributions should differ because of the dissimilar neutron energies: Neutrons emitted by deexcitation have only a few MeV.However, the spread of the gamma ray and vertex resolution of detectors smear the information.Therefore, it is challenging to determine which model requires modification from the neutrino experimental data only.Our future task to be addressed is conducting non-neutrino beam experiments that can precisely measure particles around the nucleus, such as those introduced in Sec.III, and investigating models after distinguishing between FSI and deexcitation contributions.This task is also vital in validating the few deexcitation experimental data so far.As well as non-neutrino beam data, creative use of neutrino data would be helpful.The value of neutrino data can be enhanced through verification by various neutrino experiments with different detection principles or combined analysis.Promoting both approaches, non-neutrino and neutrino experiments, is expected to lead to a more precise understanding of nuclear effects in neutrino-nucleus interactions.

V. CONCLUSION AND PROSPECTS
This paper reports the construction of NucDeEx to describe the nuclear deexcitation process, which was ne-glected in most existing neutrino Monte Carlo event generators until recent studies.Notable features of NucDeEx that distinguish it from similar studies are that it is opensource, can simulate deexcitation of 12 C and 16 O, and can be easily integrated into existing neutrino or nucleon decay Monte Carlo event generators.The simulator comprises two components: A nuclear simulator TALYS and a kinematics simulator based on ROOT.The quality of the simulator is evaluated by comparisons with the other theoretical predictions and hadron beam experiments.The most critical parameters, branching ratios for neutron, agree within 20%.This reproducibility is adequate for an initial implementation.The applications of NucDeEx to neutrino interactions using NEUT, NuWro, and GENIE are also demonstrated.For the predicted neutron multiplicity, a considerable impact equal to or greater than that of FSI is pointed out.This result shows the importance of considering the deexcitation process and the necessity of modifying the existing generators.
The main future tasks for improving the accuracy are to compare NucDeEx prediction with more experiments Calculations from NEUT [9], NuWro [19], and GENIE [21,36] after applying NucDeEx are shown.The green dashed lines represent the calculations from NEUT using its original deexcitation process, which is only applicable to 16 O.Neutrino energy is monochromatic 1 GeV.NuWro and INCL FSI model coupled with a deexcitation code ABLA has been performed recently [15].INCL has also been used for hadronic interactions in Geant4 and has better reproducibility than other models [38].The discussion of the deexcitation process based on the INCL FSI model is of interest.In order to avoid dependence on a specific nuclear model, it is desirable to promote the development of deexcitation simulators for both TALYS and ABLA and to compare each other.
The author provides the NucDeEx code and sample codes for applying the NucDeEx to the generators on Ref. [23].Note that NucDeEx supports carbon and oxygen targets but not for argon currently.Aiming the use in liquid-argon time-projection chamber neutrino experiment, MicroBooNE [39] and DUNE [40], an extension to argon targets is possible in principle.However, because of the increased number of shell levels, implementation would be more complex and difficult than with carbon or oxygen.There is currently no concrete plan for the extension, but it might be done upon request.The author also plans to integrate the simulator into NEUT and nucleon decay event generators used in Super-Kamiokande.After that, NucDeEx is expected to be used in various physics studies performed in Super-Kamiokande and T2K experiment [41].Integrations for other generators have yet to be planned but could be easily done by the generator developers.NucDeEx is expected to improve the pre-diction accuracy of neutron multiplicity associated with neutrino-nucleus interactions and nucleon decays.The author expects that NucDeEx is capable of leading to improving various physics results using neutron multiplicity conducted in neutrino experiments in the future: Neutrino CP-phase measurements, DSNB searches, and nucleon decay searches.

FIG. 2 .
FIG. 2.A procedure of the deexcitation simulation.From the event sample of neutrino interactions or nucleon decays, information on the residual nucleus is extracted and passed to the deexcitation simulator, NucDeEx.Using the branching ratios calculated with TALYS[18] (denoted as the orange part), particle kinematics is simulated (denoted as the blue part).Combining particle information, we get a comprehensive event sample, which includes the deexcitation process.

FIG. 4 .
FIG.4.Relative branching ratios of 11 B * with excitation energy of 23.6 MeV and J π = 1/2 + calculated with TALYS[18].The text in each panel denotes the absolute branching ratio for each particle emission, and the colored lines represent the relative branching ratio as a function of the excitation energy of the corresponding daughter nucleus.

FIG. 6 .
FIG.6.Comparison of relative branching ratios of n and d/α for 11 B * with excitation energies of 16-35 MeV.The magenta histograms show the results using NucDeEx, the green histograms show the previous study from[2], and the orange histograms show the predicted results by Hu et al. using TALYS[14].The experimental data by Panin et al.[28] is shown in the blue-hatched histograms.No experimental uncertainty is available.The branching ratios shown here account for only single-step decays.

FIG. 7 .
FIG.7.Comparison of measured and predicted branching ratios of n, p, d, t, and α for 11 B * with 16-35 MeV excitation energy.The branching ratios of n are scaled by a factor of 1/2.The magenta histograms show the results using NucDeEx, the green histograms show the previous study from[2], and the orange histograms show the predicted results from Hu et al. using TALYS[14].The black histograms denote experimental data from Yosoi et al. with statistical errors, and the authors provide the prediction using CASCADE code written with the blue histograms[33].The hatched or filled histograms represent the branching ratios for single-step decays, and the open histograms represent those for multistep decays.

FIG. 8 .
FIG.8.Comparison of measured and predicted branching ratios of n, p, d, t, and α for 15 N * with 20-40 MeV excitation energy.The branching ratios of n are scaled by a factor of 1/2.The magenta histograms show the results using NucDeEx, The black histograms denote experimental data from Yosoi et al. with statistical errors, and the authors provide the prediction using CASCADE code written with the blue histograms[34].The hatched or filled histograms represent the branching ratios for single-step decays, and the open histograms represent those for multistep decays.

TABLE I .
[28]ted states and branching ratios for p 3/2 -hole states of 11 B * and 11 C * used in NucDeEx.The excitation energy and gamma ray energy are denoted as Ex and Eγ, respectively.The branching ratios (Br) for 11 B * are measured by Panin et al.[28].All gamma rays emitted from these excited states are single.
[18]E II.Excited states and branching ratios for p 3/2 -hole states of 15 N * and 15 O * used in NucDeEx.The excitation energy and gamma ray energy are denoted as Ex and Eγ, respectively.The branching ratios (Br) for 15 N * are from[11].Gamma rays emitted from these excited states can be multiple.The relative branching ratios for gamma rays (RBrγ) are based on TALYS[18].Only those with RBrγ > 0.01 are listed.From excited states with high excitation energies, a proton is emitted with 100%, having a energy written as Ep.

TABLE III .
[10]arison of branching ratios of γ with total energy of 3 MeV< Eγ,tot < 6 MeV and Eγ,tot > 6 MeV.The experimental data from Kobayashi et al. was measured using proton beams[10].Excitation energies of 16-40 MeV are selected in this comparison.

TABLE IV .
Probability (Prob.) and excitation energy (Ex) range of each hole state for12C and16O.The excitation energy ranges are set to separate peaks corresponding to each hole state ofBenhar SF

TABLE V .
Summary of mean neutron multiplicities of CCQE interactions on12C and16O before and after applying NucDeEx.The numbers in parentheses shown in16O denote results enabling NEUT's original deexcitation process.andusethem to investigate nuclear models.As introduced in Sec.III, only a few experiments measured the deexcitation process.The model has room for improvement, but it is challenging to investigate how it should be modified from the limited data.New experiments using12C or16O beams are worthwhile since they can detect deexcited particles with lower energy thresholds than proton beam experiments.In addition, a detailed comparison between NucDeEx and other simulators is an interesting study to be discussed in the future.It is worth mentioning again that a simulation study that combines