The NEUT Neutrino Interaction Simulation

NEUT is a neutrino-nucleus interaction simulation. It can be used to simulate interactions for neutrinos with between 100 MeV and a few TeV of energy. NEUT is also capable of simulating hadron interactions within a nucleus and is used to model nucleon decay and hadron--nucleus interactions for particle propagation in detector simulations. This article describes the range of interactions modelled and how each is implemented.


Introduction
NEUT is primarily a neutrino-nucleus scattering simulation program library and provides a complete model capable of predicting the observations for a wide range of neutrino scattering experiments. NEUT is capable of simulating neutrino-nucleon and coherent neutrino-nucleus interactions in a number of reaction channels over a neutrino energy range from 100 MeV to a few TeV. Additionally, NEUT incorporates initial-and final state nuclear effects for interactions with nuclei from boron to lead. One of the most important nuclear effects is the re-scattering of hadrons, which are produced in the primary neutrino-nucleon interaction, as they propagate out of the nuclear medium. This re-scattering can result in hadron absorption, extra hadron production or knock-out, or distortion of the nuclear-leaving particle kinematic spectra. The NEUT hadron re-scattering model has also been used to simulate low-energy pion-nucleus scattering both to tune the model to extant data [1] and to simulate pion propagation in neutrino-scattering experimental simulations. Finally, NEUT can also simulate various nucleon decay channels to support experimental searches for the process.
NEUT has a long rich history, originally developed in the 1980s as a tool to study atmospheric neutrinos and nucleon decay in the Kamiokande experiment [2], and some of the original FORTRAN77 code is still in use. NEUT continues to be predominantly developed and maintained by members of the Kamiokande series of experiments (Super-Kamiokande, T2K, Hyper-Kamiokande) and many source files contain comments messages from the numerous physicists who have contributed to the simulation over the past 35 years-including those working on the Nobel prize-winning a e-mail: hayato@suketto.icrr.u-tokyo.ac.jp arXiv:2106.15809v2 [hep-ph] 25 Aug 2021 Super-Kamiokande (SK) analysis [3]. Recent development has targeted the improvements most-needed for precise and robust analyses of SK and T2K neutrino oscillation data and neutrino cross-section measurements. At T2K energies, Charged-current quasi-elastic interactions dominate. It has become clear over the last decade that such interactions can only be precisely predicted by incorporating detailed models of relevant nuclear dynamics. For the multi-GeV samples used in SK analyses, shallow and deep inelastic scattering channels are critical for modelling the expected rate of multi-ring events seen in the detector. The transition region between resonance excitation and deep inelastic scattering has proven particularly difficult to model well. Because of the in-house nature of NEUT development and analysis usage, it is not yet open source. We do not yet have the resources to migrate to an open source model, but, access to the code and usage instructions are available upon request.

Motivating the use of an Interaction Simulation
The primary goal of atmospheric and long-baseline neutrino-scattering experiments is to study neutrino oscillation, which occurs as a function of neutrino energy and flavor. As neutrinos are neutral particles, their properties (including their energy) can only be inferred from observable secondary particles produced when they interact with matter in our particle detectors. Interaction simulations are used to predict the probability of neutrinos of a given flavor and energy to interact with a given target, the observable secondary particle spectra produced in the neutrino-nucleus interaction, and event selection efficiencies and purities. The common use of nuclear targets in neutrino scattering experiments further complicates the problem as secondary hadrons produced in the primary neutrino interaction can be absorbed or lose energy before leaving the target nucleus, obfuscating the details of the primary interaction. As a result, interaction simulations are a critical tool in the analysis of neutrino oscillation data. The neutrino energy range modeled most carefully by NEUT is between 0.1 and 10 GeV. This range is motivated by the typical energies for oscillation features expected in SK and T2K data, as seen in Fig. 1 and Fig. 2.
To test oscillation hypotheses, observable distributions that correlate strongly with the neutrino energy distribution are sampled and compared to simulated predictions to extract oscillation parameter constraints. An example of such an observable is where M N,i , M N,f , and M are the mass of the initial-state nucleon, final-state nucleon, and final-state charged lepton respectively; E , p , and θ are the energy, threemomentum, and angle of the final-state charged-lepton respectively. E Rec qe is an unbiased neutrino energy estimator for the charged-current quasi-elastic (CCQE) reaction off a free neutron, ν + n → − + p, or a free proton,ν + p → + + n, which is only a function of the kinematics of the final-state charged lepton. All operating neutrinoscattering experiments use nuclear targets, where CCQE reactions occur predominantly with bound nucleons. For such interactions, E Rec qe is biased by the potential energy associated with the nucleon binding and is smeared by the Fermi motion of the nucleon confinement 1 . For interactions that exhibit strong nuclear-effects or produce extra hadrons, here referred to as CC non-QE, E Rec qe is in general a significant underestimate of the true neutrino energy. Fig. 3 (top) shows the NEUT prediction of Oscillation probabilities for atmospheric νµ to νe (top) andνµ toνe (bottom) as a function of neutrino energy and the cosine of the zenith angle (acting as a proxy for propagation distance through the Earth) assuming a normal neutrino mass hierarchy. Matter effects in the Earth produce distortions in oscillations of νµ between two and ten GeV, which are not seen forνµ oscillations. Figures reproduced from Ref. [4] The oscillation parameters used are: ∆m 2 32 = 2.5 × 10 −3 eV 2 , sin 2 θ23 = 0.5, sin 2 θ13 = 0.0219, and δCP = 0.
the oscillated rate of charged-current neutrino-oxygen interactions that produce no observable pions, a CC0π event topology. The CCQE and CC non-QE components are shown separately to highlight how the shape of their contributions differ between the E ν and E Rec qe projections. For non-QE interactions, the shape of the oscillation is largely smeared away. The CC non-QE component of the CC0π topology can be further broken down into interactions that either do not produce secondary pions, or do produce secondary pions that are subsequently re-absorbed before leaving the struck  [5,6] (energy spectral shape) and overlaid oscillation probabilities for the νµ to νµ (top) and νµ to νe and νµ to νe (bottom) as a function of neutrino energy. The black (blue) curves correspond to minimal (maximal) values of the CP-violating phase δcp. The oscillation parameters used are: ∆m 2 32 = 2.5×10 −3 eV 2 , sin 2 (θ23) = 0.5, sin 2 (2θ13) = 0.1.
nucleus. Fig. 3 (middle and bottom) illustrates the NEUT-predicted evolution of E Rec qe for a relevant range of discrete neutrino energies. It can be seen that E Rec qe exhibits significant reconstructed energy feed down above about 1 GeV, where secondary meson production is common. For higher energies, calorimetric estimators that account for hadronic energy are more accurate. However, current neutrino interaction models are often only predictive in lepton kinematics and significant uncertainties should be assigned to the predictions of hadronic particle spectra, especially for baryons and heavy mesons.
To accurately interpret observable distributions, interaction simulations are relied upon to predict the rate and observable projections for neutrino interactions over a range of energies and for a number of different target nuclei as well as the migration of these primary interactions into observable topologies. The next section details the nuclear dynamics, primary neutrino interaction, and hadronic re-scattering physics models implemented in NEUT.
The most important proton decay channels, predicted by various grand unification theories, produce a lepton and a meson, e.g. p → e + π 0 , p → µ + π 0 , or p →νK + . Neutrinos also often produce a lepton and one or more mesons when interacting with a nucleon. There are plenty of atmospheric neutrinos, which have sufficient energy to undergo such interactions, providing a significant background for proton decay searches. The majority of neutrino interactions also produce a baryon in addition to a lepton and a meson, but these may not be visible in detectors such as SK. The typical predicted momentum of mesons from both proton decays and atmospheric neutrino interactions are below several hundred MeV/c and thus, the probabilities for them to scatter before leaving the nucleus is high. Accurately predicting neutrinoinduced backgrounds and hadronic re-scattering within the nucleus is critical for the sensitivity of proton decay searches.

Simulating an Interaction
In general, NEUT factorizes the simulation of an interaction of a neutrino with flavour, , and energy, E ν , into four discrete steps. First, a specific interaction channel is chosen randomly with probability, P = σ i is the total cross section and σ i T (E ν ) is the cross section for the specific target nuclei, T , and channel, i, where i is an integer that identifies the interaction process and is defined in Table 1 (charged current) and Table 2 (neutral current). For neutrinonucleon interaction channels, the nuclear-target cross section is usually constructed as σ i T = Zσ i p +(A−Z)σ i n , where A and Z are the nucleon number and the proton number of the target nuclei and σ i p and σ i n are the bound proton and bound neutron cross sections. For historical reasons, free protons can be added to nuclear targets to build simple molecular targets such as H 2 O and CH. Fig. 4 shows the NEUT water-target cross-section predictions separated into classes of interaction channel.
Second, the primary neutrino interaction, or hard scatter, is simulated. For the majority of channels, this step involves choosing a bound nucleon from an initial-state nuclear model, then choosing interaction kinematics according to the specific interaction model, and finally choosing any remaining particle kinematics not specified by the model. This step is performed under the impulse approximation [8], which treats the target bound nucleon and the remnant nucleus as evolving independently during and after the hard scatter. This further factorizes the simulation as, to first order, the sampling of the nuclear model does not depend on the interaction kinematics chosen.
For the coherent pion-production channels (Enum 16 and 36), the interaction occurs coherently between the neutrino and the target nucleus and as a result no bound nucleon target is chosen and this is considered the final step of the simulation. For other channels, the final state hadrons are then passed on to the third step, the nucleon and meson intra-nuclear re-scattering simulation, where hadrons can elastically scatter, exchange charge with a nucleon in the nucleus, or be produced or absorbed as they are stepped out of the nuclear medium.
Finally, for oxygen targets only, the final state nuclear remnant can be left in an excited state after the interaction and a number of nuclear de-excitations, producing low energy photons (O (1 − 10) MeV), are modeled following Ref. [9]. Careful treat-  Fig. 4. The NEUT-predicted muon neutrino-water cross sections overlaid on the T2K muon neutrino flux [6], with a typical oscillation (top), and the T2K and upward atmospheric muon neutrino fluxes [7] multiplied by the charged-current inclusive total cross section (bottom). The flux multiplied by the cross section is proportional to the expected interaction rate. Above 4 GeV, the expected number of interactions in SK arising from the T2K beam falls significantly faster than from atmospheric neutrinos. n.b. The cross sections presented in the top pane are divided by the neutrino energy, whereas in the bottom pane, they are not. This is to emphasise the saturation of the interaction channels associated with lower fourmomentum transfer at SK energies and the sharp turn-on seen over T2K flux distribution. Table 1. Neutrino charged-current scattering reactions simulated by NEUT, where p refers to a proton, n to a neutron, N to a neutron or proton, A X to an entire nucleus, and W to the invariant mass of the final state hadronic system. The νµ andνµ cross sections per nucleon for each process at 600 MeV and 10 GeV are included, where B.T. stands for Below Threshold and highlights kinematically disallowed channels for 600 MeV neutrinos. For reference, the NEUT reaction identifier enumeration is included. CCQE (1p1h) crosssections are for bound nucleon in oxygen and calculated using the spectral function model. Coherent pion production cross-sections are also for oxygen. σ νµ /10 −38 cm 2 /N Channel Name Reaction 0.6 GeV 10 GeV Enum where W > 2.0 GeV ment of the de-excitation oxygen is important for precisely simulating interactions in the sensitive SK detector.
For the majority of particles produced in the hard scatter and subsequent rescattering, NEUT stores their properties in an event vector file that can be used as input to further experiment simulation processes. The only exceptions are tau and omega particles, which are decayed during the NEUT simulation by TAUOLA [10] and a custom omega decay simulation code, respectively. The final states of such decays are written to the event vector.

The Hard Scatter
This section details the neutrino scattering physics models implemented in NEUT, roughly ordered by the degree of inelasticity of the interaction. Table 2. Neutrino neutral-current scattering reactions simulated by NEUT, where p refers to a proton, n to a neutron, N to a neutron or proton, A X to an entire nucleus, and W to the invariant mass of the final state hadronic system. The νµ andνµ cross sections per nucleon for each process at 600 MeV and 10 GeV are included, where B.T. stands for Below Threshold and highlights kinematically disallowed channels for 600 MeV neutrinos. For reference, the NEUT reaction identifier enumeration is included. NCEL (1p1h) crosssections are for bound nucleon in oxygen and calculated using the spectral function model. Coherent pion production cross-sections are also for oxygen. σ νµ /10 −38 cm 2 /N Channel Name Reaction 0.6 GeV 10 GeV Enum

Quasi-Elastic Scattering
The cross-section of neutrino-nucleon charged-current quasi-elastic (CCQE) interaction was formalized by Llewellyn Smith [11]. In some of the literature, an equivalent channel is referred to as 1p1h, for one-particle one-hole. Three different nuclear models are implemented in NEUT for simulating CCQE interactions, the relativistic global Fermi-gas (GRFG) model, the local Fermi-gas (LFG) model and the spectral function (SF) model. The NEUT implementation of the GRFG coss-section follows the prescription by Smith and Moniz [12]. The LFG implementation uses the model by Nieves et al. [13] and includes an updated removal energy treatment implemented by Bourguille et al. [14]. Simple Fermi-gas models tend to over-predict the cross section for forward going leptons, as a result this model takes into account long and short range correlation of nucleons using the random phase approximation, suppressing the cross-section for low 4-momentum transfer [15]. The NEUT SF uses the spectral function by Benhar et al. [16] and the implementation is based on the one in NuWro [17] with additional improvements by Furmanski [18].
These three nuclear models differ in their treatment of the bound nucleon momentum and removal energy distributions and whether there are correlations between them. Fig. 5 shows the projections of interactions simulated with each of the three models into missing momentum (p miss ) and missing energy (E miss ), which are observable quantities for analogous measurements of electron-nuclear scattering. For the ν + 16 O → − + p + 15 O interaction, p miss and E miss are defined as where T15 O is the reconstructed kinetic energy and M15 O the ground-state mass of the nuclear remnant, which in this case is an unstable isotope but is long-lived on the timescale of the impulse approximation. These quantities are of interest because E miss is approximately the energy lost to the nuclear response during the interaction and p miss is the momentum of the struck nucleon in the lab frame. It is clear from Fig. 5 that the different nuclear models make different predictions about their distribution and correlation, notably the 16 O shell structure is visible in the SF predictions. The SF model is tuned to exclusive electron-nuclear scattering data and is only available for a subset of target nuclei, and in NEUT is only implemented for 12 C, 16 O, and 56 Fe, the most important nuclei for SK and T2K analyses. For interactions with other nuclei, such as 40 Ar, only the GRFG or the LFG nuclear models are implemented.
Beyond the choice of nuclear models, the vector and axial-vector nucleon form factors control the strength and shape of quasi-elastic interactions. There are two vector form factors of nucleon implemented in NEUT, the simple dipole form factor and BBBA05 [19]. By default, BBBA05 is used as this form factor was developed to reproduce experimental electron-scattering data.
There are four axial nucleon form factor models implemented, the dipole form factor, BBBA07 [20], the Z-expansion model [21] and the 3-component model [24]. The 3-component fit model is inspired by the 2-component fit model [23] and was created to provide additional shape freedom by expanding the 2-component model, which quickly decays with four-momentum transfer squared, Q 2 . This model has the freedom to vary the gradient of the form factor at low-Q 2 leaving the free parameters in the 2-component model to set the shape at higher momentum transfer. The 3component model is able to be continuously varied between the shape of both the 2-component and the simple dipole model. Fig. 5 presents the shape of three of the axial form factor models with associated uncertainties derived fits to hydrogen and deuterium bubble chamber data. The dipole and 3-component model uncertainties are reproduced from Stowell [24], and the Z-expansion model was fit by Meyer et. al. [22].  For neutral current elastic (NCEL) scattering, the treatment and available model components are equivalent to CCQE, except that an LFG initial-state model is not implemented for NCEL.

Charged Current Multi-Nucleon Scattering
A number of accelerator-based neutrino oscillation experiments, K2K, MINOS, and MiniBooNE, started taking high-statistics neutrino-nuclear scattering data in the 2000s. These experiments found that the number of the observed CCQE-like (equivalent to the CC0π topology introduced earlier) events were a few tens of percent larger than predicted by the models, but with a relative deficit of very forward-going muons [25,26,27]. One of the sources of these discrepancies was thought to be coming from neutrino-nucleus interaction channels, which were not implemented in the simulations used. The most-probable candidate is now believe to be the so-called multi-nucleon interaction, of which a similar process is known to exist in electronnucleus scattering. Inclusion of this interaction into neutrino-nucleus simulations was discussed by Marteau in 1999 [28].
In NEUT, the Valencia model by Nieves et al. [13] is implemented. This model considers an interaction involving the production of two nucleons and two holes in a ground-state nuclear target (often called a 2p2h interaction in the literature). Their model includes processes involving the exchange of mesons between two nucleons and thus, sometimes, it is often referred to as a Meson Exchange Current, or MEC, model. The model is not applicable for large momentum transfer and thus the threemomentum transfer to the nucleus is limited to q 3 < 1.2 GeV/c. This model does not predict how the four-momentum is distributed between the two final-state nucleons, NEUT follows the implementation in NuWro [17]. The directions of the outgoing nucleons are selected to be uniformly distributed in the center of mass frame of the nucleons. A separation energy (sometimes ambiguously referred to as the binding energy) is subtracted from the energy transfer from the lepton system. The outgoing nucleon momenta are required to be larger than the local Fermi surface at the interaction position within the nucleus. The model for the binding energy is described in Ref. [14] and depends on the interaction position in the nucleus. For oxygen, typically between ∼50 MeV and ∼75 MeV of energy is lost to the nuclear response.
Neutral current multi-nucleon scattering is not yet implemented in NEUT.

Single Meson and Gamma Productions
Single pion production is one of the dominant neutrino interaction channels in the few-GeV energy region. Therefore, it is important to understand these interactions to study neutrino oscillations using the atmospheric or the accelerator neutrinos. The particles from the single meson productions are similar to the ones from nucleon decay, and thus, these interactions are also important in the nucleon decay searches. Single pion production in NEUT is implemented following the model by Rein and Sehgal [29]. An improved model, which takes into account the lepton mass correction, by Berger and Sehgal [30] is also implemented. Both of the models simulate these interactions in two steps. First, a neutrino excites the nucleon and produce intermediate baryon resonance state, which then decays into a single meson or gamma and baryon. The range of the invariant mass of the intermediate state, W , is between the pion-production threshold and smaller than 2 GeV/c 2 . In these two models, two kinds of form factors are implemented. The first is a simple dipole form, as used in the original publications by Rein and Sehgal [29]. The axial vector mass of the dipole type form factor was selected to be 1.21GeV/c 2 based on fits to K2K data. The second form was formalized by Graczyk and Sobczyk [31] based on the Rarita-Schwinger formalization. The free parameters of Graczyk-Sobczyk form were extracted by fits to hydrogen and deuterium bubble chamber data. The direction of the pions in the resonance rest frame, the so-called Adler frame, can be determined using the full prescription by Rein [32] for all implemented resonances. Alternatively, the Rein calculation can be only used for interactions producing an intermediate ∆(1232) while, for higher-order resonances, final-state pions are distributed uniformly in the Adler frame. For nuclear-target interactions, Pauli-blocking is considered and final-state nucleon produced in the resonance decay is required to have the momentum larger than the local Fermi surface momentum modelled by an LFG. The overall effect is small, typically less than a few percent of the events are rejected. Single Kaon, Eta and γ productions are simulated using the same framework for the single pion production. The main differences are the decay probabilities (branching ratios) of each simulated resonance to the relevant final state. The production of other mesons, such as the omega, is not simulated with the model described in this section. However, such particles can be produced during the hadronization simulation in the Deep Inelastic Scattering reactions, described in the next section.
For neutrino interactions with nucleon bound within a nucleus, produced hadrons undergo the final state interactions as described in Sec. 3.3. Additionally, formation zone effects are taken into account as described in Sec. 3.2.6.

Shallow and Deep Inelastic Scattering
Shallow and deep inelastic scattering processes are separated as shown in Table 3 to avoid double counting single-and multi-meson production channels.
Interactions producing more than one meson is simulated with a custom multipion-production model. When W is larger than 2 GeV/c 2 , all the interactions, which produced at least one meson are simulated by PYTHIA v5.72 [33], included in CERNLIB 2005. For both cases, same parton distribution functions (PDFs) are used. The default PDF is a modified version of GRV98 [34] that is based on the Bodek and Yang model [35]. The multi-pion production cross section is a function of the pion multiplicity, which itself is a function of W , that models the probability to produce more than one pion, with W < 2 GeV/c 2 region. The enforcement of multiple finalstate pions avoids overlap with the resonance single meson production channels. The mean multiplicity of charged hadrons is estimated from the results of Fermilab 15-foot hydrogen bubble chamber experiment [36] and the Big European Bubble Chamber (BEBC) experiment [37]. Two parameter sets are available. The first uses the data from [36] and predicts n ch = 0.09 + 1.83 × ln(W 2 ), for all channels. The second uses the fit results from both [36] and [37]. This time, the parameter sets obtained by Bronner [38] treat each combination of ν + p, ν + n,ν + p andν + n separately by fitting the measured average charged hadron multiplicity  6. Results of the fit of the averaged charged hadron multiplicities as a function of W using the bubble chamber data sets [36,37]. Figure reproduced from Ref. [38]. data for each channel independently. The Bronner model for the ν + p channel is reproduced in Fig. 6. The obtained parameter sets are n ch = 0.58 + 1.35 × ln(W 2 ) for ν + p, (6) n ch = 0.35 + 1.24 × ln(W 2 ) for ν + n, (7) n ch = 0.41 + 1.18 × ln(W 2 ) forν + p, (8) n ch = 0.80 + 0.94 × ln(W 2 ) forν + n.
For such interaction, we assume KNO scaling to determine the value of W [39]. The forward-backward asymmetry of the charged hadron multiplicity in the hadronic center of mass system is modelled as which was derived from BEBC data [40]. The cross-section for the events with W larger than 2 GeV/c 2 is calculated without the multiplicity factor because there is no overlap with other implemented models. The kinematics of the produced particles are determined by PYTHIA.
For neutrino interactions with a nucleon bound within a nucleus, final state hadron interactions and formation zone effects are taken into account as for single meson production.

Coherent and Diffractive Pion Productions
There are two coherent neutrino-nucleus pion production models implemented in NEUT. The default model is based on the prescription by Berger and Sehgal [41]. This is an update to the previously publication model by Rein and Sehgal [42], which is also implemented in NEUT. The new model uses improved elastic pion-carbon crosssection data and lepton-mass effects are properly taken into account. With these improvements, the model is applicable to neutrinos with energies below several GeV.
A similar process, diffractive pion production, involves coherent-like pion production, but with a single proton. NEUT implements the prescription by Rein [43].
Pions produced through these channels are not affected by the final state rescattering, which is described in the following section 3.3.1.

Formation zone
The idea of a formation zone is implemented and the production positions of hadrons in nucleus are shifted from the initial neutrino interaction position. The implemented model in NEUT is based on SKAT data [44,45]. The production points of the hadrons for those interactions are shifted using the formation length (L F Z ), where p is the momentum of the hadron and µ = 0.08(GeV/c 2 ). The actual size of the shift is determined as L F Z × (− log(rand[0, 1]), where rand[0,1] is a random number from 0 to 1. The distribution of secondary hadron production position is shifted further from the center of the nucleus, the produced position of hadrons shifts to the outer region of the nucleus and thus reduce the probability of FSI.

Final-state Hadronic Re-scattering
NEUT simulates the interactions of pions, kaons, etas, omegas, protons and neutrons, produced via neutrino interactions or nucleon decay, within the nucleus. In order to simulate these hadron interactions, a custom semi-classical intranuclear cascade (INC) model is used, as in most other neutrino-nucleus simulations. In NEUT, a hadron produced in the nucleus is tracked step by step from the production point until the particle escapes from the nucleus. The size of each step is fixed at 0.2 fm. At each step, it is decided whether the particle has interacted or not using the mean free paths for the modelled interaction channels. A Woods-Saxon nucleon density function and the local Fermi-gas model are used to determine the interaction positions and kinematics of the initial and final states.

Meson interactions in nucleus
Among the modelled mesons, pion is the most important in the analyses of SK and T2K. The mean free paths for the pion interactions in the Delta region are calculated following the prescriptions by Salcedo et al. [46]. Their model takes into account the in-medium correction of ∆ self energy and uses the local Fermi-gas model. Therefore, the obtained mean free paths are expressed as functions of nuclear density and pion momentum. This model is applicable for the pions with momenta smaller than 500 MeV/c. The mean free paths for pions with momenta larger 500 MeV/c was  [47,48,49,50,51] and ABS (red) and CX (black). This figure is taken from Ref. [1].
motivated by fits to pion-nucleon scattering data and are expressed as functions of pion momentum only.
In total six types of pion scattering channels are defined: low momentum pion quasi-elastic scattering, low momentum pion charge exchange interaction, pion absorption, high momentum pion quasi-elastic scattering, high momentum pion charge exchange interaction, high momentum pion inelastic scattering (pion production). Energy-independent normalization factors were defined for each of these six channels and fit to pion-nucleus scattering data from various experiments by Pinzon Guerra et al. [1], covering pion energies up to 2GeV. Fig. 7 shows the comparison between the pion scattering data and the simulated results from NEUT and the other simulation. The kinematics of particles after the interaction are determined using the results of the phase shift analysis [52] with medium corrections as suggested by Seki et al. [53].
The interactions of kaons, etas and omegas are treated similarly as high momentum pions. The mean free paths are provided as functions of momentum for each meson and the local Fermi-gas model is used to determine the kinematics when an interaction occurs.
For kaons, (quasi-)elastic scattering and charge exchange interactions are simulated. The mean free paths of these kaon interactions are extracted from the results of the phase shift analysis of kaon scattering data [54,55]. These results are also used to determine the kinematics of outgoing particles.
For etas, the ηN → ηN, ηN → πN , and ηN → ππN processes are implemented. These interactions are simulated assuming that an eta produces an excited nucleon state, which then decays to give final state. The N (1540), N (1650) inter-mediate baryon resonances are considered. The production cross-section of baryon resonances are given by the Breit-Wigner formula: where W and M * N are the invariant mass of the intermediate baryon resonance and the mean mass of the resonance respectively, Γ tot is the total width of N * , Γ ηN is the partial width of N * → ηN , Γ X is the partial width to the final state X, where X is the final state meson, η, π or ππ. The produced particles are ejected isotropically in the intermediate resonance rest frame.

Nucleon interactions in nucleus
The implementation of nucleon scattering is also based on the INC model. Three types of interaction are considered, elastic scattering and one or two pion production. These are implemented following the work by Bertini et al.
where P surf F is the density-dependent local Fermi surface momentum. The momentum of the nucleon after an interaction is required to be larger than the local P surf F , which is the standard Pauli-blocking procedure in the local Fermi-gas model framework. Produced pions are assumed to arise from the decay of ∆(1232) resonances produced by nucleon-nucleon scattering. To simulate the ∆(1232) production a simple isobar model [58] is used. When a pion is produced in the nucleus after the nucleon rescattering, that pion is independently tracked from the point of generation using the pion transport simulation described in section 3.3.1. The predicted cross-sections for proton-carbon scattering is shown in Fig. 8.

Known Implementation Limitations
NEUT is primarily developed to enable world-leading SK and T2K neutrino oscillation, proton decay, and neutrino cross-section measurements, for which a fully consistent description of Nature is desirable but not necessary, and as such has some known limitations beyond the imperfect predictions of the implemented physics models.
The most prominent limitation is that the modelling of the initial nuclear state is tied to the individual interaction channel being simulated. For most channels, the initial state is always modelled as a Fermi-gas, for multi-nucleon scattering an LFG is used, and only for quasi-elastic scattering is the more sophisticated spectral function an option. The SF is only designed to model QE interactions, and only for specific nuclei where relevant electron-scattering data exist, for other nuclear targets a Fermigas is used. While, as noted above, the modelling of hadron re-scattering is based on the density and momentum predictions from an LFG model. model is sometimes affectionately referred to as a Franken-model, after the fictional scientist and his Gothic horror implementation. For single meson production, nuclear effects beyond the initial-state nuclear momentum and Pauli-blocking are ignored. It might be expected that the inclusion of the initial-state nucleon removal energy would affect the predictions. As statistical uncertainties are reduced with the next generation of long-baseline oscillation experiments, addressing these inconsistencies will become a focus of future development.
NEUT does not attempt to model a number of interaction channels that are irrelevant to SK and T2K analyses, these include neutrino-electron elastic scattering (used to provide complementary neutrino flux normalization constraints for higher energy beams [59]) and inverse beta decay (important for simulating reactor neutrino experiments).

Additional Tools
So far we have focused on the implementation and physics models for simulation single neutrino-nucleus interactions. While this constitutes the core of the simulation, to be an effective tool for data analysis, NEUT provides a number of other tools and features.
neutgeom: The neutgeom tool is used to interface to a neutrino beamline simulation, perform neutrino ray-tracing, interface with a detector geometry description and correctly distribute interaction positions through the detector. Currently, neutgeom is only able to consume neutrino flux vectors from JNUBEAM, which simulation the J-PARC neutrino beam for the T2K experiment. There is some interest in distributing neutgeom as a standalone tool capable of interfacing with other neutrino beam simulations and interaction simulations.
NReWeight: Cross section reweighting is an important tool for systematic error estimation for neutrino-scattering analyses. For systematic parameters that can be effectively reweighted it enables interaction model variations to be applied at analysis time, rather than requiring re-simulation and re-analysis, reducing the computing time taken to estimate systematic uncertainties by many many orders of magnitude. At its core, the reweighting technique calculates a weight for each already-simulated interaction, where x encapsulates the kinematics of the simulated interaction and p and p fully describe the model choices and any free parameters in the already-simulated and the varied model, respectively. This procedure is only exact when the varied model, p , predicts a range of x that is the same or a strict subset of the range predicted by p.
NEUT provides the NReWeight package that implements exact cross section reweighting enabling the variation of nucleon form factors for CCQE and single meson production parameters after simulation. NReWeight also exposes reweighting for the meson and nucleon intranuclear cascade, which is exact for modest variations of the underlying meson-nucleon and nucleon-nucleon scattering probabilities. NReWeight is a critical tool for SK and T2K analyses.
Interface to GEANT3 and GEANT4: The simulation of an intranuclear cascade and of meson-nuclear and nucleon-nuclear scattering is conceptually similar, the difference is whether the probe beings the simulation inside or outside of the simulated nucleus. As described in section 3.3.1, the pion-nucleon interaction cross-sections were tuned to pion-nuclear scattering data. As a result, it is attractive for physics consistency to be able to simulate pion-nucleus scattering with the NEUT INC model. To achieve this, an interface between NEUT and GEANT was developed and integrated to the T2K near and SK detector simulation programs.

Future direction
Current and future neutrino experiments require precise neutrino interaction simulations to achieve their ambitious physics goals. To meet these requirements, NEUT will be continue to be developed and improved. Current modelling improvements include the implementation of the state-of-theart CCQE and multi-nucleon model by Amaro et al. [60], and single pion production models by Kabirnezhad [61] and by Sato et al. [62].The implementation of a rudimentary QE-only electron-scattering simulation in NEUT is underway. This will enable improved extended validations of the implemented physics that is most critical to T2K analyses.
The dependence on CERNLIB is problematic, as the library is no longer maintained. Building and distributing the library for modern compilers and operating systems takes time away from more important development efforts. NEUT currently relies on CERNLIB for reading configuration files, random number generation and some common mathematical operations, and PYTHIA v5.72 for simulating SIS and DIS interactions. Resolving the dependence on PYTHIA v5 is not simple as changing to a newer version (v6 or v8) would require the implementation of the particles' kinematics determination (vector generation) functions in NEUT, as the PYTHIA implementation that we rely on has been dropped or does not cover the entire kinematic region. Removing the dependency on CERNLIB is a high priority for near-future maintenance.
The current closed source nature of NEUT is undesirable. Exposure to more users and use cases will result in code, interface, and physics improvements. However, the lack of human resources render it difficult to support NEUT as a more general tool. Work has begun, in collaboration with other neutrino interaction simulation stakeholders, to define, test, and implement a new community-designed event format and event generation API [63]. These critical future developments will be implemented in NEUT as they become defined and mature.

Summary
NEUT is a general purpose neutrino interaction simulation used and improved by members of the SK and the T2K collaborations, with critical additional contributions from interaction theory groups. Future development will target the physics requirements of Super-Kamiokande, T2K, and Hyper-Kamiokande, and software integration improvements to the NEUT API and data formats.